Matlab 语法记录(I)——特征值排序问题

1. 矩阵的特征值分解

矩阵的特征值分解是非常重要的数学工具。在matlab中一般使用eig()函数完成矩阵特征值和特征向量的提取，其用法如下

A = magic(4);
[V,D] = eig(A);

结果如下：

A =
16 2 3 13
5 11 10 8
9 7 6 12
4 14 15 1
V =
-0.5000 -0.8236 0.3764 -0.2236
-0.5000 0.4236 0.0236 -0.6708
-0.5000 0.0236 0.4236 0.6708
D =
34.0000 0 0 0
0 8.9443 0 0
0 0 -8.9443 0
0 0 0 -0.0000

2. eig() 和 eigs()

显然eig()就是一般意义上的计算矩阵的特征值和特征向量

E = eig(A) 返回方阵A的所有特征值，构成列向量E。
[V,D] = eig(A) 返回方阵A的特征值和特征向量，其中特征值排满对角阵D的对角元素，对应的特征向量排列为方阵V的每一列。

而eigs()也能求取矩阵的特征值和特征向量，不过其返回的方阵幅值最大的6个特征值和特征向量，用法和eig()类似。不过eigs()通过迭代的方式来求解特征值，所以其在加快运算速度的同时降低了准确度。另外，一般eigs()处理的大型稀疏矩阵。

[V,D] = eigs(A) 返回方阵A的前6个最大特征特征值和特征向量。
[V,D] = eigs(A,k) 返回前k个最大特征值和特征向量。

一般情况下，eig()和eigs()返回的特征值是从大到小排列，然而这并不是一定的。经过测试，两者的特征值排序都可能为乱序，所以，在对顺序有要求的情况下，需要通过sort()函数来自动排序。

3. 使用sort()函数解决特征值排序问题

如下

A = magic(6);
[V,D] = eig(A);
[D_sort,index] = sort(diag(D),'descend');
D_sort = D_sort(index);
V_sort = V(:,index);

按特征值大小排序结果如下：

A =
35 1 6 26 19 24
3 32 7 21 23 25
31 9 2 22 27 20
8 28 33 17 10 15
30 5 34 12 14 16
4 36 29 13 18 11
D =
111.0000 0 0 0 0 0
0 27.0000 0 0 0 0
0 0 -27.0000 0 0 0
0 0 0 9.7980 0 0
0 0 0 0 0.0000 0
0 0 0 0 0 -9.7980
D_sort =
111.0000
27.0000
-27.0000
-9.7980
9.7980
0.0000

4. Another Solution

John D’Errico设计了一个eigenshuffle.m函数能够得到排序后的特征值和特征向量。该方法排序方式为特征值大小降序排列。
速度测试：

%a test by Min Qin
A = magic(10);
iterate = 10000;
tic;
for i = 1:iterate
    [V,D] = eig(A);
    [D_sort,index] = sort(diag(D),'descend');
    V_sort = V(:,index);
end
toc;
tic;
for i = 1:iterate
    [V_sort,D_sort] = eigenshuffle(A);
end
toc;

结果：

Elapsed time is 0.325471 seconds.
Elapsed time is 0.447886 seconds.

显然eigenshuffle函数的速度比传统方法略低。

参考： https://cn.mathworks.com/matlabcentral/fileexchange/22885-eigenshuffle
eigenshuffle.m代码：

function [Vseq,Dseq] = eigenshuffle(Asequence)
% eigenshuffle: Consistent sorting for an eigenvalue/vector sequence
% [Vseq,Dseq] = eigenshuffle(Asequence)
%
% Includes munkres.m (by gracious permission from Yi Cao)
% to choose the appropriate permutation. This greatly`这里写代码片`
% enhances the speed of eigenshuffle over my previous
% release.
%
% http://www.mathworks.com/matlabcentral/fileexchange/20652
%
% Arguments: (input)
%  Asequence - an array of eigenvalue problems. If
%      Asequence is a 3-d numeric array, then each
%      plane of Asequence must contain a square
%      matrix that will be used to call eig.
%
%      Eig will be called on each of these matrices
%      to produce a series of eigenvalues/vectors,
%      one such set for each eigenvalue problem.
%
% Arguments: (Output)
%  Vseq - a 3-d array (pxpxn) of eigenvectors. Each
%      plane of the array will be sorted into a
%      consistent order with the other eigenvalue
%      problems. The ordering chosen will be one
%      that maximizes the energy of the consecutive
%      eigensystems relative to each other.
%
%  Dseq - pxn array of eigen values, sorted in order
%      to be consistent with each other and with the
%      eigenvectors in Vseq.
%
% Example:
%  Efun = @(t) [1 2*t+1 t^2 t^3;2*t+1 2-t t^2 1-t^3; ...
%               t