前言

在上一篇的PCA中使用了LinearAlgebra中封装的函数求取的特征值和特征向量。在此学习一下特征值的数值求解算法，主要参考[1]。矩阵特征值的数值求解方法中，一般有求部分特征值和特征向量的幂法和反幂法，以及求取全部特征值的QR分解方法。在此关注QR分解方法。注意内容较长，总体思路和方法来自参考文献[1]。

• QR分解，是把矩阵分解为正交矩阵和三角矩阵的乘积，其中关键点在于正交矩阵的求法。

• 正交矩阵的求解可以使用施密特正交化进行实现，需要理解施密特正交化的原理和思路。

最后回溯：由施密特正交化->>QR分解->>矩阵特征值->>特征向量矩阵 。

1. 施密特正交

施密特方法是对一组向量进行正交化。给定一组输入的m维向量，目的找出正交坐标系统，获取由这些向量张成的空间。
给定n个线性无关的向量，找到n个彼此垂直的单位向量。

(1) 利用施密特正交求出正交矩阵Q

假设存在nnn个线性无关的mmm维向量
A1,A2,⋯,AnA_{1},A_{2},\cdots,A_{n}A1​,A2​,⋯,An​

• 第111步：选定初始向量，归一化。

y1=A1y_{1}=A_{1}y1​=A1​

p1=y1∣∣y1∣∣2p_{1}=\dfrac{y_{1}}{||y_{1}||_{2}}p1​=∣∣y1​∣∣2​y1​​

• 第222步：下一个向量例如A2A_{2}A2​，减去该向量在所有已知的正交向量集(p1，pi,⋯)(p_{1}，p_{i},\cdots)(p1​，pi​,⋯)上的投影分量，再进行归一化得到当前的正交向量。

投影分量如何计算？利用向量的内积cos(θ)=<a,b>∣∣a∣∣2∣∣b∣∣2cos(\theta)=\dfrac{<a,b>}{||a||_{2}||b||_{2}}cos(θ)=∣∣a∣∣2​∣∣b∣∣2​<a,b>​，此处符号<a,b><a,b><a,b>为向量内积。

由于正交向量pip_{i}pi​的模长，也即是2范数为1，此时AjA_{j}Aj​在向量pip_{i}pi​上的投影系数为<Aj,pi>∣∣Aj∣∣2\dfrac{<A_{j},p_{i}>}{||A_{j}||_{2}}∣∣Aj​∣∣2​<Aj​,pi​>​,也即是piTAjp_{i}^{T}A_{j}piT​Aj​，注意此为一个数值标量。进一步可得到，AjA_{j}Aj​在向量pip_{i}pi​上的投影分量为pi(piTAj)p_{i}(p_{i}^{T}A_{j})pi​(piT​Aj​)。

y2=A2−p1(p1TA2)y_{2}=A_{2}-p_{1}(p_{1}^{T}A_{2})y2​=A2​−p1​(p1T​A2​)

p2=y2∣∣y2∣∣2p_{2}=\dfrac{y_{2}}{||y_{2}||_{2}}p2​=∣∣y2​∣∣2​y2​​

• 以此类推，第kkk步，

yk=Ak−p1(p1TAk)−⋯−pk−1(pk−1TAk)y_{k}=A_{k}-p_{1}(p_{1}^{T}A_{k})-\cdots -p_{k-1}(p_{k-1}^{T}A_{k})yk​=Ak​−p1​(p1T​Ak​)−⋯−pk−1​(pk−1T​Ak​)

pk=yk∣∣yk∣∣kp_{k}=\dfrac{y_{k}}{||y_{k}||_{k}}pk​=∣∣yk​∣∣k​yk​​

• 最终得出正交矩阵为Q=(p1，pi,⋯,pn)Q=(p_{1}，p_{i},\cdots,p_{n})Q=(p1​，pi​,⋯,pn​)，注意此时的维度为m∗nm*nm∗n，即mmm行nnn列。

(2) 求出上三角矩阵R

由以上求正交矩阵的过程，反推：
A1=p1∗∣∣y1∣∣2A_{1}=p_{1}*||y_{1}||_{2}A1​=p1​∗∣∣y1​∣∣2​

A2=p2∗∣∣y2∣∣2+p1∗(p1TA2)A_{2}=p_{2}*||y_{2}||_{2}+p_{1}*(p_{1}^{T}A_{2})A2​=p2​∗∣∣y2​∣∣2​+p1​∗(p1T​A2​)

Ak=pk∗∣∣yk∣∣2+p1∗(p1TAk)+⋯+pk−1∗(pk−1TAk)A_{k}=p_{k}*||y_{k}||_{2}+p_{1}*(p_{1}^{T}A_{k})+\cdots+p_{k-1}*(p_{k-1}^{T}A_{k})Ak​=pk​∗∣∣yk​∣∣2​+p1​∗(p1T​Ak​)+⋯+pk−1​∗(pk−1T​Ak​)

使用rijr_{ij}rij​代表矩阵R中的元素，令主对角线rii=∣∣yi∣∣2r_{ii}=||y_{i}||_{2}rii​=∣∣yi​∣∣2​，令rik=piTAk,i<=k−1r_{ik}=p_{i}^{T}A_{k},i<=k-1rik​=piT​Ak​,i<=k−1。则矩阵
A=(A1,Ai,⋯,An)=(p1,pi,⋯,pn)∗[r11r12⋯r1n0r22⋯r2n00⋯rin⋯⋯⋯⋯00⋯rnn]A=(A_{1},A_{i},\cdots,A_{n})=(p_{1},p_{i},\cdots,p_{n})* \begin{bmatrix} r_{11}& r_{12}&\cdots & r_{1n}\\ 0 & r_{22}&\cdots & r_{2n} \\ 0 & 0 &\cdots & r_{in} \\ \cdots & \cdots &\cdots & \cdots \\ 0 & 0 &\cdots & r_{nn} \\ \end{bmatrix}A=(A1​,Ai​,⋯,An​)=(p1​,pi​,⋯,pn​)∗⎣⎢⎢⎢⎢⎡​r11​00⋯0​r12​r22​0⋯0​⋯⋯⋯⋯⋯​r1n​r2n​rin​⋯rnn​​⎦⎥⎥⎥⎥⎤​

也即是 A=QRA=QRA=QR 此时的方法被称为消减QR分解。

#原始的消减QR分解,A_m*n=Q_m*n·R_n*n
function schmidit_QR(A::Array)m,n=size(A)Q=zeros(m,n)R=zeros(n,n)for j in 1:ny=A[:,j]if j!=1for i in 1:j-1R[i,j]=Q[:,i]'*A[:,j]y=y-R[i,j].*Q[:,i]endendR[j,j]=sqrt(sum(y.^2))Q[:,j]=y./R[j,j]endreturn Q,R
end

(3) 改进的消减QR分解

改进的消减QR分解的思路在于：改进矩阵R中的rijr_{ij}rij​。

此时计算的不再是当前向量AjA_{j}Aj​在向量pip_{i}pi​上的投影系数(piTAj)(p_{i}^{T}A_{j})(piT​Aj​)。而是，当前向量AjA_{j}Aj​减去之前所有在其他向量上的投影分量得到向量yyy，即是中间每一迭代的消减向量yyy，然后再对向量yyy在当前正交向量上求出投影系数作为rijr_{ij}rij​。

代码如下，仅在第10进行更改。

#改进的消减QR分解
function schmidit_QRc(A::Array)m,n=size(A)Q=zeros(m,n)R=zeros(n,n)for j in 1:ny=A[:,j]if j!=1for i in 1:j-1R[i,j]=Q[:,i]'*y #注意此处发生更改y=y-R[i,j].*Q[:,i]endendR[j,j]=sqrt(sum(y.^2))Q[:,j]=y./R[j,j]endreturn Q,R
end

测试
对比MATLAB中的消减qr函数qr(A,0)。

A=[1 -4; 2 3; 2 2]
Qc,Rc=schmidit_QRc(A)

2. 完全QR分解

以上的消减QR分解得出的是，Am∗n=Qm∗n∗Rn∗nA_{m*n}=Q_{m*n}*R_{n*n}Am∗n​=Qm∗n​∗Rn∗n​,当n<mn<mn<m时，矩阵Qm∗nQ_{m*n}Qm∗n​不是方阵。

完全QR分解即是将以上矩阵分解为Am∗n=Qm∗m∗Rm∗nA_{m*n}=Q_{m*m}*R_{m*n}Am∗n​=Qm∗m​∗Rm∗n​的形式，以更好的利用正交方阵的性质Qm∗mT=Qm∗m−1Q_{m*m}^{T}=Q_{m*m}^{-1}Qm∗mT​=Qm∗m−1​，即单位正交矩阵的转置即是其逆矩阵。

• 基于以下事实：
  (A1,Ai,⋯,An)=(p1,pi,⋯,pn,⋯,pm)∗[r11r12⋯r1n0r22⋯r2n00⋯rin⋯⋯⋯⋯00⋯rnn00⋯000⋯0](A_{1},A_{i},\cdots,A_{n})=(p_{1},p_{i},\cdots,p_{n},\cdots,p_{m})* \begin{bmatrix} r_{11}& r_{12}&\cdots & r_{1n}\\ 0 & r_{22}&\cdots & r_{2n} \\ 0 & 0 &\cdots & r_{in} \\ \cdots & \cdots &\cdots & \cdots \\ 0 & 0 &\cdots & r_{nn} \\ 0 & 0 &\cdots & 0 \\ 0 & 0 &\cdots & 0 \\ \end{bmatrix}(A1​,Ai​,⋯,An​)=(p1​,pi​,⋯,pn​,⋯,pm​)∗⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡​r11​00⋯000​r12​r22​0⋯000​⋯⋯⋯⋯⋯⋯⋯​r1n​r2n​rin​⋯rnn​00​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤​ 此时的RRR矩阵最后具有m−nm-nm−n行全0向量。
  根据以上的表达式，关键点是，在矩阵QQQ中如何求出额外的m−nm-nm−n个正交的单位列向量，而前n个列向量同消解QR分解得出的列向量一致。

• 解决思路是，在输入的矩阵(A1,Ai,⋯,An)(A_{1},A_{i},\cdots,A_{n})(A1​,Ai​,⋯,An​)后在补m−nm-nm−n个与之前向量都线性无关的列向量，把其补全为m∗mm*mm∗m的方阵。

• 如何确定向量是线性无关的。对于方针而言，求方阵的行列式，如果行列式不为零，则所有列向量线性无关。当然确定矩阵线性无关的方法，可以使用初等变换，或者利用最小二乘，在此使用行列式。

#矩阵的行列式的计算，利用递归和划分矩阵方法。
function mydet(A)m,n=size(A)if m!=nthrow("Error: the input Array must be a square matrix")elseif m==2return A[1,1]*A[2,2]-A[1,2]*A[2,1]elsesum=0for i in 1:m#对矩阵进行按列拆分del_row=A[setdiff(1:end,i),:] #集合操作，删除第i行 #经典cut_mat=del_row[:,setdiff(1:end,1)] #删除第1列sum += (-1)^((1+i)%2)*A[i,1]*mydet(cut_mat)endreturn sumendend
end

#将矩阵补全为方阵，并利用行列式 判断列向量之间是否线性无关
function repair_mat(A::Array)m,n=size(A)if m==nreturn Aelseres=zeros(m,m)while truerA=[A randn(m,abs(m-n))] #矩阵拼接if abs(mydet(rA))>eps(Float64)res=rAbreakendendreturn resend
end

#改进的完全QR分解,此时 A_m*n=Q_m*m·R_m*n
function schmidit_QRfull(B::Array)A=repair_mat(B)m,n=size(A)Q=zeros(m,n)R=zeros(n,n)for j in 1:ny=A[:,j]if j!=1for i in 1:j-1R[i,j]=Q[:,i]'*yy=y-R[i,j].*Q[:,i]endendR[j,j]=sqrt(sum(y.^2))Q[:,j]=y./R[j,j]endmb,nb=size(B)return Q,R[:,1:nb] #此时的R矩阵只需要前nb列即可
end

验证
MATLAB中使用qr(A)函数进行完全QR分解。

A=[1 -4; 2 3; 2 2]
Qf,Rf=schmidit_QRfull(A)

3. 矩阵QR分解的作用

(1) 方便求解线性方程组Ax=bAx=bAx=b

完全QR分解后：QRx=bQRx=bQRx=b -->> Rx=QTbRx=Q^{T}bRx=QTb。
由于矩阵RRR是上三角矩阵，此时可以回代得出所有的xix_{i}xi​。

(2) 求解最小二乘 min∣∣Ax−b∣∣2min ||Ax-b||_{2}min∣∣Ax−b∣∣2​

由于对于单位正交矩阵存在以下性质：
∣∣Qx∣∣2=∣∣x∣∣2||Qx||_{2}=||x||_{2}∣∣Qx∣∣2​=∣∣x∣∣2​，向量乘以单位正交矩阵模长不变。
∣∣Qx∣∣2=(Qx)T(Qx)=xTQTQx=xTx=∣∣x∣∣2||Qx||_{2}=(Qx)^{T}(Qx)=x^{T}Q^{T}Qx=x^{T}x=||x||_{2}∣∣Qx∣∣2​=(Qx)T(Qx)=xTQTQx=xTx=∣∣x∣∣2​

此时，
min∣∣Ax−b∣∣2=min∣∣QRx−b∣∣2=min∣∣Rx−QTb∣∣2min ||Ax-b||_{2}=min||QRx-b||_{2}=min||Rx-Q^{T}b||_{2}min∣∣Ax−b∣∣2​=min∣∣QRx−b∣∣2​=min∣∣Rx−QTb∣∣2​

注意到：
Rx−QTb=[r11r12⋯r1n0r22⋯r2n00⋯rin⋯⋯⋯⋯00⋯rnn00⋯000⋯0][x1⋯xn]−QT[b1⋯bn⋯bm]Rx-Q^{T}b= \begin{bmatrix}r_{11}& r_{12}&\cdots & r_{1n}\\ 0 & r_{22}&\cdots & r_{2n} \\ 0 & 0 &\cdots & r_{in} \\ \cdots & \cdots &\cdots & \cdots \\ 0 & 0 &\cdots & r_{nn} \\ 0 & 0 &\cdots & 0 \\ 0 & 0 &\cdots & 0 \\ \end{bmatrix} \begin{bmatrix}x_{1}\\ \cdots \\ x_{n}\end{bmatrix}-Q^{T} \begin{bmatrix}b_{1}\\ \cdots \\ b_{n}\\ \cdots\\ b_{m}\end{bmatrix}Rx−QTb=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡​r11​00⋯000​r12​r22​0⋯000​⋯⋯⋯⋯⋯⋯⋯​r1n​r2n​rin​⋯rnn​00​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤​⎣⎡​x1​⋯xn​​⎦⎤​−QT⎣⎢⎢⎢⎢⎡​b1​⋯bn​⋯bm​​⎦⎥⎥⎥⎥⎤​

由于Rx−QTbRx-Q^{T}bRx−QTb的前nnn项与自变量xix_{i}xi​有关，而后m−nm-nm−n项为确定性的常数项，与xix_{i}xi​无关。
令：
[d1⋮dm]=QT[b1⋮bm]\begin{bmatrix}d_{1}\\ \vdots\\ d_{m}\end{bmatrix}=Q^{T} \begin{bmatrix}b_{1}\\ \vdots \\ b_{m}\end{bmatrix}⎣⎢⎡​d1​⋮dm​​⎦⎥⎤​=QT⎣⎢⎡​b1​⋮bm​​⎦⎥⎤​

此时二范数最小化问题可以转化为前n项每一项都最小化的问题，也即是：
[r11r12⋯r1n0r22⋯r2n00⋯rin][x1⋮xn]−[d1⋮dn]=0\begin{bmatrix}r_{11}& r_{12}&\cdots & r_{1n}\\ 0 & r_{22}&\cdots & r_{2n} \\ 0 & 0 &\cdots & r_{in} \end{bmatrix} \begin{bmatrix}x_{1}\\ \vdots \\ x_{n}\end{bmatrix}-\begin{bmatrix}d_{1}\\ \vdots \\ d_{n}\end{bmatrix}=0⎣⎡​r11​00​r12​r22​0​⋯⋯⋯​r1n​r2n​rin​​⎦⎤​⎣⎢⎡​x1​⋮xn​​⎦⎥⎤​−⎣⎢⎡​d1​⋮dn​​⎦⎥⎤​=0

回代以上线性方程组，可以使得did_{i}di​的前nnn项为0，此时的最小二乘误差为：
∣∣e∣∣22=dn+12+⋯+dm2||e||_{2}^{2}=d_{n+1}^{2}+\cdots+d_{m}^{2}∣∣e∣∣22​=dn+12​+⋯+dm2​

(3) 求解特征值

矩阵特征值的数值求解方法中，一般有求部分特征值和特征向量的幂法和反幂法，以及求取全部特征值的QR分解方法。在此关注QR分解方法。

对于矩阵特征值的求解，QR分解中利用相似矩阵的性质。

• 如果矩阵AAA和A0A_{0}A0​相似（即存在可逆矩阵PPP，使得A=PA0P−1A=PA_{0}P^{-1}A=PA0​P−1,则矩阵AAA和A0A_{0}A0​相似），则相似矩阵具有相同的特征值。

非奇异矩阵特征值的性质。

• 对于nnn阶的实数非奇异矩阵，存在nnn个特征值。但是特征值中可能存在重根特征值，或者成对的复数特征值。

(3.1) 平移QR算法求解矩阵特征值

定理：对于实数元素的方阵AAA，存在正交矩阵QQQ和实数舒尔形式的矩阵TTT，满足A=QTTQA=Q^{T}TQA=QTTQ。其中，实数舒尔矩阵TTT是主对角线存在2 x 2的矩阵块的矩阵。例如：
[xxxxxxxxxxxxxxxxx]\begin{bmatrix} x & x & x & x & x \\ x & x & x & x & x\\ & & x & x & x\\ & & x & x & x\\ & & & & x \end{bmatrix}⎣⎢⎢⎢⎢⎡​xx​xx​xxxx​xxxx​xxxxx​⎦⎥⎥⎥⎥⎤​

这种形式的矩阵特征值是其对角块矩阵的特征值，即是对应的对角线元素（其为1x1的块）或者2x2块矩阵的特征值。矩阵AAA的舒尔分解也被称为“特征值发现分解”，如果能够实现这个分解，则能够找到相应的特征值[1]。

完整的QR算法通过一系列相似变换迭代，将任意矩阵移动到它的舒尔分解形式。

平移QR算法原理：

• 对于非奇异方阵AAA，定义A0=AA_{0}=AA0​=A，
  A0−sI=Q1R1A_{0}-sI=Q_{1}R_{1}A0​−sI=Q1​R1​ A1=R1Q1+sIA_{1}=R_{1}Q_{1}+sIA1​=R1​Q1​+sI
  其中III为单位矩阵，sss为标量系数。由于
  A1−sI=R1Q1=Q1T(A0−sI)Q1=Q1TA0Q1−sIA_{1}-sI=R_{1}Q_{1}=Q_{1}^{T}(A_{0}-sI)Q_{1}=Q_{1}^{T}A_{0}Q_{1}-sIA1​−sI=R1​Q1​=Q1T​(A0​−sI)Q1​=Q1T​A0​Q1​−sI因此A0和A1A_{0}和A_{1}A0​和A1​是相似矩阵，因而，Ak和A1A_{k}和A_{1}Ak​和A1​也是相似矩阵。

步骤

• 在每一步中使用QR步骤，将任意矩阵移动到它的实数舒尔形式，然后检查最后一行。如果最后一行除了对角线元素anna_{nn}ann​之外所有元素都很小，则这个值就是特征值。在余下的计算中去掉最后一行和最后一列进行收缩，然后在迭代寻找其他特征值 [1]。

#eyes(n)产生n阶单位矩阵，eyes(n,m)产生对角线为1的矩阵
function eyes(m::Int,index::Int=0)n=0if index==0n=melsen=indexendI=zeros(m,n)[I[i,i]=1 for i in 1:min(m,n)]return I
end# 平移QR算法求特征值，原理：逆向幂迭代和收缩技术推出平移QR
# 仅计算方阵的实数特征值，并且A中没有大小相同的特征值
# lam为返回的特征值
function shiftQR(A::Array)tol=1e-14 #精度门限m=size(A,1)lam=zeros(m,1)n=mwhile n>1#abs()函数针对标量，参数为矢量时需要在函数后面加点号while sort(abs.(A[n,1:n-1]))[end]>tol #sort(abs(A[n,1:n-1]))[end] 与max(abs(A[n,1:n-1]))等价mu=A[n,n]q,r=schmidit_QRfull(A-mu*eyes(n))A=r*q+mu*eyes(n)endlam[n]=A[n,n]n=n-1A=A[1:n,1:n]endlam[1]=A[1,1] #存在错误return lam
end#平移QR算法，能够计算方阵的实数和复数特征值
function shiftQRc(A::Array) #A为方阵tol=1e-14maxiter=500m=size(A,1)lam=zeros(m,1)n=mwhile n>1iter=0while sort(abs.(A[n,1:n-1]))[end]>tol && iter<maxiteriter=iter+1 #记录qr的次数mu=A[n,n] #平移量q,r=schmidit_QRfull(A-mu*eyes(n))A=r*q+mu*eyes(n)endif iter<maxiter #具有1x1的块lam[n]=A[n,n]n=n-1A=A[1:n,1:n]else#迭代超过最大迭代次数，认为是2x2的块#lam为实数向量，需要转化为复数向量，否则后面报错，复数不能直接复制给实数变量lam=complex(lam)disc=(A[n-1,n-1]-A[n,n])^2+4*A[n,n-1]*A[n-1,n] #一元二次求根公式 a^2-4actemp=sqrt.(disc+0*im)#需要添加0*im变为复数，才能开方lam[n]=(A[n-1,n-1]+A[n,n] +temp)/2lam[n-1]=(A[n-1,n-1]+A[n,n]-temp)/2n=n-2 #以2收缩A=A[1:n,1:n]endendif n>0lam[1]=A[1,1]endreturn lam
end

测试

(3.2) 以上平移QR分解存在的问题

以上的平移QR分解算法能够求解大多数方阵的特征值，但是对于具有重数的复数特征值，以上的算法无法将原始矩阵移动为舒尔形式。例如对于矩阵A=[0 0 0 1;0 0 -1 0;0 1 0 0;-1 0 0 0]，在MATLAB环境下，求得的特征值为[0+1∗im;0−1∗im;0+1∗im;0−1∗im][0+1*im;0-1*im;0+1*im;0-1*im][0+1∗im;0−1∗im;0+1∗im;0−1∗im]，而以上方法求出的特征值全为0。此时需要，首先对原始矩阵进行预处理，然后再使用以上算法求解。

(3.3) 解决方案

上海森伯格矩阵
在矩阵中，如果元素aij=0，i>j+1a_{ij}=0，i>j+1aij​=0，i>j+1,则m∗nm*nm∗n矩阵是上海森伯格形式。如下示例：
[xxxxxxxxxxxxx]\begin{bmatrix} x & x & x & x \\ x & x & x & x \\ & x & x & x \\ & & x & x \\ \end{bmatrix}⎣⎢⎢⎡​xx​xxx​xxxx​xxxx​⎦⎥⎥⎤​

需要注意的是，问题关键在于如何将原始矩阵转化为上海森伯格形式。
在此，使用矩阵QR分解中的House Holder反射子，将原始矩阵的每一列转化为只有第一个元素非零，其余都为零的列向量。由于我们想要的是实数舒尔形式，例如矩阵第一列，此时需要保留第一列的前两个元素，这时初始选择的向量，应当截取第一列的第二个元素到最后作为最初映射的原始向量，例如A[2:end,1]A[2:end,1]A[2:end,1]，此时第一个元素a11a_{11}a11​不参与运算。

投影矩阵和House Holder反射子
基本思想：对于任意向量xxx，我们期望把它转变成x轴上的向量，例如w=(∣∣x∣∣2,0,⋯,0)Tw=(||x||_{2},0,\cdots,0)^{T}w=(∣∣x∣∣2​,0,⋯,0)T，即仅在第一个坐标位置非零，并且模长与向量xxx相等。

假设存在矩阵HHH，使得Hx=wHx=wHx=w,向量www如上所示，此时的矩阵HHH为转换矩阵。

考虑等腰三角形，向量x,wx,wx,w可以看做是三角形的腰长，根据等腰三角形性质，其底边向量与垂线向量相互垂直，此时的底边向量可以表示为w−xw-xw−x,垂线向量表示为w+xw+xw+x,二者为正交向量。

投影矩阵
对于任意非零向量，定义矩阵PPP如下所示:
P=vvTvTvP=\dfrac{vv^{T}}{v^{T}v}P=vTvvvT​
此时的矩阵PPP是关于向量vvv的投影矩阵。

矩阵PPP具有如下性质：
P2=vvTvvT(vTv)(vTv)=vvTvTv=PP^{2}=\dfrac{vv^{T}vv^{T}}{(v^{T}v)(v^{T}v)}=\dfrac{vv^{T}}{v^{T}v}=PP2=(vTv)(vTv)vvTvvT​=vTvvvT​=P
对于向量三角形，向量xxx在底边v=w−xv=w-xv=w−x的投影为PxPxPx。根据几何关系，向量xxx减去其2倍的投影得到的向量就是向量www，即
x−2Px=(I−2P)x=wx-2Px=(I-2P)x=wx−2Px=(I−2P)x=w

因此，可以得出变换矩阵H=I−2PH=I-2PH=I−2P，此时的HHH被称为House-Holder反射子，具有以下性质。

由于矩阵HHH是对称正交矩阵，
HTH=HH=(I−2P)(I−2P)=I−4P+4P2=IH^{T}H=HH=(I-2P)(I-2P)=I-4P+4P^{2}=IHTH=HH=(I−2P)(I−2P)=I−4P+4P2=I此时可以利用House-Holder反射子进行矩阵的QR分解。

根据文献[1]478页，以5阶方阵A为例。矩阵AAA左乘House-Holder算子H1H_{1}H1​为H1AH_{1}AH1​A，其第一列的后三个元素全为0，此时的H1AH_{1}AH1​A与AAA不再相似，需要在其右侧右乘矩阵H1−1=H1TH_{1}^{-1}=H_{1}^{T}H1−1​=H1T​保持相似性，即H1AH1TH_{1}AH_{1}^{T}H1​AH1T​。依次迭代（3=5阶-2）次，即是
H3H2H1AH1TH2TH3T=(H3H2H1)A(H3H2H1)TH_{3}H_{2}H_{1}AH_{1}^{T} H_{2}^{T}H_{3}^{T}=(H_{3}H_{2}H_{1})A(H_{3} H_{2}H_{1})^{T}H3​H2​H1​AH1T​H2T​H3T​=(H3​H2​H1​)A(H3​H2​H1​)T 以上形式的迭代矩阵同原始矩阵AAA相似。
注意，一般对于一个nnn阶方阵AAA，需要(n−2)(n-2)(n−2)个House-Holder算子将其转化为上海森伯格形式。

实现代码

#将矩阵转化为上海森伯格形式
function hessen(A::Array)m,n=size(A)v=zeros(m,m)for k in 1:m-2x=A[k+1:m,k]v[1:m-k,k]=-1*sign(x[1]+eps())*sqrt(sum(x.^2))*eyes(m-k,1)-xv[1:m-k,k]=v[1:m-k,k]./sqrt(sum(v[1:m-k,k].^2)) #归一化A[k+1:m,k:m]=A[k+1:m,k:m]-2*v[1:m-k,k]*v[1:m-k,k]'*A[k+1:m,k:m] #投影求出正交向量A[1:m,k+1:m]=A[1:m,k+1:m]-2*A[:,k+1:m]*v[1:m-k,k]*v[1:m-k,k]' #在对剩余的矩阵右乘House-Holder矩阵，保持矩阵相似性endreturn A
end

感兴趣的同学，建议读一下参考文献，具有详细步骤示意图，更加直观。

测试

需要注意的是，对于编写的hessen()函数在MATLAB环境下对测试矩阵A进行运算结果正确。然而，在Julia环境下，对于整数矩阵直接求解运算，报错，但是在加上精度数据eps()之后，再运算结果正确，目前还没有搞明白原因。注意矩阵和标量的运算，需要使用点号进行广播操作。

总结

总体而言，数值方法求解矩阵特征值，使用的都是迭代思想和矩阵的相似性质。House-Holder投影算子，启发于等腰三角形，十分巧妙。
写到现在，使用QR分解，只求解出了矩阵的全部特征值，距离相应的特征向量仅有一步之遥。累了，后续休息休息，再进行更新和完善。
参考文献
[1] 裴玉茹，马庚宇译，数值分析(第二版)，机械工业出版社。

Julia 矩阵QR分解和特征值相关推荐

1. Python+numpy实现矩阵QR分解

   感谢广东东软学院计算机系赵晨杰老师的交流. 如果实(复)非奇异矩阵A能够化成正交(酉)矩阵Q与实(复)非奇异上三角矩阵R的乘积,即A=QR,则称其为A的QR分解. Python扩展库numpy实现了矩 ...

2. 双步位移求解特征值matlab,数值分析——带双步位移的QR分解求特征值算法

   C语言实现数值分析中带双步位移的QR分解求特征值算法. 数 值 分 析(B) 大 作 业(二) 1.算法设计: ①矩阵的拟上三角化: 对实矩阵A进行相似变换化为拟上三角矩阵A(n 1),其变换矩阵采用 ...

3. 【NA】基于QR分解的特征值迭代法

   Francis于1961-1962年利用矩阵的QR分解建立了计算矩阵特征值的QR方法,是计算中小型矩阵全部特征值的最有效方法之一. 本篇的主线是第一部分介绍QR分解,第二部分介绍从QR分解引出的特征值 ...

4. r语言中矩阵QR分解_R语言常用的矩阵操作

   R语言是一门非常方便的数据分析语言,它内置了许多处理矩阵的方法.下面列出一些常用的矩阵操作方法示例. 矩阵的生成 > mat <- matrix(1:16, ncol = 4, nrow ...

5. 数值分析--矩阵QR分解的三种方法

   QR分解法是目前求一般矩阵全部特征值的最有效并广泛应用的方法,一般矩阵先经过正交相似变化成为Hessenberg矩阵,然后再应用QR方法求特征值和特征向量.它是将矩阵分解成一个正规正交矩阵Q与上三角形 ...

6. QR算法的Matlab 程序,三种实现矩阵QR分解的算法与程序

   To learn, to share, to debate, then comes progress. ------------------------------------------------ ...

7. 用Givens旋转进行矩阵QR分解

   不多废话,直接贴代码 function [A_, T_, T] = my_qr_givens( A ) %利用givens旋转进行qr分解 %输出 %A_ 每次变换后的A矩阵 %T_ 对应于A_的变换 ...

8. r语言中矩阵QR分解_从零开始学R语言Day4|向量、矩阵和数组

   从零开始学R语言Day4|向量.矩阵和数组 1.1向量 1.1.1向量 在Day2中我们提及过用和c()函数来构建向量,具体实例如下. 我们还可以采用vector("类型",长度) ...

9. 施密特正交化 对矩阵QR分解

最新文章

1. 给GRUB添加新的项目
2. .net网格怎么把值插入指定列表_Python列表有什么内置函数可以使用，怎么使用这些函数...
3. 递归二分法php,PHP基于二分法实现数组查找功能示例【循环与递归算法】
4. excel 避免下拉操作
5. java基础之线程（1）
6. uwsgi怎么通过浏览器访问某个脚本_4个Shell小技巧帮你提高机器学习效率：写好脚本，事半功倍...
7. 苹果隐藏app_iOS 14的隐藏功能盘点：不知道等于白更新！
8. 剑指offer之对称的二叉树
9. 知识蒸馏 | 综述: 网络结构搜索应用
10. 联想服务器改win7系统教程,联想台式机10代cpu改win7系统详细教程
11. 计算机网络基础知识 - 物理层
12. vue学习笔记 el-dialog 固定宽度
13. string entitlement = Application.dataPath+ “/Editor/Entitle Unity工程到处iOS工程，用脚本把Push Notifications打开
14. 如何从官网下载Chrome浏览器离线安装包
15. QT绘图实现橡皮擦效果
16. java里arcsin_java编程用泰勒级数计算arcsin
17. 互联网”还是“技术”？派系分明的电子烟市场
18. 高德 php,高德地图WEB版的使用
19. 网易视频云 php接口
20. 消息队列RabbitMQ基本使用(Java代码实现)

热门文章

1. 上海亚商投顾：沪指震荡反弹涨1.2% 中国移动创历史新高
2. 《如何让你爱的人爱上你》-读书笔记
3. 免费收录外链的分类目录网站推荐
4. 【Gym - 102174J】 金色传说（观察性质+计数dp）
5. 利用jieba库和wordcloud库，进行中文词频统计并利用词云图进行数据可视化
6. Spring Security 4 整合Hibernate 实现持久化登录验证（带源码）
7. nextcloud+docker在阿里云服务器上搭建个人云存储盘（如何在服务器上搭建私有云盘）
8. 【果断收藏】阿里云盘首页
9. 多重签名地址和P2SH
10. 推荐一个程序员业余写的小说