MATLAB Answers

How to speed up this for loop containing kronecker product?

55 views (last 30 days)
DS L
DS L on 8 Apr 2021
Commented: DS L on 9 Apr 2021
I want to speed up the following for loop.
% x is n by N. S is k by N and nonnegative.
% Use random matrices for testing.
% Elapsed time of the following code is around 19 seconds.
% So, when N is large, it's very slow.
n= 800; N =100; k =10;
x = rand(n,N); S = rand(k,N); H = 0;
for i = 1: size(x,2)
X = x(:,i)*x(:,i)' ;
DW = diag( S(:,i) ) - S(:,i)*S(:,i)' ;
H = H + kron(X,DW);
end
My attempt:
kron(X, DW) = kron(x(:,i)*x(:,i)' ,diag(S(:,i))) - kron(x(:,i)*x(:,i)', S(:,i)*S(:,i)');
We can use and to rewrite the above equation.
kron(x(:,i)*x(:,i)',diag(S(:,i))) = kron(x(:,i), sqrt( diag(S(:,i))))* kron(x(:,i)*x(:,i)',diag(S(:,i)))' ; (since S is nonnegative then we can take sqrt )
kron(x(:,i)*x(:,i)', S(:,i)*S(:,i)') = kron( x(:,i), S(:,i))*kron( x(:,i), S(:,i))'.
Therefore, we only need to compute kron(x(:,i), sqrt( diag(S(:,i)))) and kron(x(:,i), S(:,i)).
Here are the codes:
n = 800; N = 100; k = 10;
x = rand(n,N); S = rand(k,N);
H1= 0; K_D= zeros(n*k, k*1, N); K_S = zeros(n*k,N);
%K_D records kron( x(:,i), sqrt (diag(S(:,i)) ) ) ,
% K_S records kron(x(:,i), S(:,i));
for i = 1:N
D_half = sqrt( diag(S(:,i)));
K_D(:,:,i) = kron( x(:,i),D_half);
K_S(:,i) = reshape (S(:,i)*x(:,i)',[],1);
end
K_D = reshape(K_D,n*k,[]);
H = K_D*K_D' - K_S*K_S';
The new codes save much time. Elapsed time is around 1 second. But I still want to speed up it.
Can someone help me speed up the above code (my attempt)? Or do someone have a new idea/method to speed up my original problem?
Thank you very much!

Answers (2)

Jan
Jan on 8 Apr 2021
I get a measureable speedup with reducing the number of sqrt() calls:
% Replace
D_half = sqrt( diag(S(:,i)));
% by
D_half = diag(sqrt(S(:, i)));
  5 Comments
DS L
DS L on 9 Apr 2021
Yes. I use this line
H = K_D * K_D' - K_S * K_S';
to replace the following sum. The sum is much slower.
for i = 1:N
DW = diag(S(:,i)) - S(:,i)*S(:,i)';
H = H + kron(x(:,i)*x(:,i)',DW);
end

Sign in to comment.


Bruno Luong
Bruno Luong on 9 Apr 2021
Edited: Bruno Luong on 9 Apr 2021
Yesterday I profile your code and the majority of time is eaten by
H = K_D*K_D' - K_S*K_S';
not by Kronecker.
The question is what is your intention of using H? If it can reduce to (matrix-vector) product (including some EIG,SVD algorithms), then you should leave the low rank decomposition (K_D/K_S), rather than storing explicitly H.
Just like using BFGS formula, no one stores the matrix B, they rather store a sequence of vectors (y,s).
Similarly, most of the time Kroneker matrix can be avoid to group the state vectors, rather using linear operators.
  5 Comments
DS L
DS L on 9 Apr 2021
Thank you. The Hessian is nonsingular because H is not the "final Hessian" and the approximation is SPD.
Because I want to try something new.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!