Prev: crating function
Next: Urgent Help Required
From: David Romero-Antequera on 5 Aug 2010 10:53
"David Romero-Antequera" <dromero_fisica(a)yahoo.es> wrote in message <i3c7i1$adm$1(a)fred.mathworks.com>...
> Hello, everyone.
>
> I need some suggestion to make the following function faster.
>
> function psi=NLF(z,y,funvector,orig_size)
> psi=zeros(size(y));
> Y=ifft2(reshape(y,orig_size));
> for fun=funvector
> fY=feval(fun{2},z,Y);
> kY=fun{1}.*fft2(fY);
> kY=kY(:);
> psi=psi+kY;
> end
> end
>
> where z is a double scalar, y is a complex matrix, funvector is a cell array with double scalars in the first column and function handles in the second, and orig_size is a 2x1 matrix.
>
> As you can see, I need to compute the inverse Fourier transform of the vector Y, then evaluate it into several user-provided functions and compute the Fourier transform and add all of that together. This is part of an ODE, so you might expect that this is going to be evaluated THOUSANDS of times, and every second counts.
>
> Any suggestions?
> Thanks in advance.

THANKS TO EVERYONE.
I must say that this discussion have been really helpful, and I reached my goal.
Hope this helps somebody else as well.
From: Walter Roberson on 5 Aug 2010 11:05
David Romero-Antequera wrote:

> I have to made some more testing, specially to see if I am not
> compromising the accuracy of the whole computation. However, up to now,
> the implementation of the modified exponential have saved me TWO THIRDS
> of the time!!!
>
> WOW!
> I mean, I know that this is probably best case scenario, but...
> WOW!!!!

Well, that sounds promising. How did you end up implementing

exp(1i*t*M)

with varying t ?
From: David Romero-Antequera on 5 Aug 2010 12:12
Walter Roberson <roberson(a)hushmail.com> wrote in message <XqA6o.1209$1v3.1198(a)newsfe20.iad>...
> David Romero-Antequera wrote:
>
> > I have to made some more testing, specially to see if I am not
> > compromising the accuracy of the whole computation. However, up to now,
> > the implementation of the modified exponential have saved me TWO THIRDS
> > of the time!!!
> >
> > WOW!
> > I mean, I know that this is probably best case scenario, but...
> > WOW!!!!
>
> Well, that sounds promising. How did you end up implementing
>
> exp(1i*t*M)
>
> with varying t ?

Well, this is a preliminary version, I still need to check accuracy and stuff. Suggestions and advice are very welcome.


function y = fexp(D,t)
persistent A t0; % Status variables
persistent fact L10 m; % Constants

if ischar(D) % Initialization
A=[]; t0=[]; fact=[]; L10=[];
else
if isempty(A) % Empty matrix: normal exponential
y=exp(1i*t*D);
fact = 1./[6.2270208e+009 479001600 39916800 3628800 ...
362880 40320 5040 720 120 24 6 2 1 1];
m = 9; % Accuracy (1e-m)
L10=abs(log(10^(-m)));
else % Fast computation
dt=t-t0;
ldt=log(abs(dt));
n=fix(L10/ldt); % # terms of expansion to achieve 1e-m
if n>m
y=exp(1i*t*D);
else
y=A.*polyval(fact(1:n+5),1i*dt*D); % Compute +4 terms
end
end
A=y;
t0=t;
end
end
From: Walter Roberson on 6 Aug 2010 00:17
David Romero-Antequera wrote:

> function y = fexp(D,t)
> persistent A t0; % Status variables
> persistent fact L10 m; % Constants
> if ischar(D) % Initialization
> A=[]; t0=[]; fact=[]; L10=[];

%At this point you have not assigned or used m, but you set up for the
%second state. You have assigned fact to be empty here but have not used
%it yet.

> else
> if isempty(A) % Empty matrix: normal exponential
%this section can only be entered once per
%re-initialization, because after the endif you always
%assign a non-empty value to A.

> y=exp(1i*t*D);
> fact = 1./[6.2270208e+009 479001600 39916800 3628800 ...
> 362880 40320 5040 720 120 24 6 2 1 1];

%this calculation of fact can only be done once per
%reinitialization, for the same reason. You do not change
%fact once it is set, so you might as well assign it
%*before* the initialization case (protected with
%isempty)) and make it persistent

> m = 9; % Accuracy (1e-m)

%likewise, this is the only assignment to m and so it might
%as well be done once at the top (protected by the
%isempty() for fact)

> L10=abs(log(10^(-m)));

%L10 is really a constant too

> else % Fast computation
> dt=t-t0;
> ldt=log(abs(dt));
> n=fix(L10/ldt); % # terms of expansion to achieve 1e-m
> if n>m
> y=exp(1i*t*D);
> else
> y=A.*polyval(fact(1:n+5),1i*dt*D); % Compute +4 terms
> end
> end
> A=y;
> t0=t;

%If I have understood the algorithm (which I might not have),
%then I would expect t0 to be updated to t only in the case
%where the prediction was found to be too imprecise and the
%the full exp(1i*t*D) is done again. I do not, though, catch
%how you are determining if the precision would be high enough
%so possibly your current version is correct here.

> end
> end
From: David Romero-Antequera on 6 Aug 2010 09:17
Walter Roberson <roberson(a)hushmail.com> wrote in message <o1M6o.2819$EF1.2017(a)newsfe14.iad>...
> David Romero-Antequera wrote:
>
> > function y = fexp(D,t)
> > persistent A t0; % Status variables
> > persistent fact L10 m; % Constants
> > if ischar(D) % Initialization
> > A=[]; t0=[]; fact=[]; L10=[];
>
> %At this point you have not assigned or used m, but you set up for the
> %second state. You have assigned fact to be empty here but have not used
> %it yet.
>
> > else
> > if isempty(A) % Empty matrix: normal exponential
> %this section can only be entered once per
> %re-initialization, because after the endif you always
> %assign a non-empty value to A.
>
> > y=exp(1i*t*D);
> > fact = 1./[6.2270208e+009 479001600 39916800 3628800 ...
> > 362880 40320 5040 720 120 24 6 2 1 1];
>
> %this calculation of fact can only be done once per
> %reinitialization, for the same reason. You do not change
> %fact once it is set, so you might as well assign it
> %*before* the initialization case (protected with
> %isempty)) and make it persistent
>
> > m = 9; % Accuracy (1e-m)
>
> %likewise, this is the only assignment to m and so it might
> %as well be done once at the top (protected by the
> %isempty() for fact)

You're totally right about these suggestions. Thanks. The thing is that in my first version, I wasn't using an initialization procedure, which is a very bad idea if you're using persistent variables (I didn't had that figured out, it is the first time with persistent variables).

>
> > L10=abs(log(10^(-m)));
>
> %L10 is really a constant too
>
> > else % Fast computation
> > dt=t-t0;
> > ldt=log(abs(dt));
> > n=fix(L10/ldt); % # terms of expansion to achieve 1e-m
> > if n>m
> > y=exp(1i*t*D);
> > else
> > y=A.*polyval(fact(1:n+5),1i*dt*D); % Compute +4 terms
> > end
> > end
> > A=y;
> > t0=t;
>
> %If I have understood the algorithm (which I might not have),
> %then I would expect t0 to be updated to t only in the case
> %where the prediction was found to be too imprecise and the
> %the full exp(1i*t*D) is done again. I do not, though, catch
> %how you are determining if the precision would be high enough
> %so possibly your current version is correct here.
>
> > end
> > end

For the exponential series, the error goes as (dt^(m+1))/factorial(m+1). Here I used m as the accuracy goal (in this case, 1e-9), and computed the number of terms needed to achieve AT LEAST this accuracy using n=fix(L10/ldt). This is by no means a good approximation of the real formula, but it is, however, always higher than the real one, which is good. On the other hand, I needed an upper number of Taylor terms, and I just decided to use m as well.

The real intention here was to optimize time, and then I will improve accuracy from there. However, it seems that this is not a good option, because Taylor series converge really slow, and I will need a lot of terms to achieve a good accuracy, which means an slower algorithm.

I'm trying some other things, but if anyone else have a better idea, I surely would like to try it. :)
 |  Next  |  Last
Pages: 1 2
Prev: crating function
Next: Urgent Help Required