疎行列で逆冪乗法

• 逆冪乗法では、逆行列を出して、初期ベクトルに対して、それを何度か掛けて近似解を出す
• 逆行列の計算が重いので、逆行列を先に計算して、それを繰り返し使うというのが、普通のやり方
• ただし、疎行列の場合には、それをやると、逆行列が、疎行列ではなくなるので、逆行列を作って「持って」おくことができない
• その代り、逆行列を陽に作らずに、Ax=bのxを算出する演算を逆行列ライブラリの関数にやらせることで、可能にする
• 以下、my.inv.pow()は逆行列を出すパターン、my.inv.pow2()はsolve()関数に任せて逆行列を陽に作らないパターン

# 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))))
	}
}

my.inv.pow.2 <- 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 <- solve(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))))
	}
}
tmp.Es.v <- list()
for(i in 1:4){
	tmp.Es.v[[i]] <- sparseVector(rnorm(16),i=1:16,length=16)
}

N <- 20
X <- matrix(rnorm(N^2),N,N)
b <- rnorm(N)
tmp.outs <- my.inv.pow(X,n.iter=100,b=b,log=TRUE)
tmp.outs.2 <- my.inv.pow(X,n.iter=100,b=b,log=TRUE)