Efficiently multiplying diagonal and general matrices

I wish to find the most efficient way to implement the following equation

M'*D*M

where M is a m*n dense rectangular matrix (with no specific structure), and D is a m*m diagonal matrix with all positive elements. In addition, m >> n, and M is constant throughout the course of the algorithm, with only the elements of D changing.

I know there are tricks for a related problem (D*M*D) to reduce the number of operations considerably, but is there one for this problem? Ideally is there a way to factorize / rearrange this so I can compute M' * M offline (or something similar), and update D at each iteration?

Thanks!

Answer 1

[Teja Muppirala]

The best solution is going to depend on what your m and n actually are (if you know representative values of them, you should include those in your problem statement).

Thinking generally though:

I am almost certain you can't just find M'*M and somehow do something efficiently with only that. But you can do something similar. Notice how this expression is linear in the entries of D.

You can express D as a sum of elementary basis functions

D = d1*e1 + d2*e2 + ... + dm*em

where dk, a scalar, is the kth diagonal entry of D, and ek is a [m x m] matrix with all zeros except for a 1 in the kth position along the diagonal.

Then:

M'*D*M

= M'*(d1*e1 + d2*e2 + d3*e3 + ... + dm*em)*M

= d1 * (M'*e1*M) + d2 * (M'*e2*M) + ... + dm * (M'*em*M)

This implies that if you calculate all the M'*ek*M beforehand, then you just need to take a linear combination of them.

But each M'*ek*M is simply M(k,:)'*M(:,k).

I will calculate these offline and store them in an 3-d array "J". I reshape J to an [(n^2) x m] matrix since we want to take linear combinations of its columns by postmultiplying it with the elements in D.

m = 100000;
n = 5;
M = rand(m,n);

J = zeros(n,n,m); % Preallocate J for n*n*m elements of storage

for k = 1:m
    J(:,:,k) = M(k,:)'*M(k,:);
end
J = reshape(J,n^2,m);

Now, I can use J to quickly calculate the answer for any D. We'll try all 3 methods.

d = rand(m,1); %Generate a new d (only the diagonal entries)

tic; D = sparse(1:m,1:m,d); A = M'*D*M; toc; % Method 1, direct multiplication
tic; B = bsxfun(@times,M,sqrt(d)); B = B.'*B; toc; % Method 2, using BSXFUN
tic; C = reshape(J*d,n,n); toc; % <-- Method 3, precalculating matrices.

norm(A-C)

Again, depending on what m and n actually are, the fastest method may be different (for this choice of m and n, it seems method 3 is somewhat faster). One drawback, however, is that you need to be able to store a dense [n x n x m] array, and this may not be feasible if the n and m are too large.

[Jonathan Currie]

Thanks Teja Method 3 worked out to be faster. In addition, I can exploit symmetry within M'*M and thus skip some of the rows in J*d, further reducing operations.

Answer 2

[Teja Muppirala]

M = randn(10000,10);
D = diag(randn(10000,1).^2);

tic
A = M'*D*M;
toc

tic
B = bsxfun(@times,M,sqrt(diag(D)));
B = B.'*B;
toc

[Jonathan Currie]

Thanks Teja for that, I have updated my question to reflect a further requirement which I don't think your solution completes?

[Jan]

15% faster by avoiding sqrt:

B = bsxfun(@times,M, diag(D));
B = M.' * B;