Efficiently Setting up Sparse Matrices in MATLAB for Matrix Inversion

Question

Consider the following weighted solution to the normal equation for a least squares inverse problem:

m = inv(G'*W'*W*G)*G'*W'*W*d

I would like to set up the weighting matrix W, which is a square diagonal matrix with weights on the diagonal.

As I have a large number of data points in d (10⁷), my system matrix G is also large, but only in one dimension (as I have far more data points than model parameters). In the case of 6 model parameters, G is of size (10⁷ × 6). Hence, W has to be of size (10⁷ × 10⁷). However, it is sparse, with only 10⁷ non-zero entries (the weights).

To reduce memory, I use sparse on W.

To assign weights, I do the following

d = [d1;d2];   
W = sparse(length(d),length(d)) 
w1 = length(d2)/length(d);
w2 = length(d1)/length(d);
W(1:length(d)+1:length(d)*length(d1)) = w1;
W(length(d)*length(d1)+1:length(d)+1:end) = w2;

d1 and d2 are column vectors with observations.

This will assign weights to the diagonal, but it is awfully slow.

My question:

Can I either

Somehow speed up the assignments of weights to the diagonal, or
Re-write m = inv(G'*W'*W*G)*G'*W'*W*d so that I do not have to set-up W at all?

Note 1: In the shown weights are two different constants, but in practice they will vary allowing the diagonal!

Note 2: The bottle neck of the code is indeed setting up W, not the inversion itself, as the inverted matrix is only of size (6 × 6).

What are d1 and d2? Please post runnable code. Also, d = [d1;d2] implies w1 and w2 are just ones and W is eye? Why two diagonal assignments (last two lines)? — Luis Mendo
d = [d1;d2] is a vertical concatenation of n x 1 and m x 1 vectors and as such does not imply anything on the length of d1 and d2 except length(d1)+length(d2) = length(d) — TheodorBecker
Have you tried using W = sparse([1:numel(d1) 1:numel(d2)], 1:numel(d), [w1; w2], numel(d), numel(d));, where w1 and w2 are column vectors? This replaces the intiallization W = sparse(length(d),length(d)); and the two assignment lines — Luis Mendo
@LuisMendo: Great, I did not realize such initialization of sparse matrices was possible. It's actually even given in the documentary. It is faster by several orders of magnitude. Thanks a lot! — TheodorBecker
If W is truly diagonal, then you can better do bsxfun(@times, G, [length(d2) length(d1)]/length(d)) to compute G.' · W.'. What are WG and Wd? — Rody Oldenhuis

chtz chtz · Accepted Answer · 2018-03-26T01:08:11

By the structure of W I'd say it is more efficient to split G into two halves and calculate the products with W as scalar products:

G1 = G(1:length(d1),:);
G2 = G(length(d1)+1:end,:);

m = (w1^2*G1'*G1 + w2^2*G2'*G2) \ (w1^2*G1'*d1 + w2^2*G2'*d2);

Also, generally use A \ b instead of inv(A)*b (even if A is small).

Efficiently Setting up Sparse Matrices in MATLAB for Matrix Inversion

1 Answers