逆べき乗法でさっさと最小固有値ベクトルを得る

• 離散外積代数・離散微分幾何で、大きな疎行列の計算をしたり、その「たとえばの固有値」の算出などをしている
• 全部の固有値と固有ベクトルを出すのは、重いけれど、一つの固有値と固有ベクトルを出すだけなら、「逆べき乗法」で出せばよいという
• 逆べき乗法のWikiはこちら
• 適当に大きな正方行列を作ってやってみよう
• 逆行列を計算してそれを繰り返して使うから、逆行列がある場合に使えます
• 確かに速い
• おそらく、疎行列の場合には、疎行列パッケージを使えばメモリのことも気にしなくてよいはず
• 実は繰り返しさえせず、たった一度でも、そこそこの近似解は得られるらしい（そのときの初期ベクトルは、(1,1,...)らしい？）

n <- 1000
X <- matrix(rnorm(n^2),n,n)
# 逆べき乗法
# 初期固有ベクトル候補は適当に決める
y <- rnorm(n)
# 繰り返し数（結構すぐに収束するが、ここは大きめにしておく）
k <- 100
ys <- matrix(0,k+1,n)
ys[1,] <- y
# 繰り返して使う逆行列はあらかじめ計算しておく
X.inv <- solve(X)
for(i in 1:k){
	ys[i+1,] <- X.inv %*% ys[i,]
	ys[i+1,] <- ys[i+1,] / sqrt(sum(ys[i+1,]^2))
}
# 収束の具合をプロット
matplot(ys,type="l")
# 固有値を全部出そう
eigen.out <- eigen(X)
# 固有値の絶対値が最小のものを取り出す
Mod(eigen.out[[1]])
this <- which(Mod(eigen.out[[1]])==min(Mod(eigen.out[[1]])))[1]
# 2法の結果を比べる
plot(eigen.out[[2]][,this],ys[k+1,])

• 関数化

# Aは実正方行列
my.inv.pow <- function(A,n.iter=3,b=rep(1,ncol(A)),log=FALSE){
	x <- b
	if(log){
		x.log <- matrix(0,n.iter+1,ncol(A))
		x.log[1,] <- x
	}
	#x <- x/sqrt(sum(x^2))
	A. <- solve(A)
	for(i in 1:n.iter){
		x <- A. %*% x
		x <- x/sqrt(sum(x^2))
		if(log){
			x.log[i+1,] <- x
		}
		
	}
	if(log){
		return(list(x=x,x.log=x.log))
	}else{
		return(list(x=x,x.log=matrix(0,0,ncol(A))))
	}
}