ryamadaのコンピュータ・数学メモ RSSフィード

統計学・遺伝学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
数学用語集(TOMAC)

2018-06-30 整数係数多項式だから

[][]整数係数多項式 15:51

  • グラフのゼータ関数の逆数は、整数係数多項式であることが知られている
  • 無向グラフのエッジを2本の相互に逆向きのエッジに切り替え、そのうえで、その有向エッジをたどることで、Closed backtrackless tailess primitive paths (CBTPP)というのが出来たなら、それらすべてについて¥prod_{p ¥in ¥[P¥]} (1-u^{|p|})というのがゼータ関数の逆数
  • 今、|p|はサイクルの長さなのだが、グラフのエッジの長さが整数であるなら、ゼータ関数の逆数は整数係数多項式であることになる
  • グラフのゼータ関数の値は、行列式を使って出すことができ、ゼータ関数の逆数に出てくるuの最大次数はグラフのエッジ数の2倍であることもわかっているから、適当にゼータ関数の値を発生させれば、整数係数が連立方程式で算出できる
  • やってみる
library(igraph)
#el <- rbind(c(1,2),c(2,3),c(3,4),c(4,1),c(1,3))
el <- rbind(c(1,2),c(2,3),c(3,1),c(1,4))
g <- graph.edgelist(el,directed =FALSE)
plot(g)
my.bigraph <- function(g){
  el <- as_edgelist(g)
  el2 <- rbind(el,cbind(el[,2],el[,1]))
  ret <- graph.edgelist(el2)
  return(ret)
}
g.bi <- my.bigraph(g)
plot(g.bi,edge.label=1:8,edge.curved=TRUE)
my.WE <- function(g){
  el <- as_edgelist(g)
  n.e <- length(el[,1])
  edge.mat <- matrix(0,n.e,n.e)
	for(i in 1:n.e){
		st <- el[i,1]
		ed <- el[i,2]
		tmp <- which(el[,1]==ed & el[,2]!=st)
		edge.mat[i,tmp] <- 1
	}
  return(edge.mat)
}
WE <- my.WE(g.bi)
WE
WE2 <- WE[,c(2,7,3,4,5,6,1,8)]
WE2
u <- 3
WE2. <- WE2*u
WE2.
IWE2. <- (diag(rep(1,8))-WE*u)[,c(2,7,3,4,5,6,1,8)]
IWE2.
prod(diag(IWE2.))
(-u)^3
library(complexplus) # 複素行列計算用
my.Ihara.zeta.e <- function(g,u){
  g.bi <- my.bigraph(g) # 両方向グラフにする
  we <- my.WE(g.bi) # WE 行列を作る
  return(Det(diag(rep(1,length(we[,1]))) - we * u)) # 行列式を計算して返す
}
my.sign <- function(s){
  I <- diag(rep(1,length(s)))
  return(det(I[,s]))
}
u <- 0.3 + 1i * 0.2 # 適当な複素数
my.Ihara.zeta.e(g,u)
# 地道に計算するときは4つの置換
# (1,2,3,4,5,6,7,8) 置換なし。オリジナル
# (2,7,3,4,5,6,1,8) (1,2,7)
# (1,2,6,4,3,5,7,8) (3,6,5)
# (2,7,6,4,3,5,1,8) (1,2,7,3,6,5)
1 + (-u)^3 + (-u)^3 + (-u)^6
el <- rbind(c(1,2),c(2,3),c(3,4),c(4,1))
g <- graph.edgelist(el,directed =FALSE)
plot(g)
g.bi <- my.bigraph(g)
plot(g.bi,edge.label=1:10,edge.curved=TRUE)
WE <- my.WE(g.bi)
WE
u <- 3
WE2. <- WE2*u
WE2.
u <- 0.3 + 1i * 0.2 # 適当な複素数
my.Ihara.zeta.e(g,u)

1  +(-1)* 2*(-u)^4 + (-u)^8
my.Ihara.zeta.poly <- function(g){
  n.e <- length(E(g))*2
  us <- rnorm(n.e)
  U <- matrix(0,n.e,n.e)
  for(i in 1:length(us)){
    U[i,] <- (-us[i])^(1:n.e)
  }
  vs <- rep(0,n.e)
  for(i in 1:n.e){
    vs[i] <- my.Ihara.zeta.e(g,us[i])
  }
  A <- solve(U) %*% (vs-1)
  return(c(1,A))
}

my.Ihara.zeta.poly.calc <- function(g,u){
  A <- round(my.Ihara.zeta.poly(g))
  v <- sum(A * (-u)^(0:(length(A)-1)))
  return(list(v=v,A=A))
}
A <- my.Ihara.zeta.poly(g)
round(A,10)
my.Ihara.zeta.e(g,u)
my.Ihara.zeta.poly.calc(g,u)
トラックバック - http://d.hatena.ne.jp/ryamada/20180630

2018-06-23 グラフのゼータ関数(頂点ゼータ、エッジゼータ、パスゼータ)

[][]グラフのゼータ関数 その後 08:51

  • 資料
  • グラフの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])


トラックバック - http://d.hatena.ne.jp/ryamada/20180623

2018-06-19 ゼータ函数とは、のメモ

[]メモ 16:26

  • A directory of all known zeta functions
  • 数え上げる対象(自然数・素数、グラフの最短距離…)などがあったとき、その全体は「ぱぱっと」説明しにくい。そんなとき、良い感じの有理関数が定義できることがあって、それがゼータ関数
  • 関数は値を渡すと値を返す装置なので、k個の値を渡すとk個の値が返ってくる。無限個の値を渡すと無限個の値が返ってくるから、説明するために無限次元が必要な何かを体現させる装置としても有用
  • たとえば、分布関数の特性関数もそんな装置
トラックバック - http://d.hatena.ne.jp/ryamada/20180619

2018-06-17 完全グラフのゼータ関数

[][][][]完全グラフのゼータ関数 18:37

  • 完全グラフK_pはp個のノードを持つグラフで、すべてのノードペアにエッジがある(無向)グラフ
  • これのゼータ関数は
  • ¥zeta_{K_p}(u)^{-1} = (-1)^n ((p-2)u -1) (u+1)^n (u-1)^{n+1} ((p-2)u^2+u+1)^{p-1}
    • ただし、n=¥sum_{i=0}^p i -2p = ¥frac{p(p-3)}{2}
  • ちなみにp=0,1,2にも定義された式になっている。p=0,1,2の場合は項がキャンセルアウトして、結果としては、すべて¥zeta_{K_{p=0,1,2}}(u) = 1となる
my.ihara.n <- function(p){
	sum(0:p) - 2*p
}
for(i in 1:6){
  print(my.ihara.n(i))
}
[1] -1
[1] -1
[1] 0
[1] 2
[1] 5
[1] 9
  • 完全グラフ用の伊原ゼータ関数を作る関数
    • その上で、昨日の記事で書いた行列式による出力と比較する
my.ihara.zeta.Kp <- function(p){
	n <- my.ihara.n(p)
	ret <- function(u){
		1/((-1)^n * ((p-2)*u-1) * (u+1)^n * (u-1)^(n+1) * ((p-2)*u^2+u+1)^(p-1))
	}
	return(ret)
}

p <- 10
iharaKp <- my.ihara.zeta.Kp(p)

adj <- matrix(1,p,p)
diag(adj) <- 0
g <- graph.adjacency(adj,mode="undirected")

a <- runif(1) + runif(1)*1i
my.Ihara.zeta(g,a)
iharaKp(a)
トラックバック - http://d.hatena.ne.jp/ryamada/20180617

2018-06-16 グラフのゼータ関数 伊原とBartholdi

[][][][][][]グラフのゼータ関数 11:49

  • 資料
  • 伊原のゼータ関数を非正則グラフにまで拡張したのが橋本のゼータ関数で
    • ¥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のノードの次数を対角成分とする対角行列
    • ちなみに、グラフラプラシアンはL=D-Aで、これは、ゼータ関数で行列式を取っている部分のu=1に相当する(I- 1A +1^2(D-I) = I -A + D - I = D-A)。
  • Bartholdiのゼータ関数は、2つの複素数を取る関数
    • ¥zeta_G(u,t)^{-1}_{Barthroldi} = (1-(1-u)^2)t^2)^{m-n} det(I-tA + (1-u)(D-(1-u)I)t^2)
      • ここで¥zeta_G(u=0,t)^{-1}_{Barthroldi} = ¥zeta_G(t)^{-1}という関係が成り立つ
  • 以下では、複素行列の行列式の計算をするためにcomplexplusパッケージを使っている。伊原のゼータ関数(橋本版)は、Barthroldiの特殊版として書いてある
  • 検算例として、ノード数5の完全グラフK5を用いている。K5の伊原のゼータ関数の多項式表現はこちらにあるように(そしてsagemathを導入すれば引き出せるように)、(-1) ¥times (3¥times t - 1) ¥times (t + 1)^5 ¥times (t - 1)^6 ¥times (3 ¥times t^2 + t + 1)^4
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))
}
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)
}
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)))