Does the vectorization method can be used in the vector iteration?
2 views (last 30 days)
Show older comments
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?
2 Comments
Matt J
on 28 Aug 2015
Hmmmm. In case it matters, what manipulations are you planning to do with M0 once you've built it?
Answers (1)
Walter Roberson
on 29 Aug 2015
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.
2 Comments
Walter Roberson
on 29 Aug 2015
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
See Also
Categories
Find more on Matrix Indexing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!