%% Rayleigh quotient iteration for an operator
% Nick Hale and Yuji Nakatsukasa, March 2017

%%
% (Chebfun example ode-eig/RayleighQuotient.m)

%%
% [revised July 2019]

%% 1. Symmetric matrix
% The Rayleigh quotient iteration (RQI) is a well-known algorithm for computing
% an eigenpair of a matrix $A$ (symmetric or nonsymmetric). Using an approximate
% $\tilde\lambda$ and approximate eigenvector $\tilde x$ (normalized by
% $\|\tilde x\|_2=1$), it updates the approximate eigenvector via one step of
% the shifted-and-inverted power method, $\tilde x:=(A-\tilde\lambda
% I)^{-1}\tilde x$, which is computed by solving a linear system (then
% normalized again). The approximate eigenvalue is then updated via the Rayleigh
% quotient $\tilde\lambda = \tilde x^\ast
% A\tilde x$. Under mild assumptions, by repeating the process, $\tilde x$
% converges to an eigenvector, usually to the one corresponding to the
% eigenvalue that the initial $\tilde\lambda$ is closest to. The convergence of RQI
% is known to be asymptotically cubic when $A$ is symmetric [1, Sec. 4.6], and
% otherwise quadratic.

%%
% Here is an example showing cubic convergence for a random symmetric $10\times
% 10$ matrix. (Note that we sacrifice some efficiency in computing $Au$ in the
% residual computation for the sake of clarity. This could easily be avoided.)

rng(10)                                      % choose a random number seed
tol = 1e-10;
n = 10;
A = randn(n); A = A'+A;                      % symmetric matrix
I = eye(size(A));                            % identity

% Initial guesses:
disp('lam:')
lam = A(end,end); disp(lam)                  % seek eigenvalue near lam
u = rand(n,1);  u = u/norm(u);               % random guess
res = norm(A*u - lam*u)/norm(A*u);           % initial residual

% Rayleigh quotient iteration:
while ( res(end) > tol )
    u = (A - lam*I)\u; u = u/norm(u);        % core RQI
    lam = u'*A*u; disp(lam)                  % update Rayleigh quotient
    res = [res; norm(A*u-lam*u)/norm(A*u)];  % store residual
end
res

%% 2. Nonsymmetric matrix
% Next let's try a nonsymmetric matrix. The convergence becomes quadratic. 

A = randn(n);                                 % nonsymmetric matrix
% Initial guesses:
disp('lam:')
lam = A(end,end); disp(lam)                   % seek eigenvalue near lam
u = rand(n,1)+1i*randn(n,1); u = u/norm(u);   % random guess (eigval can be complex)
res2 = norm(A*u - lam*u)/norm(A*u);           % initial residual

% Rayleigh quotient iteration:
while ( res2(end) > tol )
    u = (A - lam*I)\u; u = u/norm(u);         % core RQI
    lam = u'*A*u; disp(lam)                   % update Rayleigh quotient
    res2 = [res2; norm(A*u-lam*u)/norm(A*u)]; % store residual
end
res2

%%
% We can recover the cubic convergence by using the (conjugate) transpose and
% running a two-sided Rayleigh quotient iteration. Note that the algorithm is
% equivalent to RQI when $A$ is symmetric.

A = randn(n);                                 % nonsymmetric matrix
% Initial guesses:
disp('lam:')
lam = A(end,end); disp(lam)                   % seek eigenvalue near lam
u = rand(n,1)+1i*randn(n,1); u = u/norm(u);   % random guess for right evec
v = rand(n,1);  v = v/norm(v);                % random guess for left evec
res3 = norm(A*u - lam*u)/norm(A*u);           % initial residual

% two-sided Rayleigh quotient iteration: 
while ( res3(end) > tol )
    u = (A -lam*I)\u; u = u/norm(u);          % core RQI
    v = (A'-lam*I)\v; v = v/norm(v);          % core RQI for left eigvec    
    lam = (v'*A*u)/(v'*u); disp(lam)          % update Rayleigh quotient
    res3 = [res3; norm(A*u-lam*u)/norm(A*u)]; % store residual
end
res3

% Plot convergence rates:
semilogy(res, 'b-o'), hold on
text(length(res) + .1, res(end), 'symm', 'color', 'b')
semilogy(res2, 'r--x' )
text(length(res2) - .9, res2(end-1), 'nonsymm', 'color', 'r')
semilogy(res3, 'm--^' )
text(length(res3) - .9, res3(end-1), 'nonsymm two-sided', 'color', 'r')
xlabel('iteration'), ylabel('residual')
grid on, hold off

%% 3. Selfadjoint linear operator
% Now we explore the use of RQI for a linear operator represented by a chebop.
% Let us consider the selfadjoint operator $Au = -u''$ with selfadjoint
% Dirichlet boundary conditions $u(-\pi/2) = u(\pi/2) = 0$. As an initial guess
% we use a random function generated by the randnfun command. Note that the
% code for the RQI is almost identical to the matrix case above! (This is
% enabled by a Chebfun feature that allows additions of chebops, along with
% eye(A) for the identity operator on the domain of A.)

dom = [-pi/2, pi/2];
A = chebop(@(u) -diff(u, 2), dom, 0);       % self-adjoint operator
I = eye(A);                                 % identity on dom; equivalent
                                            % to I = chebop(@(u)u,dom)

% Initial guesses:
disp('lam:')
lam = 3.8; disp(lam)                        % seek eigenvalue near lam 
u = randnfun(.1, dom);  u = u/norm(u);      % random guess for eigenfunction
res = norm(A*u-lam*u)/norm(A*u);            % initial residual

% Rayleigh quotient iteration:
while ( res(end) > tol )
    u = (A-lam*I)\u; u = u/norm(u);         % core RQI
    lam = u'*A*u; disp(lam)                 % update Rayleigh quotient
    res = [res; norm(A*u-lam*u)/norm(A*u)]; % store residual
end

%%
% From the residual output, we can guess that the convergence is cubic (it is
% actually too fast to verify).

res

%% 4. Non-selfadjoint linear operator
% Now consider the non-selfadjoint operator $Au = -u'' + u' + u$, again with
% zero Dirichlet boundary conditions:

dom = [-pi/2, pi/2];
x = chebfun('x',dom);
A = chebop(@(x,u) -diff(u,2) + diff(u) + u, dom, 0);
I = eye(A);

% Initial guesses:
disp('lam:')
lam = 1; disp(lam)
u = randnfun(.1, dom); u = u/norm(u);
res2 = norm(A*u-lam*u)/norm(A*u);             % initial residual

% Rayleigh quotient iteration:
while ( res2(end) > tol)
    u = (A-lam*I)\u; u = u/norm(u);           % core RQI
    lam = u'*A*u; disp(lam)                   % update Rayleigh quotient
    res2 = [res2; norm(A*u-lam*u)/norm(A*u)]; % store residual
end

%%
% From the residual output, we can see that here the convergence is quadratic:

res2

%%
% As before, let's try to improve the convergence to cubic. This involves the
% adjoint, which plays the role of the conjugate transpose in the matrix case.

% Initial guesses:
disp('lam:')
lam = 1; disp(lam)
u = randnfun(.1, dom); u = u/norm(u);         % random guess for right eigfunc
v = randnfun(.1, dom); v = v/norm(v);         % random guess for left eigfunc
res3 = norm(A*u-lam*u) / norm(A*u);

% Rayleigh quotient iteration:
while ( res3(end) > tol)
    u = (A -lam*I)\u; u = u/norm(u);          % core RQI
    v = (A'-lam*I)\v; v = v/norm(v);          % adjoint RQI    
    lam = v'*(A*u)/(v'*u); disp(lam)          % update Rayleigh quotient
    res3 = [res3; norm(A*u-lam*u)/norm(A*u)]; % store residual   
end

res3

%%
% RQI appears to have computed eigenpairs for both operators, selfadjoint and
% non-selfadjoint. Now let's examine the convergence rates by plotting the
% residual convergence.

semilogy(res, 'b-o'), hold on
text(length(res) + .1, res(end), 'selfadj', 'color', 'b')
semilogy(res2, 'r--x' )
text(length(res2) - .9, res2(end-1), 'non-selfadj', 'color', 'r')   
semilogy(res3, 'm--^' )
text(length(res3) - .9, res3(end-1), 'non-selfadj two-sided', 'color', 'm')   
xlabel('iteration'), ylabel('residual')
grid on

%% 5. References
%
% 1. B. N. Parlett, _The Symmetric Eigenvalue Problem_, SIAM, 1996.
