From: Raul on
Hi all,

Is anyone in the forum aware of any issues concerning the loss of precision when computing a matrix with a loop, as opposed to computing it with vectorized code?

I vectorized a nested loop in my original code (Code 1 below) to make things run faster, but when I do the vectorization, it turns out that some entries are the same, but some others are slightly off.

Here I attach two pieces of code, the first one done with 3 loops, and the second one with 2 loops and vectorized code for variable aa. Things like valFuncSigma, M_0v, M2_S and so on used to calculate the matrix M are not functions, but matrices that have previously been computed.

Thanks in advance!
--Raul

%%% Code 1 %%%%%%
N_S = 97;
N_v = 17;
M = zeros (N_S*N_v,N_S*N_v);
bb = 1:N_v;

for jj = 1:N_S
SS = (jj-1)/2^JS;
for kk = 1:N_S
for aa = 1:N_v
vv = (aa-1)/2^Jv;
M((jj-1)*N_v+aa,(kk-1)*N_v+bb) = 0.5*valFuncSigma(aa)*SS^2*M_2S(jj,kk)*M_0v(aa,bb) ...
+ rho*sqrt(abs(valFuncSigma(aa)*valFuncVoVol(aa)))* SS* M_1S(jj,kk)*M_1v(aa,bb) ...
+ (r-div)*SS*M_1S(jj,kk)*M_0v(aa,bb) ...
+ 0.5*funcVoVol(vv)* M_0S(jj,kk)*M_2v(aa,bb) ...
+ valFuncVolDrift(aa)*M_0S(jj,kk)*M_1v(aa,bb) ...
- r * M_0S(jj,kk)*M_0v(aa,bb);
end
end
end

%%% Code 2 %%%%%
N_S = 2^JS*Smax+1;
N_v = 2^Jv*vmax +1;
M4 = zeros (N_S*N_v,N_S*N_v);
bb = 1:N_v;
aa = 1:N_v;
vv = (aa-1)/2^Jv;
tic
for jj = 1:N_S
SS = (jj-1)/2^JS;
for kk = 1:N_S

M((jj-1)*N_v+aa,(kk-1)*N_v+bb) = 0.5*repmat(valFuncSigma(aa),1,N_v)'.*(SS^2*M_2S(jj,kk)*M_0v(aa,bb)) ...
+ rho*sqrt(abs(repmat(valFuncSigma(aa),1,N_v)'.*repmat(valFuncVoVol(aa),1,N_v)')).* (SS* M_1S(jj,kk)*M_1v(aa,bb)) ...
+ (r-div)*SS*M_1S(jj,kk)*M_0v(aa,bb) ...
+ 0.5*(repmat(funcVoVol(vv),N_v,1))'.*(M_0S(jj,kk)*M_2v(aa,bb)) ...
+ repmat(valFuncVolDrift(aa),1,N_v)'.*(M_0S(jj,kk)*M_1v(aa,bb)) ...
- r * M_0S(jj,kk)*M_0v(aa,bb);

end
end