%%
clear
clc
rng(6)

%%
% We illustrate a few steps of the Arnoldi iteration for a small matrix.
np = 6;
A = magic(np)

%%
% The seed vector determines the first member of the orthonormal basis.
u = randn(np, 1);
Q = u/norm(u)

%%
% Multiplication by A gives us a new vector in K_2. 
Aq = A*Q(:, 1)

%% 
% We subtract off its projection in the previous direction. The remainder
% is rescaled to give us the next orthonormal column.
v = Aq - (Q(:,1)'*Aq)*Q(:,1);
Q(:, 2) = v/norm(v);

%%
Q
norm(Q'*Q - eye(2))

%%
% On the next pass, we have to subtract off the projections in two previous
% directions.
Aq = A*Q(:, 2);
v = Aq - (Q(:,1)'*Aq)*Q(:,1) - (Q(:,2)'*Aq)*Q(:,2);
Q(:, 3) = v/norm(v);
Q

%%
% At every step, Q_m is an ONC matrix.
norm(Q'*Q - eye(3))