Efficiently multiplying diagonal and general matrices
28 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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!
0 commentaires
Réponse acceptée
Teja Muppirala
le 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.
Plus de réponses (1)
Teja Muppirala
le 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 commentaires
Jan
le 20 Sep 2013
15% faster by avoiding sqrt:
B = bsxfun(@times,M, diag(D));
B = M.' * B;
Voir également
Catégories
En savoir plus sur Operating on Diagonal Matrices dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!