回転してきれいにする

• こちらで正方行列の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)
	}
	# これからQR分解のRを作る
	R<-t(X*L)
	# QR分解のQを作る
	Q<-M%*%solve(R)
	# Q%*%R=M
	# diag(L)%*%X=t(R)
	# Q%*%diag(L)%*%X=M
	# X%*%t(X)=T
	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
>