How to calculate this in matlab without for loop

Question

Since matlab is slow when executing for loop, I usually avoid for loop for all my codes, and turn those into matrix calculation, which would be way fast. But here is a problem I can not find a smart way:

I have a n x n matrix

A=[a1,a2,a3,...,an], 

here a1,a2,a3....an are columns of the matrix.

Another n x n matrix

B=[b1,b2,b3,...,bn], 

similarly b1,b2,b3... are also columns of B.

And also a n x n matrix M.

I want to calculate the n x n matrix

C=[c1,c2,c3,...,cn], 

thus (M+diag(ai))*ci = bi.

namely

ci = (M+diag(ai))\bi.

I know one way without for loop is:

C(:)=( blkdiag(M)+diag(A(:)) )\B(:).

But this will do too much calculation than needed.

Any smart solutions? you can assume there is no singularity problem in the calculation.

Answer 1

The statement "for loops are slow in Matlab" is no longer generally true since Matlab...euhm, R2008a? (someone please fill me in on this :)

Anyway, try this:

clc
clear all

M = rand(50);
a = rand(50);
b = rand(50);

% simple loop approach
tic
c = zeros(size(b));
for ii = 1:size(b,2)
    c(:,ii) = ( M+diag(a(:,ii)) ) \ b(:,ii);
end
toc

% not-so-simple vectorized approach
tic
MM = repmat({M}, 50,1);
c2 = (blkdiag(MM{:})+diag(a(:)))\b(:);
toc

norm(c(:)-c2(:))

Results:

Elapsed time is 0.011226 seconds.  % loop
Elapsed time is 1.049130 seconds.  % no-loop
ans =
     5.091221148787843e-10         % results are indeed "equal"

There might be a better way to vectorize the operation, but I doubt it will be much faster than the JIT'ed loop version.

Some problems are just not suited for vectorization. I think this is one.