How to calculate matrix entries efficently using Matlab

Question

I have a cell array myBasis of sparse matricies B_1,...,B_n.

I want to evaluate with Matlab the matrix Q(i,j) = trace (B^T_i * B_j).

Therefore, I wrote the following code:

for i=1:n
    for j=1:n
        B=myBasis{i};        
        C=myBasis{j};   
        Q(i,j)=trace(B'*C);
    end
end

Which takes already 68 seconds when n=1226 and B_i has 50 rows, and 50 colums.

Is there any chance to speed this up? Usually I exclude for-loops from my matlab code in a c++ file - but I have no experience how to handle a sparse cell array in C++.

You could try a parfor and run the loops in parallel. http://www.mathworks.se/help/distcomp/parfor.html — PureW, Jul 24 '14 at 13:16
Actually, this might be a nice exercise to implement in C++. — Shai, Jul 24 '14 at 13:48
@patrik Actually I want to run my script on 12800000000 matricies where each matrix has 400 rows, 400 columns, and at most 405 entries. Therefore I think I should use sparse matricies. — Adam, Jul 24 '14 at 13:55
Even storing 12 800 000 000 8 byte numbers will require 95 TB of storage so that is nonsense in my opinion — EJG89, Jul 24 '14 at 14:00
@patrik I just noticed that I made a calculation error instead of 12800000000 it are only 80000. Sorry about that. — Adam, Jul 24 '14 at 14:01

score 4 · Accepted Answer · edited May 23 '17 at 12:05

4

As noted by Inox Q is symmetric and therefore you only need to explicitly compute half the entries.
Computing trace( B.'*C ) is equivalent to B(:).'*C(:):
trace(B.'*C) = sum_i [B.'*C]_ii = sum_i sum_j B_ij * C_ij
which is the sum of element-wise products and therefore equivalent to B(:).'*C(:).
When explicitly computing trace( B.'*C ) you are actually pre-computing all k-by-k entries of B.'*C only to use the diagonal later on. AFAIK, Matlab does not optimize its calculation to save it from computing all the entries.

Here's a way

for ii = 1:n
    B = myBasis{ii};  
    for jj = ii:n
        C = myBasis{jj};
        t = full( B(:).'*C(:) ); % equivalent to trace(B'*C)!
        Q(ii,jj) = t;
        Q(jj,ii) = t;
    end
end

PS,
It is best not to use i and j as variable names in Matlab.

PPS,
You should notice that ' operator in Matlab is not matrix transpose, but hermitian conjugate, for actual transpose you need to use .'. In most cases complex numbers are not involved and there is no difference between the two operators, but once complex data is introduced, confusing between the two operators makes debugging quite a mess...

edited May 23 '17 at 12:05

Community

1
1

answered Jul 24 '14 at 13:28

Shai

111,146
38
238
371

You must use `full(t)` to make it equal to the matlab implementation of trace. `trace` converts `t` to a full matrix again. This is probably a time thief as well. @Adam are you sure sparse matrices are needed here? 1226*50*50 is 3 million elements. For double it means a little less than 30Mb. That is not so bad. – patrik Jul 24 '14 at 13:39
Yes, but `myBasis` is only a 1D array, temporary storing the other elements will not do a serious damage. So I think it is 1226*50*50. Correct me if I miss something. – patrik Jul 24 '14 at 13:47
@patrik I take it back. thanks for the `full` comment. you might comment the question directly about converting it all to full. – Shai Jul 24 '14 at 13:49
@Shai are you sure this works for sparse matricies as well? I dont understand why I need the full command. Sry - I get it now. – Adam Jul 24 '14 at 13:49
1

@Adam you only `full` the result - which is a scalar. a representation thing. you can try without it. – Shai Jul 24 '14 at 13:59

score 1 · Answer 2 · answered Jul 24 '14 at 13:23

Well, a couple of thoughts

1) Basic stuff: A'*B = (B'*A)' and trace(A) = trace(A'). Well, only this trick cut your calculations by almost 50%. Your Q(i,j) matrix is symmetric, and you only need to calculate n(n+1)/2 terms (and not n²)

2) To calculate the trace you don't need to calculate every term of B'*C, just the diagonal. Nevertheless, I don't know if it's easy to create a script in Matlab that is actually faster then just calculating B'*C (MatLab is pretty fast with matrix operations).

But I would definitely implement (1)

How to calculate matrix entries efficently using Matlab

2 Answers2