5

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.

Min Lin
  • 3,177
  • 2
  • 19
  • 32
  • 5
    Have you tried to benchmark a for loop vs what you already have? sometimes it is faster than vectorized code... – bla Dec 03 '12 at 07:06
  • What do you mean by `diag(ai)`? `ai` is a scalar, right? – Tim Dec 03 '12 at 08:51
  • 1
    It is generally much slower to sit waiting for someone (including yourself) to devise a smart and fast loop-free code than to use the loop-based equivalent; Matlab is far more about optimising your time than optimising the CPU's time. And your belief that vectorised codes always (or generally) outperform their loop-based equivalents is much less true than it used to be. Follow @natan's advice. – High Performance Mark Dec 03 '12 at 09:15
  • 1
    You are solving linear systems of equations, possibly using dense matrices. That is an expensive computation. There is no way the loop overhead will be visible for even small systems. Why don't you benchmark your code before optimizing it? – angainor Dec 03 '12 at 10:13
  • @TimN: I think `ai` is the `i`-th column of `a`, so `diag(ai)` would be an all-zero matrix, except for `ai` on the main diagonal – Rody Oldenhuis Dec 03 '12 at 11:27

1 Answers1

2

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.

Rody Oldenhuis
  • 37,726
  • 7
  • 50
  • 96
  • It is expected that this vectorized version can be slow, because there are a lot more calculations (inverse of a 2500x2500 matrix). I know that for loop would be nothing compared to inverse, I'm just curious whether there is a way. – Min Lin Dec 03 '12 at 12:39