ベクトル積とクロネッカー積

• 行列のベクトル化 VectorizationのWiki記事はこちら
  • 行列を行ベクトルにすること。行列の列ベクトルをタンデムにつないで作る
• ３つの行列の積を、行列のクロネッカー積と行列をベクトル化したものとの積に切り替えることができる
• ３つの行列A,X,Bがあって、行列の積があったとする
• ここで、XとYとをそれぞれベクトル化してと書くことにする
• という関係がある
• ここではクロネッカー積
  • 行列方程式と言うらしい(こちら)
• なんのことやらわからないのでRでやってみる
• 

A <- matrix(runif(3*4),3,4)
X <- matrix(runif(4*5),4,5)
B <- matrix(runif(5*6),5,6)
A %*% X %*% B -> Y

• Vectorization

vecX <- c(X)
vecY <- c(Y)

• クロネッカー積

tB <- t(B)
tB %x% A
(tB %x% A) %*% vecX - vecY # 0ベクトル

• 以下は、上記と個人的にはつながりがある話だけれど、筆者本人以外には何のメリットもない話→自分用のメモ
  • 今、複数の正方行列のクロネッカー積をとりたい。実は、その個々の正方行列のある行は無視して、非正方行列のクロネッカー積もとりたい。そのとき、正方行列の積を取りつつ、最終的に、非正方行列のクロネッカー積を、正方行列の正方行列から行を抜き取る、という手順で作成したい、という要請がある。それをどうするか、という話
  • 諸事情あって、正方行列は行列のサイズを指定してある関数で作成されるようなものとし、無視したい行は必ず最後の１行とする…

# 正単体に関連する、ちょっとした関数
CategoryVector<-function (d = 3){
	df <- d - 1
	diagval <- 1:d
	diagval <- sqrt((d)/df) * sqrt((d - diagval)/(d - diagval + 1))
	others <- -diagval/(d - (1:d))
	m <- matrix(rep(others, d), nrow = d, byrow = TRUE)
	diag(m) <- diagval
	m[upper.tri(m)] <- 0
	as.matrix(m[, 1:df])
}
make.simplex <- function(k){
	cv <- CategoryVector(k)
	rbind(t(cv*sqrt(1-1/k)),rep(1/sqrt(k),k))
}
# 関数のサイズを並べたベクトルrをとって、左側にだんだんクロネッカー積を重ねていく
make.simplex.multi <- function(r){
	X <- make.simplex(r[1])
	k <- length(r)
	if(k > 1){
		for(i in 2:k){
			X <- make.simplex(r[i]) %x% X
		}
	}
	X
}
# 採用する行を1、省く行を0として返す
make.vector0 <- function(r){
	tmp <- list()
	for(i in 1:length(r)){
		tmp[[i]] <- c(rep(1,r[i]-1),0)
	}
	apply(expand.grid(tmp),1,prod)
}

# r はテーブルの形
# クロネッカー積、省く行を表すベクトル、省いた行列、行列の逆行列を返す
make.X.vector0 <- function(r){
	X <- make.simplex.multi(r)
	X.inv <- solve(X)
	vector0 <- make.vector0(r)
	non.zero <- which(vector0==1)
	return(list(X=X,X.sub=X[non.zero,],X.inv=X.inv,X.inv.sub=X.inv[,non.zero],vector0=vector0))
}

D <- matrix(c(1,2,-3,-1,-2,3),3,2)
r <- dim(D)
out <- make.X.vector0(r)
print(F <- out$X %*% c(D))
print(F2 <- out$X.sub %*% c(D))
out$X.inv %*% F
out$X.inv.sub %*% F2

> print(F <- out$X %*% c(D))
              [,1]
[1,]  1.732051e+00
[2,]  5.000000e+00
[3,]  0.000000e+00
[4,]  0.000000e+00
[5,] -2.220446e-16
[6,]  0.000000e+00
> print(F2 <- out$X.sub %*% c(D))
         [,1]
[1,] 1.732051
[2,] 5.000000
> out$X.inv %*% F
     [,1]
[1,]    1
[2,]    2
[3,]   -3
[4,]   -1
[5,]   -2
[6,]    3
> out$X.inv.sub %*% F2
     [,1]
[1,]    1
[2,]    2
[3,]   -3
[4,]   -1
[5,]   -2
[6,]    3