Performing inline 3D matrix calculations

조회 수: 1(최근 30일)
Varoujan 2015년 5월 22일
편집: Varoujan 2015년 5월 22일
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.
댓글 수: 2표시 이전 댓글 수: 1숨기기 이전 댓글 수: 1
Varoujan 2015년 5월 22일
편집: Varoujan 2015년 5월 22일
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.

댓글을 달려면 로그인하십시오.

답변(1개)

Walter Roberson 2015년 5월 22일
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.
댓글 수: 1표시 없음숨기기 없음
Varoujan 2015년 5월 22일
편집: Varoujan 2015년 5월 22일
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.

댓글을 달려면 로그인하십시오.

범주

Find more on Matrix Indexing in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by