グラフのゼータ関数 その後

  • 資料
  • グラフのvertex zeta 関数はノード数xノード数の行列を使って以下のように表せる
    • \zeta_G(u)^{-1} = (1-u^2)^{m-n} det(I-uA + u^2(D-I))
      • ただし、uは複素数、mとnはグラフGのエッジ数とノード数、det(M)はMの行列式で、Iはノード数xノード数の単位行列、AはGの隣接行列、DはGのノードの次数を対角成分とする対角行列
  • また、同じvertex zeta 関数はエッジ数xエッジ数の行列Wを使って
    • \zeta_G(u)^{-1} = det(I-W u)
      • ただし、uは複素数det(M)はMの行列式で、Wは(エッジ数x2)x(エッジ数x2)の行列である。グラフは無向グラフであり、そのエッジに2方向を考慮して2エッジを対象とする。今、あるエッジの終点があるエッジの始点であるとき、その2エッジのペアに対応するWのセルに1を立て、他は0とする。ただし、始点と終点を取り換えただけのエッジペアに対応するセルには0を立てる
    • さらに。エッジ長がすべて1のときは、上記でよいが、エッジ長が1とは限らないときには工夫が必要となる。エッジ長がすべて整数であるときには、その長さに応じてノードを追加し、すべてのエッジが1のグラフに変換すればよい
  • ちなみに、エッジの行列の行列式でさらっとvertex zeta 関数になるのは:
    • エッジ連結について、全部1が立っているとパスはそれを掛け合わせると、相変わらず1で、それが「存在する」という情報になるから。したがって、そこにひとしなみに複素数uを掛けるとu^kの項がエッジ数だけ冪乗されて現れる。これは、行列Wの要素をuとすることに同じ
    • さらに、エッジに長さ重みを入れたいときには、この原理を使って、行列要素をu^{{L(a)+L(b)}/2}とエッジごとに変えることで実現できる。ちなみにL(a),L(b)はエッジペア(a,b)のそれぞれのエッジの重み
  • さらに小さな行列の行列式計算に変形することもできる。それがpath zeta 関数
    • グラフに全域木を取り出し、全域木に含まれないエッジを取り出す。この全域木に含まれないエッジ数の2倍の二乗のサイズの行列の行列式計算にできる
    • 補全域木グラフの両方向を考えたエッジの集合を考える
    • それらのペアに対して、自身の逆方向にあたるエッジとのペアのときには、対応する行列成ー分を0にする
    • それ以外の場合には、ゼータ関数のuに重みを考慮した値を与える
    • その値は、補全域木グラフの2つのエッジが接続した木の2点間について、木の上のパス(最短)について、最短距離を取る。edge zeta 関数ではつながるエッジペアについて、2エッジの長さの平均を用いた。それと同様にエッジのつながりについて足し合わせて行きたいが、エッジペアをタンデムに足し合わせていくことは、結局そのグラフ距離を計算することのようなもの
    • こちらの9章にmathematica のコードgあるのでそれを読むのが早そう
  • 以下は、ノード数xノード数の行列でのゼータ関数、(エッジ数x2)x(エッジ数x2)の行列でのゼータ関数、エッジ長に合わせてノードを補う関数である
library(igraph)
library(complexplus)
my.Ihara.zeta.elem <- function(g){
	A <- as.matrix(as_adjacency_matrix(g))
	n.v <- length(V(g))
	n.e <- length(E(g))
	r <- n.e - n.v
	D <- diag(degree(g))
	H <- diag(degree(g)-1) 
	I <- diag(rep(1,n.v))
	return(list(r=r,A=A,H=H,I=I,D=D,n.v=n.v,n.e=n.e))
}
my.Ihara.zeta.ori <- function(g,u){
	elem <- my.Ihara.zeta.elem(g)
	1/((1-u^2)^elem$r * Det(elem$I-u*elem$A + u^2*elem$H))
}

my.Bartholdi.zeta <- function(g,u,t){
	elem <- my.Ihara.zeta.elem(g)
	tmp <- (1-(1-u)^2*t^2)^elem$r * Det(elem$I-t*elem$A+(1-u)*(elem$D-(1-u)*elem$I)*t^2)
	return(1/tmp)
}

my.Ihara.zeta <- function(g,u){
	my.Bartholdi.zeta(g,0,u)
}
# (エッジ長x2)^2 行列でのゼータ関数
my.Ihara.zeta.e <- function(g,u){
	el <- as_edgelist(g)
	el2 <- rbind(el,cbind(el[,2],el[,1]))
	n.e <- length(el[,1])
	edge.mat <- matrix(0,n.e*2,n.e*2)
	for(i in 1:(n.e*2)){
		st <- el2[i,1]
		ed <- el2[i,2]
		tmp <- which(el2[,1]==ed & el2[,2]!=st)
		edge.mat[i,tmp] <- 1
	}
	I <- diag(rep(1,n.e*2))
	tmpmat <- I - edge.mat * u
	return(1/Det(tmpmat))
}

# エッジに整数長さがあるときに、すべてのエッジ長が1となるように
# ノードを加えたグラフに変換する
# E(g)$weight <- sample(1:20,n.e)

my.graph.inflation.weight <- function(g){
	w <- E(g)$weight
	el <- as_edgelist(g)
	A <- as.matrix(as_adjacency_matrix(g))
	for(i in 1:length(w)){
		if(w[i] > 1){
			nv <- w[i]-1
			st <- el[i,1]
			ed <- el[i,2]
			A[st,ed] <- A[ed,st] <- 0
			n <- length(A[,1])
			newA <- matrix(0,n+nv,n+nv)
			newA[1:n,1:n] <- A
			newA[n+1,st] <- newA[st,n+1] <- 1
			newA[n+nv,ed] <- newA[ed,n+nv] <- 1
			if(nv > 1){
				for(j in 1:(nv-1)){
					newA[n+j,n+j+1] <- newA[n+j+1,n+j] <-1
				}
			}
			A <- newA
		}
	}
	newg <- graph.adjacency(A,mode="undirected")
	return(newg)
}

E(g)$weight <- rep(1,n.e)
newg <- my.graph.inflation.weight(g)


fu <- function(u){
	tmp <- (-1)*(3*u-1)*(u+1)^5*(u-1)^6*(3*u^2+u+1)^4
	return(1/tmp)
}

#g <- sample_gnp(10, 2/10)
n <- 5
adj <- matrix(1,n,n)
diag(adj) <- 0
g <- graph.adjacency(adj,mode="undirected")

x <- y <- seq(from=-0.8,to=0.8,length=100)

xy <- expand.grid(x,y)
xy. <- xy[,1] + 1i*xy[,2]

vs <- rep(0,length(xy.))
vs2 <- vs
for(i in 1:length(vs)){
	vs[i] <- my.Ihara.zeta(g,xy.[i])
	vs2[i] <- fu(xy.[i])
}
my.Ihara.zeta(g,vs[55])
fu(vs[55])
plot(xy,col=abs(Mod(vs)))


fu2 <- function(u){
	tmp <- (1-u^10)^5*(1-3*u^5)*(1-u^5)*(1+u^5+3*u^10)
	return(1/tmp)
}


fu(vs[515]^5)
fu2(vs[515])