数値解析 共立出版 森正武先生


#QR法 p.94  (森正武著 共立出版「数値解析」)
#固有値(A_matの対角成分)と固有ベクトル(Q:Aが対象行列の場合)
A=matrix(sample(1:9,9,replace=T),ncol=3)
A=A+t(A)
A_mat=A
ite=1000
Q=diag(rep(1,nrow(A_mat)))
for(l in 1:ite){
q_mat=array(0,dim=c(nrow(A_mat),ncol(A_mat)))
r_mat=array(0,dim=c(nrow(A_mat),ncol(A_mat)))
q1=A_mat[,1]/sqrt(sum(A_mat[,1]^2))
q_mat[,1]=q1
r_mat[1,1]=sqrt(sum(A_mat[,1]^2))
for(j in 2:ncol(A_mat)){
u=A_mat[,j]-apply(t(t(q_mat[,1:(j-1)])*c(t(q_mat[,1:(j-1)])%*%A_mat[,j])),1,sum)
r_mat[j,j]=sqrt(sum(u^2))
q_mat[,j]=u/sqrt(sum(u^2))
}
r_mat2=t(A_mat)%*%q_mat
r_mat2[upper.tri(r_mat2)]=0
diag(r_mat2)=diag(r_mat)
r_mat=t(r_mat2)
A_mat=r_mat%*%q_mat
Q=Q%*%q_mat
}