Rational Krylov Method（克雷洛夫法）求解特征值问题matlab代码示例（和Arnoldi方法比较）

最基本的特征值问题分为三类：
1、标准的线性特征值问题：
$Ax = \lambda x,A \in {\mathbb C^{n*n}}$
2、普遍的线性特征值问题：
$Ax = \lambda B x,A、B \in {\mathbb C^{n*n}}$
3、普遍的艾米特正定线性特征值问题：
$Ax = \lambda B x,A、B \in {\mathbb C^{n*n}}$
$A^*=A,B^*=B > 0 \in {\mathbb C^{n*n}}$

不妨做一个程序将RKM方法和Arnoldi方法做一个比较。

% demo routine for Rational Krylov vs Arnoldi with shift-invert
clc
clear
n = 1000;
m = 40; % subspace size
A = diag([1:1:n]); % testing matrix A of eigenvalues 1, 2, .., n

% shifts: define more shift as you like
sigma1 = 101.5;
sigma2 = 121.5;


% --  Rational Krylov with shifts SIGMA
m1 = floor(m/2); m2 = m-m1;
SIGMA = [repmat(sigma1, m1, 1); repmat(sigma2, m2, 1)];
%SIGMA = repmat([sigma1; sigma2], m/2, 1); % even number m

[Q1, K, L] = rarnoldi(A, SIGMA, m);
[V1, E1] = eig(L(1:m,:), K(1:m,:));

% pick eig of relative residual norm < tol
tol = 1.0E-6;
nA = norm(A,1);
e1 = [];
for i = 1:m
    lam = E1(i,i); 
    v = Q1*(K*V1(:, i));
    res = norm(A*v-lam*v)/norm(v)/nA;
    if res<tol, e1 = [e1, lam]; end
end

% -- shift invert Arnoldi with shift sigma1 
sigma = sigma1;
[LL, UU, PP] = lu(A-sigma*speye(n));
Afun = @(x) UU \ (LL \ (PP*x));
[Q2, H] = arnoldi(Afun, n, m);
[V2, E2] = eig(H(1:m,:)); E2 = 1./E2 + sigma;

% pick eig of relative residual norm < tol
tol = 1.0E-6;
nA = norm(A,1);
e2 = [];
for i = 1:m
    lam = E2(i,i); 
    v = Q2(:,1:m)*V2(:, i);
    res = norm(A*v-lam*v)/norm(v)/nA;
    if res<tol, e2 = [e2, lam]; end
end

% -- display results
figure
plot(real([sigma1, sigma2]), imag([sigma1, sigma2]), 'or', 'DisplayName', 'shifts', 'MarkerSize', 8)
hold on;
plot(real(e2), imag(e2), 'xb', 'DisplayName', 'Arnoldi', 'MarkerSize', 8);
plot(real(e1), imag(e1), '+k', 'DisplayName', 'Rational Krylov', 'MarkerSize', 8);
legend show
xlabel('real'); ylabel('imag')
title('approximate eigenvalues')

% END

用到的一些子函数就一起列在下面了：

function [V, H] = arnoldi(Afun, n, m)
% function [V, H] = arnoldi(Afun, n, m) produce an Arnoldi decomposition
%   of order m for a square matri A. 
%       Afun(v) = A*v
%       n = dim of A
%
v = randn(n,1);
v = v / norm(v);

V = v;
H = zeros(m+1, m);

for i = 1:m
    w = Afun(v);
    h = V'*w;
    w = w - V*h;
    gamma = norm(w);

    if gamma==0
        return
    end

    v = w/gamma;
    V = [V, v];
    H(1:i,i) = h;
    H(1+i, i) = gamma;   
end



function [V, K, L] = rarnoldi(A, SIGMA, m)
% function [V, K, L] = rarnoldi(Afun, n, m) produce an rational Krylov
%   decomposition of order m for a square matri A. 
%       SIGMA: length m vector containing shifts
%

n = size(A,1);
v = randn(n,1);
v = v / norm(v);

V = v;
K = zeros(m+1, m);
L = zeros(m+1, m);

for i = 1:m
    sigma = SIGMA(i);
    if sigma ~= inf
        w = (A-sigma*speye(n))\v;
        h = V'*w;
        w = w - V*h;
        gamma = norm(w);
        v = w/gamma;

        V = [V, v];
        K(1:i, i) = h;
        K(i+1, i) = gamma;
        L(1:i+1, i) = sigma*K(1:i+1, i);
        L(i, i) = L(i, i) + 1;
    else
        w = A*v;
        h = V'*w;
        w = w - V*h;
        gamma = norm(w);
        v = w/gamma;

        V = [V, v];
        K(i, i) = 1;
        L(1:i, i) = h;
        L(i+1, i) = gamma;
    end
    % breakdown case L(i+1, i) = K(i+1,i) = 0 is not treated.
end