2021-03-07

Matrix Balancing、Sinkhornの定理

Matrix Balancing Sinkhornの定理二重確率行列置換行列 Birkhoff polytope

二重確率行列と言う行列がある
正方行列であって、すべての行和とすべての列和が１であり、かつ、行列の全成分が０以上(０より大)
二重確率行列を確率推行列とする確率推移では全要素値が均一状態に収束する
二重確率行列は、置換行列の線形和 $\sum_{i=1}^k p_i P_i$ と表せて、 $\sum_{i=1}^k p_i = 1; p_i >0$ となるという
もちろん、n!通りある、置換行列のすべてについて足し合わせてこの条件を満足することもできるが、実際には、n!よりもっと少ない置換行列の線形和で表せることが知られている
$\sum_{i=1}^k p_i P_i$ からわかるように、二重確率行列は置換行列が張る凸包の表面に相当する
二重確率行列の成分数は $n^2$ であるところ、置換行列の数はn!となっており、 $2^n << n!$ であるから、n!通りの置換行列を使わずとも、限定的な個数の置換行列の線形和で表せることになっている
さて。
任意の要素　正　なる正方行列　Aは $D A E$ と二つの正方行列でサンドイッチすることで、二重確率行列にすることができる。Sinkhornの定理

そして、D,Eは対角行列であって、対角成分は正だという
正成分の正方行列を二重確率行列にすることをmatrix balancingと言う
D,Eは単純な繰り返し計算により、近似的に求まることも知られている
そのアルゴリズムは複数、知られているが、以下の方法は、単純なアルゴリズムである
Aの成分 $a_{ij}$ は、i行の和とj列の和とを参考にして変化させたい値である
もし、i行の和だけで標準化するだけでよければ、i行の和で割れば、行和１を全行について達成できる
列の方だけを気にするなら、列和で割ればよい
今は行も列も整えたいので、 $a_{ij}$ 成分は、i行和とj列和の平方根の積で割ると、「良い感じ」
それを繰り返す

n <- 100
A <- matrix(runif(n^2),n,n)

As <- Ds <- Es <- list()
As[[1]] <- A
Ds[[1]] <- diag(rep(1,n))
Es[[1]] <- diag(rep(1,n))

# 繰り返し回数
n.iter <- 30

err <- rep(0,n.iter+1)
err[1] <- sum((apply(As[[1]],1,sum)-1)^2) + sum((apply(As[[1]],2,sum)-1)^2)

for(i in 1:n.iter){
	As[[i+1]] <- Ds[[i]] %*% As[[1]] %*% Es[[i]]
	tmp.r <- apply(As[[i+1]],1,sum) # 行和
	tmp.c <- apply(As[[i+1]],2,sum) # 列和
	tmp.d <- 1/sqrt(tmp.r) # 行和の平方根
	tmp.e <- 1/sqrt(tmp.c) # 列和の平方根
	Ds[[i+1]] <- diag(tmp.d) %*% Ds[[i]] # 行和に関して調整するのは対角行列を左からかけること
	Es[[i+1]] <- diag(tmp.e) %*% Es[[i]] # 列和に関して調整するのは対角行列を右からかけること
	# 行和、列和の１からのずれの二乗和を繰り返しごとに計算しておく
	err[i+1] <- sum((apply(As[[i+1]],1,sum)-1)^2) + sum((apply(As[[i+1]],2,sum)-1)^2)

}
# 結果の確認
apply(As[[n.iter+1]],1,sum)
apply(As[[n.iter+1]],2,sum)

plot(err)

# 関数化
my.sinkhorn.matbalance <- function(M,n.iter = 30){
	n <- length(M[,1])
	As <- Ds <- Es <- list()
	As[[1]] <- M
	Ds[[1]] <- diag(rep(1,n))
	Es[[1]] <- diag(rep(1,n))
	err <- rep(0,n.iter+1)
	err[1] <- sum((apply(As[[1]],1,sum)-1)^2) + sum((apply(As[[1]],2,sum)-1)^2)

	for(i in 1:n.iter){
		As[[i+1]] <- Ds[[i]] %*% As[[1]] %*% Es[[i]]
		tmp.r <- apply(As[[i+1]],1,sum) # 行和
		tmp.c <- apply(As[[i+1]],2,sum) # 列和
		tmp.d <- 1/sqrt(tmp.r) # 行和の平方根
		tmp.e <- 1/sqrt(tmp.c) # 列和の平方根
		Ds[[i+1]] <- diag(tmp.d) %*% Ds[[i]] # 行和に関して調整するのは対角行列を左からかけること
		Es[[i+1]] <- diag(tmp.e) %*% Es[[i]] # 列和に関して調整するのは対角行列を右からかけること
		# 行和、列和の１からのずれの二乗和を繰り返しごとに計算しておく
		err[i+1] <- sum((apply(As[[i+1]],1,sum)-1)^2) + sum((apply(As[[i+1]],2,sum)-1)^2)
	}
	return(list(double.st.mat=As[[n.iter+1]],D=Ds[[n.iter+1]],E=Es[[n.iter+1]],As=As,Ds=Ds,Es=Es,err=err))
}

M <- matrix(runif(n^2),n,n)

out <- my.sinkhorn.matbalance(M)
apply(out$double.st.mat,1,sum)
apply(out$double.st.mat,2,sum)

> apply(out$double.st.mat,1,sum)
  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [45] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [89] 1 1 1 1 1 1 1 1 1 1 1 1
> apply(out$double.st.mat,2,sum)
  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [45] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [89] 1 1 1 1 1 1 1 1 1 1 1 1