1

I'm completely new to MATLAB, and I can't figure out how to do the following: I have

A=[-1,0.4,0.8;1,0,0;0,1,0];
b=[0;0.3;6];

What I want to do is define, for $N$ fixed but arbitrary, a $3\times N$ matrix C by having the $i$th column of C be $A^{(N-i)}b$. Without doing the calculations separately beforehand, how can I have MATLAB compute and set C as this matrix?

littleO
  • 51,938
Jeh
  • 414

2 Answers2

2
A=[-1,0.4,0.8;1,0,0;0,1,0];
b=[0;0.3;6];
N=10;
Ai=eye(3); % Identity 3 x 3 matrix
C=zeros(3,N);
for i=1:N
    C(:,N-i+1)=Ai*b;
    Ai=Ai*A;
end

Note that the first column of C will be equal to $A^{N-1}b$, while the last column will be equal to $b$. If you wanted them to be $A^Nb$ and $Ab$ respectively then replace the fourth line of the code with

Ai=A;

Edit: The code in littleO's answer is much faster. Mine is an example of what you shouldn't do (at least when you care about efficiency)!

posilon
  • 2,253
  • The matrix multiplication at each iteration is expensive! This should be avoided. – littleO Oct 26 '13 at 02:18
  • For the given example where A is just 3 x 3, the matrix multiplication is not expensive. So for large N (number of iterations) my code runs actually faster than that of Nasser. If you are interested in doing the calculation for matrices of bigger dimension then my code is slower. I don't know what algorithms MATLAB uses for matrix multiplication or exponentiation but it seems that the best choice depends on whether you are optimizing with respect to N or the size of A. How can A^i be calculated faster than a multiplication of two n x n matrices? – posilon Oct 26 '13 at 02:54
  • I wrote an answer that uses only one matrix-vector multiplication at each iteration. Your code is certainly faster than nasser's (about 6.5 times faster in this example), but nasser's code is very inefficient. – littleO Oct 26 '13 at 02:56
  • You are right, your code is much faster! But my code was actually slower than Nasser's for bigger dimensions. Which means that A^i can be computed faster than a matrix multiplication. Do you know how? – posilon Oct 26 '13 at 03:00
  • I'm observing that your code is much faster than Nasser's for larger dimensions. – littleO Oct 26 '13 at 03:13
  • Well, I ran it again. For small N you are right. For bigger (say N=200, dimension = 100) it gets worse! – posilon Oct 26 '13 at 03:25
  • hmm, I just ran an example where $A$ is $100 \times 100$ and $M = 200$, and your code took .064 seconds, whereas Nasser's took .345 seconds. (I think Nasser took down his code, which did have the advantage of being a single line.) – littleO Oct 26 '13 at 03:30
  • Well, I would take down mine since it has no advantage over yours. But then I thought it has pedagogical value. It teaches you what not to do! – posilon Oct 26 '13 at 03:51
  • haha, no you should leave it up. I think this has been an interesting example. (And I think somebody else might have more to say about computing Krylov matrices stably.) – littleO Oct 26 '13 at 03:59
2

This code avoids having a matrix-matrix multiplication at each iteration.

A = [-1,0.4,0.8;1,0,0;0,1,0]; 
b = [0;0.3;6];
N = 5;

C = zeros(3,N);
vec = b;

for i = N:-1:1

    C(:,i) = vec;
    vec = A*vec;

end

Here's another option in Matlab:

C = fliplr(gallery('krylov',A,b,N));
littleO
  • 51,938