More about exponents

My previous post raised some interest (i.e., one person asked me about it), so I decided to follow up a little bit on that. In a nutshell, I suggested to use logarithms and exponents to raise large matrices to some non-integer power because it's faster. That is, exp(log(A)*x) is faster than A.^x, subject to the constraint that A>0 everywhere. Referring to this, I can confirm that it works on Octave as well as in Matlab. I wouldn't know about other programming languages though, if you feel like trying and letting me know in the comments, I'll create a new post about it.

To that, I wanted to add something else. If we want to raise a large matrix to some integer power, it's worth using the exponentiation by squaring. The following Matlab code is significantly faster than the .^ operator already when the exponent is 4.

function y = pow(x,p)

y = ones(size(x));

while (p > 0)
  if (mod(p,2) == 1)
    y = y.*x;
  end
 
  p = floor(p/2);
  x = x.*x;
end

Interestingly, this does not work for Octave, which probably implements a similar algorithm when it realises that the exponent is an integer. So, if you use Octave keep going with your .^ operator.

In Matlab the code doesn't look great, probably because most of its commands are native for double precision floating-point arithmetic, while here we're dividing the exponent which is an integer. Somehow, the same function written in C looks a little bit more neat.

float powN(float x, unsigned int p) {
  float y = 1.0;
 
  while (p > 0) {
    if (p & 1)
      y *= x;
  
    x *= x;
    p >>= 1;
  }
 
  return y;
}

Here we can use directly bitwise operators in the condition and in place of the division by 2. That said, this is a very low-level operation which is likely to be implemented more efficiently in whichever API you're using (except for Matlab, it seems). So I'd recommend you to double-check first, because it might very well be that you can use a vectorised version that outperforms this one.

.^ operator vs. logs in Matlab

Quick comparison test. Say you have a large matrix A, what's faster, A.^2.4 or exp(log(A)*2.4)?

Obviously, the fact that I'm asking should kind of give away the answer and the latter is faster. There's however a catch in all this, so it really depends on what you're doing, because you can't take log(0). But, for what I want to do this isn't a problem.

Let's say that our big matrix is actually an image in sRGB colour space, and I want to convert it to linear RGB, that is, I want to remove the gamma correction. From the Wikipedia page I find the formula, and the straightforward implementation is:

linear = srgb;
small = linear < 0.04045;
linear(small) = srgb(small) / 12.92;
linear(~small) = ((srgb(~small) + 0.055)/1.055) .^ 2.4;

But if the image is big (or your computer is slow) we can save seconds using logs.

small = srgb < 0.04045;
linear = exp(log((max(srgb,1e-3)+0.055)/1.055) * 2.4);
linear(small) = srgb(small) / 12.92;

This code, if in a function, is about three times as fast as the previous.

In this case, clipping the image to 0.001 works perfectly fine whether you use 8 or 16 bit precision images, because all those small values are only transformed linearly by this gamma function. Obviously, if you're not doing image processing this may not be good enough for you and you might have to stick to the .^ operator.

Anyhow, this works also for applying the gamma to a linear RGB image. The code is the following:

small = linear < 0.0031308;
srgb = exp(log(1.055) + log(max(linear,1e-3)) / 2.4) - 0.055;
srgb(small) = linear(small) * 12.92;