- こちらで正方行列のQR分解とは『正方行列を回転行列Qと上三角行列Rとに分解すること』であって、『上三角行列Rの列ベクトルについては、個々の列ベクトルのノルムと、列ベクトルのペアのなす角とが与えられたときに、作り上げられるもの』という『意味』を書いた
- 逆に言えば、任意の正方行列の、列ベクトルについて、そのノルムと、列ベクトルペアのなす角に関する情報を取り出した上で、上三角行列Rを作り、それを受けて、回転行列Qを作ることもできることがわかる
- この処理は、もちろん、計算的には重いのだが、やってみることにしよう
my.QR<-function(M){
k<-length(M[,1])
L<-sqrt(apply(M^2,2,sum))
M2<-t(t(M)/L)
T<-t(M2)%*%(M2)
X<-matrix(0,k,k)
X[1,1]<-1
for(i in 2:k){
tmp<-X[1:(i-1),1:(i-1)]
tmp2<-solve(tmp)
y<-tmp2%*%T[i,1:(i-1)]
y2<-sqrt(1-sum(y^2))
X[i,1:i]<-c(y,y2)
}
R<-t(X*L)
Q<-M%*%solve(R)
return(list(Q=Q,R=R,X=X,L=L,T=T,M=M))
}
k<-4
M<-matrix(rnorm(k^2),k,k)
my.QR.out<-my.QR(M)
Q<-my.QR.out$Q
R<-my.QR.out$R
L<-my.QR.out$L
X<-my.QR.out$X
T<-my.QR.out$T
Q%*%R-M
diag(L)%*%X-t(R)
Q%*%t(diag(L)%*%X)-M
X%*%t(X)-T
> Q%*%R-M
[,1] [,2] [,3] [,4]
[1,] 2.220446e-16 -1.110223e-16 0.000000e+00 5.551115e-17
[2,] 0.000000e+00 -5.551115e-17 0.000000e+00 1.387779e-16
[3,] 0.000000e+00 0.000000e+00 0.000000e+00 5.551115e-16
[4,] 0.000000e+00 0.000000e+00 1.110223e-16 4.440892e-16
> diag(L)%*%X-t(R)
[,1] [,2] [,3] [,4]
[1,] 0 0 0 0
[2,] 0 0 0 0
[3,] 0 0 0 0
[4,] 0 0 0 0
> Q%*%t(diag(L)%*%X)-M
[,1] [,2] [,3] [,4]
[1,] 2.220446e-16 -1.110223e-16 0.000000e+00 5.551115e-17
[2,] 0.000000e+00 -5.551115e-17 0.000000e+00 1.387779e-16
[3,] 0.000000e+00 0.000000e+00 0.000000e+00 5.551115e-16
[4,] 0.000000e+00 0.000000e+00 1.110223e-16 4.440892e-16
> X%*%t(X)-T
[,1] [,2] [,3] [,4]
[1,] 0 0.000000e+00 0.000000e+00 0.000000e+00
[2,] 0 1.110223e-16 0.000000e+00 1.110223e-16
[3,] 0 0.000000e+00 1.110223e-16 0.000000e+00
[4,] 0 1.110223e-16 0.000000e+00 0.000000e+00
>