Prev: Blockwise Matrix Expansion
Next: help with matlab code
From: Irl on 19 Jul 2010 13:49
MATLAB seems to implement x^2 as a multiply but x^3, or higher power, via something much slower (exp of number times log?). For large enough n, exp n log will be faster than multiplying n times, but the internal code generator seems to err wildly on the side of using the slow method, whatever it is, rather than just repeated multiplies. But it does do x^2, apparently, as x*x. Is there a way to make it use repeated multiplies for higher exponents? I wrote a Zernike polynomial program and it's very tiresome to go change all the powers to multiply strings.
Example program, including comments giving the execution times:
function TestIntOpt(NumEl)
if nargin<1; NumEl = 1e7; end
% Commented results are for 1e7 as input argument
tt=rand(1,NumEl);
tic
x = tt; %zero ms
T1=toc;fprintf('Simple assign, %.0f elements: %.0f ms\n',NumEl,1000*T1);
tic
x = tt.^2; %takes 70 msec
T1=toc;fprintf('Second power (square): %.0f ms\n',1000*T1);
tic
x = tt.^3; %takes 1200 msec
T1=toc;fprintf('Third power (cube): %.0f ms\n',1000*T1);
tic
x = tt.^2.*tt; %takes 100 msec
T1=toc;fprintf('Cube via square times first power: %.0f ms\n',1000*T1);
tic
x = tt.^10; %takes 1200 msec
T1=toc;fprintf('Tenth power: %.0f ms\n',1000*T1);
tic
x = tt.*tt.*tt.*tt.*tt.*tt.*tt.*tt.*tt.*tt; %takes 170 msec
T1=toc;fprintf('Tenth power via nine multiplies: %.0f ms\n',1000*T1);
tic
x = exp(10*log(tt)); %takes 660 msec
T1=toc;fprintf('Tenth power via exp 10 log: %.0f ms\n',1000*T1);
From: Walter Roberson on 19 Jul 2010 14:25
Irl wrote:
> MATLAB seems to implement x^2 as a multiply but x^3, or higher power,
> via something much slower (exp of number times log?). For large enough
> n, exp n log will be faster than multiplying n times, but the internal
> code generator seems to err wildly on the side of using the slow method,
> whatever it is, rather than just repeated multiplies. But it does do
> x^2, apparently, as x*x. Is there a way to make it use repeated
> multiplies for higher exponents?

What you say about the speed agrees with my tests of 2007a. When I noticed
this, I implemented a routine that tested the power and some conditions on the
base value and used the most efficient routine (e.g., sometimes
multiplication, sometimes realpow(), falling back to .^ for cases uncommon
enough not to be worth special handling in my usage.)

However, after a while of using that, I encountered some subtle accuracy
issues that took a while to track down (partly because I "knew" the "right
answer"). In time I abandoned this approach and went back to .^ when I realized:

For floating point values, repeated multiplication to form x^n is not the same
as exp(n*log(x)) due to differences in the propagation of errors (with
repeated multiplication being less accurate.)
From: James Tursa on 19 Jul 2010 16:24
Walter Roberson <roberson(a)hushmail.com> wrote in message <i225i6$anh$1(a)canopus.cc.umanitoba.ca>...
> Irl wrote:
> > MATLAB seems to implement x^2 as a multiply but x^3, or higher power,
> > via something much slower (exp of number times log?). For large enough
> > n, exp n log will be faster than multiplying n times, but the internal
> > code generator seems to err wildly on the side of using the slow method,
> > whatever it is, rather than just repeated multiplies. But it does do
> > x^2, apparently, as x*x. Is there a way to make it use repeated
> > multiplies for higher exponents?
>
> What you say about the speed agrees with my tests of 2007a. When I noticed
> this, I implemented a routine that tested the power and some conditions on the
> base value and used the most efficient routine (e.g., sometimes
> multiplication, sometimes realpow(), falling back to .^ for cases uncommon
> enough not to be worth special handling in my usage.)
>
> However, after a while of using that, I encountered some subtle accuracy
> issues that took a while to track down (partly because I "knew" the "right
> answer"). In time I abandoned this approach and went back to .^ when I realized:
>
> For floating point values, repeated multiplication to form x^n is not the same
> as exp(n*log(x)) due to differences in the propagation of errors (with
> repeated multiplication being less accurate.)

For low powers, e.g. 2 or 3 etc., this is likely not an issue. The direct multiplication will be faster and just as accurate. For higher powers, the propagation of errors could probably be controlled with a smarter multiplication scheme (e.g., binary decomposition of the power), but I don't know at what point the speed becomes a wash. I would have to run some tests. e.g., a modified version of this FEX submission to use times instead of mtimes could do it somewhat optimally:

http://www.mathworks.com/matlabcentral/fileexchange/25782-mpower2-a-faster-matrix-power-function

James Tursa
From: Bruno Luong on 19 Jul 2010 17:56
"Irl " <irl_smith(a)raytheonRemoveTheseWords.com> wrote in message <i2236g$oco$1(a)fred.mathworks.com>...
> MATLAB seems to implement x^2 as a multiply but x^3, or higher power, via something much slower (exp of number times log?). For large enough n, exp n log will be faster than multiplying n times, but the internal code generator seems to err wildly on the side of using the slow method, whatever it is, rather than just repeated multiplies. But it does do x^2, apparently, as x*x. Is there a way to make it use repeated multiplies for higher exponents? I wrote a Zernike polynomial program and it's very tiresome to go change all the powers to multiply strings.

Usually when using Zernike polynomials, not only a specific power is needed, but all powers in the range are needed. In this case, you might compute them at once (for each variable). This function can accomplish that, using the same idea of binary decomposition mentioned by James. Put the following in the same mfile, then call this example in command line:

p = allpow(1.1,100);
1.1^100 / p(end)

function powers = allpow(x, pmax)
% function powers = allpow(x, pmax)
% compute all powers of order(0,1,...pmax) of scalar x, return the
% result in powers (having length=pmax+1)

powers= zeros(1, pmax+1);
powers(1) = 0;
powers(2) = x;
for k=2:pmax
[p2 r] = bindecomp(k);
powers(k+1) = powers(p2+1)*powers(r+1);
end
end %% allpow

function [p2 r] = bindecomp(n)
% [p2 r] = bindecomp(n)
% Given interger n>1, find p2 and r such that:
% p2 is power-of-two
% p2+r=n
% p2<n and r<n
previouspower2 = floorlog2(n);
if previouspower2==n
p2 = n/2;
else
p2 = previouspower2;
end
r = n-p2;
end % bindecomp

function previouspower2 = floorlog2(n)
% function previouspower2 = floorlog2(n)
% Given n integer >= 1
% Find the larger power-of-two smaller or equal to n
previouspower2 = 1;
while n>1
n = floor(n/2);
previouspower2 = 2*previouspower2;
end
end

% Bruno
From: Steve Amphlett on 20 Jul 2010 02:46
"Irl " <irl_smith(a)raytheonRemoveTheseWords.com> wrote in message <i2236g$oco$1(a)fred.mathworks.com>...
> MATLAB seems to implement x^2 as a multiply but x^3, or higher power, via something much slower (exp of number times log?). For large enough n, exp n log will be faster than multiplying n times, but the internal code generator seems to err wildly on the side of using the slow method, whatever it is, rather than just repeated multiplies. But it does do x^2, apparently, as x*x. Is there a way to make it use repeated multiplies for higher exponents? I wrote a Zernike polynomial program and it's very tiresome to go change all the powers to multiply strings.
> Example program, including comments giving the execution times:
> function TestIntOpt(NumEl)
> if nargin<1; NumEl = 1e7; end
> % Commented results are for 1e7 as input argument
> tt=rand(1,NumEl);
> tic
> x = tt; %zero ms
> T1=toc;fprintf('Simple assign, %.0f elements: %.0f ms\n',NumEl,1000*T1);
> tic
> x = tt.^2; %takes 70 msec
> T1=toc;fprintf('Second power (square): %.0f ms\n',1000*T1);
> tic
> x = tt.^3; %takes 1200 msec
> T1=toc;fprintf('Third power (cube): %.0f ms\n',1000*T1);
> tic
> x = tt.^2.*tt; %takes 100 msec
> T1=toc;fprintf('Cube via square times first power: %.0f ms\n',1000*T1);
> tic
> x = tt.^10; %takes 1200 msec
> T1=toc;fprintf('Tenth power: %.0f ms\n',1000*T1);
> tic
> x = tt.*tt.*tt.*tt.*tt.*tt.*tt.*tt.*tt.*tt; %takes 170 msec
> T1=toc;fprintf('Tenth power via nine multiplies: %.0f ms\n',1000*T1);
> tic
> x = exp(10*log(tt)); %takes 660 msec
> T1=toc;fprintf('Tenth power via exp 10 log: %.0f ms\n',1000*T1);

Related...

If you are writing MEX files (or just normal programs/functions) that use powers, consider switching to C++ and use std::pow() in place of C's very restrictive pow() and powf() functions. It has overloads for integral exponents, which can choose the optimum method (multiple multiplications vs logarithms vs ???).

I cannot comment on what FORTRAN and its modern incarnations do for intergral exponents.

I'm guessing that Matlab _probably_ uses C's pow() and powf() internally.
 |  Next  |  Last
Pages: 1 2
Prev: Blockwise Matrix Expansion
Next: help with matlab code