Performing inline 3D matrix calculations

All:

I have the following code:

% f, tau, jonesx, jonesy are (1, numParams) vectors
% sxyt = 1e6; numFrames = 10; numParams = 18;
A = single(zeros(sxyt, numFrames));
B = single(zeros(sxyt, numFrames));
osc = single(zeros(sxyt, numFrames, numParams));
mag = single(zeros(sxyt, numFrames, numParams));
DxDyOut = single(zeros(sxyt, numFrames, 2));
for k = 1:numParams
    osc(:,:,k) = exp((2.*pi.*f(k).*t(:,:) + randomPhase(k)).*1i);
    mag(:,:,k) = a(k) .* (k1(k) + k2(k).*exp(-t(:,:)/(tau(k) + 0.0001))); 
    A = A + real( jonesx(k) .* osc(:,:,k) .* mag(:,:,k) );
    B = B + real( jonesy(k) .* osc(:,:,k) .* mag(:,:,k) );
end
DxDyOut(:,:,1) = A;
DxDyOut(:,:,2) = B;

I would like to perform these operations inline rather than in the loop. This code gets called many times and for each time it takes ~7.5s for first equation and 2.5s for 2nd thru 4th equations (for sxyt = 1e6; numFrames = 10; numParams = 18), so roughly 400 ms and 140 ms per execution of a line.

Is there any way one can eliminate the loop? I figure it will be about 10x faster inline.

I looked at mtimes and mtimesx (in user submissions) but can't quite figure out how either may be used for this. I thought I could use mtimesx submission as follows:

• Convert the (1, numParams) vectors to (numFrames,numFrames,numParams) size - each value in the (:,:,i) would be the same as the (1,i) element.
• Convert the (sxyt, numFrames) vectors to (sxyt,numFrames, numParams) size

Any help or suggestion on how to use mtimesx or any other method is much appreciated.

[Walter Roberson]

I would still like to know the speed difference between doing this as single precision compared to doing it in double precision.

[Varoujan]

I will try to get the answer to that. However, the person I inherited the simulation code from may have done this to save memory - the osc and mag matrices are 1.5 GB and 0.8GB in single.

Answer 1

[Walter Roberson]

Pull values out of the loop when practical. For example 2.*pi.*f(k) can be calculated ahead of the loop as can k2(k)/(tau(k)+0.0001). You can also factor out a(k)

f_2pi = single(2.*pi.*f);
k2_tau = single(f2 ./ (tau + 0.0001));
a_k1 = a .* k1;

for k = 1:numParams
    osc(:,:,k) = exp(f2pi(k).*t(:,:) + randomPhase(k)).*1i);
    mag(:,:,k) = a_k1(k) + a(k) .* k2_tau(k).*exp(-t(:,:)); 
    A = A + real( jonesx(k) .* osc(:,:,k) .* mag(:,:,k) );
    B = B + real( jonesy(k) .* osc(:,:,k) .* mag(:,:,k) );
end

and then clearly a(k) can be factored into k2_tau(k):

f_2pi = single(2.*pi.*f);
a_k2_tau = single(a .* f2 ./ (tau + 0.0001));
a_k1 = single(a .* k1);

for k = 1:numParams
    osc(:,:,k) = exp(f2pi(k).*t(:,:) + randomPhase(k)).*1i);
    mag(:,:,k) = a_k1(k) + a_k2_tau(k).*exp(-t(:,:)); 
    A = A + real( jonesx(k) .* osc(:,:,k) .* mag(:,:,k) );
    B = B + real( jonesy(k) .* osc(:,:,k) .* mag(:,:,k) );
end

Also experiment,

for k = 1:numParams
    osck = single(exp(f2pi(k).*t + randomPhase(k)).*1i));
    magk = single(a_k1(k) + a_k2_tau(k).*exp(-t));
    osc_magk = osck .* magk;
    A = A + real( jonesx(k) .* osc_magk );
    B = B + real( jonesy(k) .* osc_magk );
    osc(:,:,k) = osck;
    magk(:,:,k) = magk;
end

The next question would be whether jonesx and jonesy consist entirely of real values. If they do, then you can move the real():

for k = 1:numParams
    osck = single(exp(f2pi(k).*t + randomPhase(k)).*1i));
    magk = single(a_k1(k) + a_k2_tau(k).*exp(-t));
    osc_magk = single(real(osck .* magk));
    A = A + jonesx(k) .* osc_magk;
    B = B + jonesy(k) .* osc_magk;
    osc(:,:,k) = osck;
    magk(:,:,k) = magk;
end

Possibly some of the single() calls can be omitted.

[Varoujan]

Thanks. Let me study your answer and try to implement. I'll comment back after I try it out.

Looking at it quickly though, I see that some of your suggestions aren't applicable since you mis-read the equations for osc and mag. For instance, you can't pull tau out since it's inside the exponential.

Nevertheless, your general point about pulling as much of the repeating computation out of the loop is a good one. I will do that soon and find out how much benefit I get.