replacing for loop with more efficient code

1 Ansicht (letzte 30 Tage)
Lukas Müller
Lukas Müller am 10 Jan. 2019
Kommentiert: Anusha Sridharan am 11 Jan. 2019
Hi, given symmetric matrices für i=1,...,m, positiv definite and symmetric and ,
I want to compute with .
I solved this with the following for loop:
for i=1:m
M_row=X*A(:,(i-1)*n+1:i*n)*inv_S;
for j=i:m %since M is symmetric
M(i,j)=trace(M_row*A(:,(j-1)*n+1:j*n));
end
end
Is there a more efficient way to solve this in matlab? with vectorization?
  2 Kommentare
Bruno Luong
Bruno Luong am 11 Jan. 2019
OP wants to compute
M(i,j) = Trace(X * A(:,:,i) * inv(S) * A(:,:,j))
where A(:,:,i), X and S are n x n symmetric, matrices, S and X are definite positive. i=1,...m
His original code that needs to be optimized is
M = zeros(m,m);
inv_S = inv(S);
for i=1:m
M_row=X*A(:,(i-1)*n+1:i*n)*inv_S;
for j=i:m %since M is symmetric
M(i,j)=trace(M_row*A(:,(j-1)*n+1:j*n));
end
end
Anusha Sridharan
Anusha Sridharan am 11 Jan. 2019
[Answers Dev] Restored Edits

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Bruno Luong
Bruno Luong am 10 Jan. 2019
% Generate test matrices
m = 5;
n = 10;
S = rand(n,n);
S = 0.5*(S + S.');
X = rand(n,n);
X = 0.5*(X + X.');
A = zeros(n,n,m);
for k=1:m
Ak = rand(n,n);
Ak = 0.5*(Ak + Ak.');
A(:,:,k) = Ak;
end
AA = reshape(A,[n n m]);
M = zeros(m,m);
for i=1:m
Ti = S\(AA(:,:,i)*X);
M(i,:) = Ti(:).'*reshape(AA(:,:,:),[],m);
end
  1 Kommentar
Bruno Luong
Bruno Luong am 11 Jan. 2019
The loop can be replaced by vectorized code if you use this mtimesxmtimesx FEX
AA = reshape(A,[n n m]);
Ti = mtimesx(mtimesx(inv(S),AA),X);
M = reshape(Ti,[],m).' *reshape(AA(:,:,:),[],m)

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (1)

Lukas Müller
Lukas Müller am 10 Jan. 2019
Thanks for you fast answer. So I guess it doesnt work without at least one for loop. Still way better then before thanks a lot :)

Kategorien

Mehr zu Loops and Conditional Statements finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by