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.