加速子空间迭代法（Accelerated Subspace Iteration）求特征值问题matlab程序

最基本的特征值问题分为三类：
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}}$

下面给出用加速子空间迭代法求解特征值问题的matlab代码：

% Demo routines for the subspace iteration with Rayleigh-Ritz projection.

rng(0);

load barbell.mat;
%   A matrix pencil arising from solving the Laplacian eigenvalue 
%   problem in a barbell shaped domain.

n = size(A,1); % matrix size
J = 5; % subspace size

%FILTER = 'Exact spectral projector';
FILTER = 'Power iteration'; 
%FILTER = 'Multi-step power iteration';

% --- Define filters 
switch FILTER
    case 'Exact spectral projector'
        % 1. rhoM(Q) = exact projector
        [V, E] = eigs(A, B, J);
        L = chol(V' * B * V); 
        V = V / L'; % B orthonormalization
        rhoM = @(Q) V*(V'*(B*Q)); 
    case 'Power iteration'
        % 2. rhoM(Q) = M*Q. 
        [L, U, P] = lu(B);
        rhoM = @(Q) U\(L\(P*(A*Q))); 
    case 'Multi-step power iteration'
        % 3. rhoM(Q) = M*M*Q.
        [L, U, P] = lu(B);
        rhoM0 = @(Q) U\(L\(P*(A*Q))); 
        rhoM = @(Q) rhoM0( rhoM0(Q) );
    otherwise
        error('TESTING CASE NOT FOUND')
end


% Subspace iteration with Rayleigh-Ritz projection
maxit = 100;
resd = zeros(maxit,J);
Q = randn(n,J);
nA = norm(A,1);
nB = norm(B,1);    
for j=1:maxit
    Y = rhoM(Q);
    AA = Y' * (A * Y);
    BB = Y' * (B * Y);
    [VV, EE] = eig(AA, BB, 'chol');
    BBB = VV' * BB * VV; BBB = (BBB+BBB')/2;
    LL = chol(BBB); 
    VV = VV / LL'; % B orthonormalization
    Q = Y * VV;

    % evaluate relative residual norm of each Ritz pair
    for jj = 1:J
        r = A*Q(:,jj) - B*Q(:,jj) * EE(jj,jj);
        res = norm(r) / norm(Q(:,jj)) / (nA + nB*abs(EE(jj,jj)));
        resd(j,jj) = res;
    end

end

% residual plot
figure
semilogy(sort(resd')', '+')
xlabel('iteration')
ylabel('relative residual norms')
title(FILTER);

% eigenvalue plot
figure
ee = eig(full(A),full(B));
plot(real(ee), imag(ee), 'or', 'DisplayName', 'exact eigval');
hold on;
eec = diag(EE);
plot(real(eec), imag(eec), '+b', 'DisplayName', 'computed eigval');
legend show;
title(FILTER);

%END
% 

A是输入的大规模矩阵，求出来的EE和VV蕴含了部分的特征值和特征向量。程序中的所使用的FILTER可以改变。barbell.mat包含的数据如下形式：

运行结果示例：