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

最基本的特征值问题分为三类：

1、标准的线性特征值问题：

Ax=λx,A∈Cn∗n

2、普遍的线性特征值问题：

Ax=λBx,A、B∈Cn∗n

3、普遍的艾米特正定线性特征值问题：

Ax=λBx,A、B∈Cn∗n

A∗=A,B∗=B>0∈Cn∗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包含的数据如下形式：



运行结果示例：