55 views (last 30 days)

Show older comments

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!

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)));

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.

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

Start Hunting!