Householder和Lanczos双对角化算法用R语言实现

目录

目录

基础知识

Householder变换
H变换双对角化代码实现

未完待续

基础知识

复数：z=a+bi，其中a，b是实数，i是虚数单位。

a=Re(z)称为实部，实部为零时，z为纯虚数。

b=Im(z)称为虚部，虚部b等于零时，z可以视为实数。

共轭复数：实部相等，虚部互为相反数的复数。

例如：x+yi与x-yi称为共轭复数。

共轭复数有些有趣的性质：

|x+yi|=|x−yi| (x+yi)×(x−yi)=x2+y2

共轭转置：把行和列互换再共轭一下，记为AH

一般A上面加个横线 表示复共轭（共轭复数）。有时候共轭转置用*符号。

性质：

(AB)∗=B∗A∗ (A∗)∗=A

Hermite阵：也称共轭矩阵，与自己共轭转置。A=AH （实数范围称对称矩阵, A=AT）

酋矩阵：实数域上称为正交矩阵 AHA=AAH=E。

正规矩阵: 任意正规矩阵都可在经过一个酉变换后变为对角矩阵，反过来所有可在经过一个酉变换后变为对角矩阵的矩阵都是正规矩阵。

逆矩阵：我们称B是A的逆矩阵，AB=BA=E

奇异矩阵：行列式=0的矩阵，也就是rank(A)<n。若A为非奇异矩阵，则A满秩，rank(A)=n。

Householder变换

Householder变换又称为反射变换或镜像变换，有明显的几何意义。

在R3中，给定一个向量α，令β表示α关于平面π（以ω为法向量）的反射变换所得得像，如图所示，记：

ω=α−β|α−β|∈R3H(ω)=I−2ωωT

则

H(ω)α=β



即：该变换将向量α变成了以ω为法向量的平面的对称向量β。

H变换双对角化代码实现

Householder<-function(A)
{
#列变换
H1<-function(A,n)
{
if(n>=nrow(A)){return (diag(1,nrow =nrow(A),ncol = nrow(A)))}
a<-A[n:nrow(A),n];m<-length(a)
I<-diag(x=1L,nrow = m,ncol = m)
e<-numeric(length=m);e[1]<-1
w<-a-norm(as.matrix(a),type = "f")*e
w<-w/norm(as.matrix(w),type = "f")
h<-I-2*w%*%t(w)
h<-as.matrix(bdiag(diag(n-1),h))
h<-round(h,digits = 15)
return(h)
}
#行变换
H2<-function(A,n)
{
if((n+1)>=ncol(A)){return (diag(1,nrow =nrow(A),ncol = nrow(A)))}
a<-A[n,(n+1):ncol(A)];m<-length(a)
e<-numeric(length=m); e[1]<-1
w<-a-norm(as.matrix(a),type = "f")*e
w<-w/norm(as.matrix(w),type = "f")
I<-diag(x=1L,nrow = m,ncol = m)
h<-I-2*w%*%t(w)
h<-as.matrix(bdiag(diag(n),h))
h<-round(h,digits = 15)
return(h)
}
#主函数
a<-min(nrow(A),ncol(A))
for(n in 1:a)
{
A<-H1(A,n)%*%A
A<-A%*%H2(A,n)
}
A<-round(A,digits = 10)#四舍五入
return(A)
}

未完待续