近似計算

• 参考(こちら)
• 正方行列を行と列の両方についてstochastic(足して１)になるように補正する処理
• こちら
• これを使うと数独も解けるという(こちら)

# sinkhorn

sinkhorn<-function(A,epsilon=0.00001){
	B<-A
	
	row_sum<-apply(A,1,sum)
	while(max(abs(t(row_sum)-1))>epsilon){
		B<-B/row_sum
		col_sum<-apply(B,2,sum)
		B<-t(t(B)/col_sum)
		row_sum <- apply(B,1,sum)
	}
	B
}
N<-3
A<-matrix(runif(N^2),N,N)

epsilon<-10^(-10)
shA<-sinkhorn(A,epsilon)
apply(shA,1,sum)
apply(shA,2,sum)

• Huber/LawファクターとMincファクター
  • Mincファクターは
  • それの漸近近似がHuber/Law
• 近似の関数をMATLABから書き直してみた

# sinkhorn

sinkhorn<-function(A,epsilon=0.00001){
	B<-A
	scaleX<-scaleY<-rep(1,length(A[,1]))
	row_sum<-apply(A,1,sum)
	while(max(abs(t(row_sum)-1))>epsilon){
		B<-B/row_sum
		scaleX<-scaleX/row_sum
		col_sum<-apply(B,2,sum)
		scaleY<-scaleY/col_sum
		B<-t(t(B)/col_sum)
		row_sum <- apply(B,1,sum)
	}
	list(mat=B,scaleX=scaleX,scaleY=scaleY)
}
N<-4
A<-matrix(runif(N^2),N,N)

epsilon<-10^(-10)
shA<-sinkhorn(A,epsilon)
apply(shA[[1]],1,sum)
apply(shA[[1]],2,sum)

# Huber/Law factor

HuberLaw<-function(x){
	x<-x
	x[which(x==0)]<-x[which(x==0)]+0.5
	hl<-x
	mores<-which(hl>1)
	unMores<-which(hl<=1)
	hl[mores]<-hl[mores]+0.5*log(hl[mores])+exp(1)-1
	hl[unMores]<-hl[unMores]*(exp(1)-1)+1
	hl/(exp(1))
}

Minc<-function(x){
	toosmall<-which(x<10^(-10))
	x2<-x
	x2[toosmall]<-10^(-10)
	exp(lgamma(x2+1)*(1/x2))
}

a<-seq(from=0.5,to=1.1,length=100)
#a<-0:1000
HL<-HuberLaw(a)
Mc<-Minc(a)
plot(a,HL/Mc)

# 
permanentEstimation<-function(A,iter=10000,epsilon=10^(-4),factor="hl"){
	n<-length(A[1,])
	shOut<-sinkhorn(A,epsilon)
	#print(shOut)
	B<-shOut[[1]]
	x<-shOut[[2]]
	y<-shOut[[3]]
	row_scale<-1/apply(B,1,max)
	C<-B*row_scale
	save_C<-C
	
	number_of_successes<-0
	
	for(i in 1:iter){
		column<-1
		C<-save_C
		row_sums<-apply(C,1,sum)
		while(column<=n){
			#print(column)
			if(factor=="hl"){
				hl<-prod(HuberLaw(row_sums))
			}else if(factor=="minc"){
				hl<-prod(Minc(row_sums))
			}
			if(factor=="hl"){
				tmphl2<-HuberLaw(row_sums-C[,column])
			}else if(factor=="minc"){
				tmphl2<-Minc(row_sums-C[,column])
			}
			hl2<-prod(tmphl2)
			#print(hl)
			#print(hl2)
			#print(tmphl2)
			row_probs<-C[,column]/tmphl2*hl2/hl
			#print(row_probs)
			# 重みづけサンプリング
			cumsumprob<-cumsum(row_probs)
			rrand<-runif(1)
			row_pick<-sum(cumsumprob<rrand)+1
			#print(row_pick)
			if(row_pick==(n+1)){
				column<-n+2
			}else{
				row_sums<-row_sums-C[,column]
				C[row_pick,]<-rep(0,n)
				column<-column+1
			}
			row_sums[row_pick]<-0
		}
		if(column==(n+1)){
			number_of_successes=number_of_successes+1
		}
	}
	C<-save_C
	row_sums<-apply(C,1,sum)
	
	if(factor=="hl"){
		hl_C<-prod(HuberLaw(row_sums))
	}else if(factor=="minc"){
		hl_C<-prod(Minc(row_sums))
	}
	per_estimate<-hl_C*number_of_successes/iter
	per_estimate<-per_estimate/prod(row_scale)/prod(x)/prod(y)
	
	per_estimate
}
PermanentBeta<-function(A){
	n<-length(A[,1])
	perms<-permutations(n,n)
	p<-0
	for(i in 1:length(perms[,1])){
		t<-1
		for(j in 1:n){
			t<-t*A[j,perms[i,j]]
		}
		p<-p+t
	}
	p
}
N<-100
A<-matrix(runif(N^2),N,N)

Niter<-1000

EstP<-permanentEstimation(A,iter=Niter)
EstPMinc<-permanentEstimation(A,iter=Niter,factor="minc")
#RyserP<-PermanentRyser(A)
#FullP<-PermanentBeta(A)

print(EstP)
print(EstPMinc)
#print(RyserP)
#print(FullP)