Hermitian矩阵特征分解学习笔记

1. Hermitian 矩阵实对称化

​ Hermitian 矩阵是复数方阵，$R=(r_{ij})$满足$R=R^H$，其中$R^H$表示共轭转置，相当于$r_{ij}=r_{ji}^{\ast }$ ，例如对于一个四阶矩阵应满足：
$$R=\left[\begin{array}{cccc}{r}_{11}& r_{21}^{\ast }& r_{31}^{\ast }& r_{41}^{\ast }\\ r_{21}& r_{22}& r_{32}^{\ast }& r_{31}^{\ast }\\ r_{31}& r_{32}& r_{22}& r_{21}^{\ast }\\ r_{41}& r_{31}& r_{21}& r_{11}\end{array}\right]$$
​ 在针对此复数矩阵的特征分解计算中可以将其转化为实对称矩阵再进行计算
$$\mathbf{RX}=\lambda \mathbf{X}$$
其中 $\mathbf{X}$ 是协方差矩阵的特征向量，$\lambda$ 是协方差矩阵的特征值。将 $\mathbf{R}$ 和 $\mathbf{X}$ 以实部和虚部的形式表示：
$$\begin{gathered} \mathbf{R}={\mathbf{R}}_{\mathbf{r}}+i{\mathbf{R}}_{\mathbf{i}} \\ \mathbf{X}=\mathbf{U}+i\mathbf{V} \end{gathered}$$
将上式表达成矩阵形式为：
$$\left[\begin{array}{cc}{\mathbf{R}}_r& -{\mathbf{R}}_i\\ {\mathbf{R}}_i& {\mathbf{R}}_r\end{array}\right]\left[\begin{array}{c}\mathbf{U}\\ \mathbf{V}\end{array}\right]=\lambda \left[\begin{array}{c}\mathbf{U}\\ \mathbf{V}\end{array}\right]$$
将 $M\times M$ 阶的复 Hermitian 矩阵 $\mathbf{R}$ 的特征值分解等价于求解由实矩阵 ${\mathbf{R}}_r$、${\mathbf{R}}_i$ 组成的 $2M\times 2M$ 阶分块矩阵 ${\mathbf{R}}^{\prime }$ 的特征值分解。

Matlab仿真验证

1. 首先构建Hermitian 矩阵

图1. 构建Hermitian矩阵

2. 构建实对称矩阵

图2. 构建实对称矩阵

3. 对实对称矩阵进行Jacobi特征分解

图3. 实对称矩阵Jacobi特征分解

function [V,D,iter]=Jacobi_classical(A,maxIter,tol)
n = size(A, 1); % 矩阵的大小
V = eye(n); % 初始化特征向量矩阵为单位矩阵
iter = 0; % 初始化迭代次数
% 设置最大迭代次数和误差精度
if nargin < 3 || isempty(tol)
    tol = 1e-9; % 默认误差精度
end
if nargin < 2 || isempty(maxIter)
    maxIter = 1000; % 默认最大迭代次数
end
while(iter < maxIter)
    iter=iter+1;
    D=A;
    n=size(D,1);
    p=1;q=2;
    for i=1:n
        for j=i+1:n
            if(abs(D(i,j))>abs(D(p,q)))%找到对称矩阵的上三角矩阵中最大的元素的下标
                p=i;q=j;
            end
        end
    end
    if(abs(D(p,q))<tol)
        break;
    end
    if(A(p,q)~=0)
        d=(A(q,q)-A(p,p))/(2*A(p,q));
        if(d>0)
            t=1/(d+sqrt(d^2+1));
        else
            t=-1/(-d+sqrt(d^2+1));
        end
        c=1/sqrt(t^2+1);s=c*t;
    else
        c=1;s=0;
    end
    R=[c s;-s c];
    A([p,q],:)=R'*A([p,q],:);
    A(:,[p,q])=A(:,[p,q])*R;
    V(:, [p, q]) = V(:,[p,q])*R;
end

4. 观察特征值对角矩阵可得，每个特征值都有2重，所以取一半特征值，同时根据

$$\mathbf{X}=\left[\begin{array}{c}\mathbf{U}\\ \mathbf{V}\end{array}\right]$$

重新组合特征向量，可得

图4. 重新组合特征值、特征向量

R=[real(Rxx),-imag(Rxx);imag(Rxx),real(Rxx);];%化为实对称矩阵
[VV,DD,iter]=Jocobi_classical(R);%求特征值、特征向量
EV=VV(1:M, 1:M)+1j*VV(M+1:2*M,1:M);D=DD(1:M, 1:M);%重新组合为复数特征向量X

5. 使用eig函数进行比对

图5. eig函数结果

由于特征向量不唯一，可根据定义$\mathbf{RX}=\lambda \mathbf{X}$验证，验证代码如下

function verificationResult = verifyEigen(A, D, V)
% verifyEigen 函数验证给定的矩阵 A 的特征值和特征向量
% 输入：
%   A - 原始矩阵
%   D - 特征值矩阵，对角线上是特征值，其他元素为0
%   V - 特征向量矩阵，每一列是对应的特征向量
% 输出：
%   verificationResult - 逻辑值，true 如果所有特征值和特征向量验证通过，否则为 false
n = size(A, 1); % 获取矩阵 A 的大小
verificationResult = true; % 初始化验证结果为 true
% 检查 V 的列是否为 A 的特征向量
for i = 1:n
    % 取出第 i 个特征值和对应的特征向量
    lambda = D(i,i);
    eigenvector = V(:,i);
    % 计算 A * eigenvector
    Av = A * eigenvector;
    % 计算 lambda * eigenvector
    lambda_v = lambda * eigenvector;
    % 检查 A * eigenvector 是否等于 lambda * eigenvector
    % 使用 norm 函数和 inf 范数来检查两个向量的差是否在容忍误差内
    tolerance = 1e-10; % 设置一个容忍误差
    if norm(Av - lambda_v, inf) > tolerance
        verificationResult = false; % 如果不等，则验证失败
        fprintf('验证失败：特征值 %f 的特征向量\n', lambda);
        break; % 跳出循环
    end
end
% 输出最终验证结果
if verificationResult
    disp('所有特征值和特征向量验证成功。');
else
    disp('存在不满足特征值和特征向量条件的案例。');
end
end

图6. 特征向量验证

2. 二阶主子阵实数化法与Jacobi方法

​ 上一种方法将 $M\times M$ 阶的复 Hermitian 矩阵 $\mathbf{R}$ 的特征值分解等价于求解由实矩阵 ${\mathbf{R}}_r$、${\mathbf{R}}_i$ 组成的 $2M\times 2M$ 阶分块矩阵 ${\mathbf{R}}^{\prime }$ 的特征值分解，提高了阶数也就提高了计算时间，以下方法是变换过程中将元素实数化再应用Jacobi方法求解¹。

​ 将二阶主子阵实数化法²与 Jacobi 方法相结合，算法的主要思想是在每个迭代中，将矩阵的一个二阶主子阵用酉对角阵相似变换成二阶的实对称矩阵，然后再利用 Jacobi 方法将此二阶实对称矩阵对角化。通过不断循环，不断选取不同的二阶主子阵，应用 Jacobi 方法进行对角化，直到整个矩阵最终近似为对角阵。

首先用酉对角阵相似变换变换成二阶的实对称矩阵，设在 Hermitian 矩阵 $\mathbf{A}=\mathbf{B}+i\mathbf{C}$ 中, $a_{pq}=b_{pq}+i{c}_{pq}\ne 0,(p<q)$ , 则接下来讨论二阶主子阵的西对角化过程,

$${\mathbf{A}}_{pq}=\left[\begin{array}{ll}{a}_{pp}& a_{pq}\\ a_{qp}& a_{qq}\end{array}\right]=\left[\begin{array}{ll}{b}_{pp}& b_{pq}\\ b_{qp}& b_{qq}\end{array}\right]+i\left[\begin{array}{ll}{c}_{pp}& c_{pq}\\ c_{qp}& c_{qq}\end{array}\right]={\mathbf{B}}_{pq}+i{\mathbf{C}}_{pq}\text{.}$$

根据反对角线上的元素作西对角阵, 有

$${\mathbf{T}}_{pq}=\left[\begin{array}{cc}1& 0\\ 0& \frac{{\bar{a}}_{pq}}{|a_{pq}|}\end{array}\right]=\left[\begin{array}{cc}1& 0\\ 0& \frac{b_{pq}-i{c}_{pq}}{\sqrt{b_{pq}^2+c_{pq}^2}}\end{array}\right]=\left[\begin{array}{cc}1& 0\\ 0& {\mathrm{e}}^{-j\phi }\end{array}\right].$$

在上式 中, $\phi =\mathrm{angle}\left(a_{pq}\right)$ , 其中 ${\mathrm{e}}^{-j\phi }=\cos \phi -j\sin \phi$ , 故有

$$\cos \phi =\frac{b_{pq}}{\sqrt{b_{pq}^2+c_{pq}^2}},\sin \phi =\frac{c_{pq}}{\sqrt{b_{pq}^2+c_{pq}^2}}.$$

因而通过西对角化, 有

$${\mathbf{A}}_{pq}^{(1)}={\mathbf{T}}_{pq}^H{\mathbf{A}}_{pq}{\mathbf{T}}_{pq}=\left[\begin{array}{cc}{a}_{pp}& \left|a_{pq}\right|\\ \left|a_{pq}\right|& a_{qq}\end{array}\right]={\mathbf{B}}_{pq}^{(1)}.$$
通过得到的 ${\mathbf{B}}_{pq}^{(1)}$ 为二阶实对称矩阵. 然后再将其通过 Jacobi 旋转变换之后得到对角阵. 在本特征分解算法中, 使用的平面旋转矩阵的模型为

$${\mathbf{J}}_{pq}=\left[\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right],$$

则有

$${\mathbf{B}}_{pq}^{(2)}={\mathbf{J}}_{pq}^T{\mathbf{B}}_{pq}^{(1)}{\mathbf{J}}_{pq}=\left[\begin{array}{cc}{a}_{pp}^{(2)}& 0\\ 0& a_{qq}^{(2)}\end{array}\right].$$

式 中, ${\mathbf{B}}_{pq}^{(2)}$ 为对角阵.通过分析可得，${\mathbf{T}}_{pq}^H\mathbf{A}{\mathbf{T}}_{pq}$ 会改变 $\mathbf{A}$ 的 $q$ 行 $q$ 列
$${\mathrm{e}}^{-j\phi }\ast {\left(\mathbf{A}\right)}_{kq}\to \left({\mathbf{A}}^{\prime }\right)$$

$${\mathrm{e}}^{j\phi }\ast {\left({\mathbf{A}}^{\prime }\right)}_{qk}\to \left({\mathbf{A}}^{\prime \prime }\right)$$

这样 ${\mathbf{A}}^{\prime \prime }$ 就部分实数化了再进行后续的 Givens 变换，其中旋转矩阵中的 $\theta$ 满足下式：
$$\tan 2\theta =\frac{2\tan \theta }{1-{\tan }^2\theta }=\frac{2\mid a_{pq}^{(k)}\mid }{a_{pp}^{(k)}-a_{qq}^{(k)}}.$$
通过不断迭代就可以对角化得到特征值，记录迭代过程的旋转矩阵就可以得到特征向量

实现代码如下：

function [V,D,iter]=Jacobi_Hermitian(A,maxIter,tol)
%厄米特矩阵实数化后jacobi迭代求解
n = size(A, 1); % 矩阵的大小
V = eye(n); % 初始化特征向量矩阵为单位矩阵
iter = 0; % 初始化迭代次数
% 设置最大迭代次数和误差精度
if nargin < 3 || isempty(tol)
    tol = 1e-9; % 默认误差精度
end
if nargin < 2 || isempty(maxIter)
    maxIter = 1000; % 默认最大迭代次数
end
while(iter < maxIter)
    iter=iter+1;
    D=A;
    p=1;q=2;
    for i=1:n
        for j=i+1:n
            if(abs(A(i,j))>abs(A(p,q)))%找到对称矩阵的上三角矩阵中最大的元素的下标
                p=i;q=j;
            end
        end
    end
    if(abs(D(p,q))<tol)
        break;
    end
    phi=angle(A(p,q));%实数化角度
    A(:,q)=exp(-1i*phi).*A(:,q);%实数化涉及到的列变换
    A(q,:)=exp(1i*phi).*A(q,:);%实数化涉及到的行变换
    V(:,q)=exp(-1i*phi).*V(:,q);%记录变换矩阵
    theta=0.5*atan(2*abs(A(p,q))/(A(p,p)-A(q,q)));%Givens变换角度
    c=cos(theta);s=sin(theta);
    R=[c s;-s c];
    A([p,q],:)=R*A([p,q],:);%Givens变换涉及到的行变换
    A(:,[p,q])=A(:,[p,q])*R';%Givens变换涉及到的列变换
    V(:,[p,q])=V(:,[p,q])*R';%记录变换矩阵
end

Matlab验证

图7. 法2 Matlab验证

3. 参考文献

---

1. 曾富红,司伟建,曲志昱.基于Hermitian矩阵的特征分解算法[J].沈阳大学学报(自然科学版),2016,28(06):512-517. ↩︎

2. 征道生.Hermite矩阵特征值问题的2阶主子阵实数化法[J].华东师范大学学报(自然科学版),1996(03):1-6. ↩︎