Does the vectorization method can be used in the vector iteration?

Hi,guys. I use the following codes to calculate a problem of vector iteration:

M0=zeros(n,m);
for k=1:m
     u=M*u;
     M0(:,k)=u;    
end

u is a column vector with n elements, M is a n*n transfer matrix. The iteration is carried out by m times, the result is saved in the matrix M0.Using the for loop, this calculation can be done easily. But when m is large, the calculation time is so long. I want use the vectorization method to improve this program and reduce the time. Is this possible?

[Matt J]

Hmmmm. In case it matters, what manipulations are you planning to do with M0 once you've built it?

[nan guaa]

M0 just a matrix I use to save the results. All I want are the m vectors, which are used to do contour plots with different coordinates.

Now, I doubt the iteration whether can be implemented by vectorization method. Because iteration use the former result to get the later one. Maybe this can not be done by vectorization essentially.

Anyone, please give me some advices, thanks.

Answer 1

[Walter Roberson]

Your values are M^k * u0 for k = 1 to m. The obvious optimization is to use SVD

[U, S, V] = svd(M);

after which M^k = U * S^k * V' and since S is diagonal, the S^k is easy to calculate.

You could calculate V' * u0 as a single matrix. The iterative multiplication by the diagonal matrix S would be multiplying the first row by S(1,1), the second row by S(2,2), the third row by S(3,3) and so on.

Now in theory you can calculate that diagonal multiplication faster than a full matrix multiplication, but unless you create some mex code for it, you would have the overhead of the MATLAB interpretation for every operation, so this is not necessarily going to be any faster than just letting MATLAB do the matrix multiplication.

[nan guaa]

Thanks for your reply.

Is there some other method to accelerate iteration?

[Walter Roberson]

   bsxfun(@power, diag(S), (1:m).')

But then you have to break up by row and do the equivalent of U * diag(therow) * V' * u0

Ah since the V' * u0 is going to have to be a vector then the diag() to expand to a matrix followed by the right multiply is going to be equivalent to using .* of the two vectors. So calculate V' * u0 to get a vector, repmat it to m rows long, use .* with the above matrix of values. Then all that would be left would be the equivalent of taking one riw at a time of that and left multiply by U. There might be a more efficient way of doing that, I would have to think about it