Matrix Balancing、Sinkhornの定理

  • 二重確率行列と言う行列がある
  • 正方行列であって、すべての行和とすべての列和が1であり、かつ、行列の全成分が0以上(0より大)
  • 二重確率行列を確率推行列とする確率推移では全要素値が均一状態に収束する
  • 二重確率行列は、置換行列の線形和 \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の定理

シンクホーンの定理 - Wikipedia

  • そして、D,Eは対角行列であって、対角成分は正だという
  • 正成分の正方行列を二重確率行列にすることをmatrix balancingと言う
  • D,Eは単純な繰り返し計算により、近似的に求まることも知られている
  • そのアルゴリズムは複数、知られているが、以下の方法は、単純なアルゴリズムである
  • Aの成分 a_{ij}は、i行の和とj列の和とを参考にして変化させたい値である
  • もし、i行の和だけで標準化するだけでよければ、i行の和で割れば、行和1を全行について達成できる
  • 列の方だけを気にするなら、列和で割ればよい
  • 今は行も列も整えたいので、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]] # 列和に関して調整するのは対角行列を右からかけること
	# 行和、列和の1からのずれの二乗和を繰り返しごとに計算しておく
	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]] # 列和に関して調整するのは対角行列を右からかけること
		# 行和、列和の1からのずれの二乗和を繰り返しごとに計算しておく
		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