Performing inline 3D matrix calculations
3 次查看(过去 30 天)
显示 更早的评论
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 个评论
Walter Roberson
2015-5-22
I would still like to know the speed difference between doing this as single precision compared to doing it in double precision.
回答(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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Function Creation 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!