Efficiently multiplying diagonal and general matrices

51 views (last 30 days)
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!

Accepted Answer

Teja Muppirala
Teja Muppirala on 19 Sep 2013
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.
  1 Comment
Jonathan Currie
Jonathan Currie on 20 Sep 2013
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.

Sign in to comment.

More Answers (1)

Teja Muppirala
Teja Muppirala on 19 Sep 2013
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
  2 Comments
Jonathan Currie
Jonathan Currie on 19 Sep 2013
Thanks Teja for that, I have updated my question to reflect a further requirement which I don't think your solution completes?
Jan
Jan on 20 Sep 2013
15% faster by avoiding sqrt:
B = bsxfun(@times,M, diag(D));
B = M.' * B;

Sign in to comment.

Categories

Find more on Operating on Diagonal Matrices in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!