Optimizing the weighted sum of sparse matrices

I welcome any help for the following code optimization problem:

I have a collection of N sparse matrices of identical sizes ([s1 s2]) stored in a cell array A and a corresponding number of scalar weights stored in an vector w. I want to compute the sum of all the matrices in A weighted by the values stored in w. Through the iterations of my program, only the values in wchange. I can therefore compute a priori the number of non-zero elements in my result and pre-allocate some memory for it using spalloc.

For the moment I have something like:

result = spalloc(s1,s1,number_of_non_zero);
for i=1:N
    result = result + w(i)*A{i};
end

I really need do optimize this part which for the moment takes most of the computing time in my program (checked with the profiling tools).

Some additional information:

-The above code runs millions of times so even minor improvements are welcome.

-The matrices in A come from a finite element code (1D or 2D)

-I have no problem moving away from the cell structure if I can save some time (like using cell2mat(A))

Thank you for any hint on how to speed up this part o开发者_C百科f the code.

A.

---

I can't say for sure if this solution will be more computationally-efficient, but it's something else to try...

If you really have to represent your matrices as sparse (i.e. full matrices take up too much memory) and the built-in sparse representation in MATLAB isn't giving you the performance you desire, then you can try representing the sparse matrices in a different way. Specifically, you can represent them as N-by-3 matrices where the first two columns contain the row and column indices into the matrix for all the non-zero values and the third column contains the non-zero values. You can convert your data to this form using the function FIND like so:

for iMatrix = 1:numel(A)
  [r,c,v] = find(A{iMatrix});
  A{iMatrix} = [r c v];
end

Each time you need to compute the weighted sum of these matrices, you first need to multiply the values by the weights:

B = A;  %# Store a temporary copy of A
for iMatrix = 1:numel(B)
  B{iMatrix}(:,3) = w(iMatrix).*B{iMatrix}(:,3);
end

Then, you can compute the final sum using the function ACCUMARRAY:

B = vertcat(B{:});  %# Convert B from a cell array to an N-by-3 matrix
result = accumarray(B(:,1:2),B(:,3));

The variable result in this case will be a full matrix. If you need result to be a sparse matrix, you can add extra arguments in the call to ACCUMARRAY like so:

result = accumarray(B(:,1:2),B(:,3),[],[],[],true);

---

If you convert A to a matrix instead of a cell array then you could vectorize the loop using some RESHAPE and REPMAT acrobatics. Pretend we have the following data:

>> A(:,:,1) = [1 0 0; 0 0 2; 0 0 0];
>> A(:,:,2) = [3 0 1; 0 3 0; 0 1 0]

A(:,:,1) =

     1     0     0
     0     0     2
     0     0     0


A(:,:,2) =

     3     0     1
     0     3     0
     0     1     0

>> w = [2; 3]

w =

     2
     3

Reshape w so that you can do an element-by-element multiplication and then sum:

>> w = reshape(w, [1 1 length(w)])

w(:,:,1) =

     2


w(:,:,2) =

     3

>> w = repmat(w, [size(A,1) size(A,2) 1])

w(:,:,1) =

     2     2     2
     2     2     2
     2     2     2


w(:,:,2) =

     3     3     3
     3     3     3
     3     3     3

>> w .* A

ans(:,:,1) =

     2     0     0
     0     0     4
     0     0     0


ans(:,:,2) =

     9     0     3
     0     9     0
     0     3     0

>> sum(w .* A, 3)

ans =

    11     0     3
     0     9     4
     0     3     0