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!
采纳的回答
Teja Muppirala
2013-9-19
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 个）
Teja Muppirala
2013-9-19
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 个评论
Jonathan Currie
2013-9-19
Jan
2013-9-20
15% faster by avoiding sqrt:
B = bsxfun(@times,M, diag(D));
B = M.' * B;
类别
在 帮助中心 和 File Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息
2013-9-19
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
