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.
