複素行列のdeterminant

• 行列の指数関数というのがある
• 簡単に言うと、固有値を指数の肩に乗せて、それを固有ベクトルの行列でサンドイッチする

exp.m <- function(A,n=1,eigen=FALSE){
  # 固有値分解
	eigen.out<-eigen(A)
	# P=V,P^{-1}=U
	V<-eigen.out[[2]]
	U<-solve(V)
	B<-diag(exp(eigen.out[[1]]*n))
	X <- V%*%B%*%U
  if(eigen){
    return(list(matrix = X, eigen.vs <- eigen.out[[1]]))
  }else{
    return(X)
  }
}

• じゃあ、行列の対数関数は？というと、これと同じで、以下のようにすればよいらしい

log.m <- function(A,n=1,eigen=FALSE){
  # 固有値分解
	eigen.out<-eigen(A)
	# P=V,P^{-1}=U
	V<-eigen.out[[2]]
	U<-solve(V)
	B<-diag(log((0*1i+eigen.out[[1]])*n)) # ここがlog()
	X <- V%*%B%*%U
	if(eigen){
		return(list(matrix = X, eigen.vs = eigen.out[[1]]))
	}else{
		return(X)
	}
}

• 確かめておこう

# 実行列
n <- 4
A <- matrix(rnorm(n^2),n,n)
exp.m(log.m(A)) - A

# 複素行列
my.random.complex <- function(n,r.rg=c(0,1),i.rg=c(0,1)){
  runif(n,min=r.rg[1],max=r.rg[2]) + 1i * runif(n,min=i.rg[1],max=i.rg[2])
}
A <- matrix(my.random.complex(n^2),n,n)
exp.m(log.m(A)) - A

• さて
• 指数行列のdeterminantはなのだという
• じゃあ、任意の（複素）行列 B についてなるAを持ち出してきてやれば、Bのdeterminantはとすればよいではないか？

my.trace <- function(M){
  sum(diag(M))
}
my.det.complex <- function(B){
	# 行列の対数関数をとって、そのトレースのexpを返してもよいし
	# 行列の固有値の対数をとって、その和のexpを返してもよい
	A <- log.m(B)
	#eigen.out<-eigen(B)
	#s <- sum(log(0*1i+eigen.out[[1]]))
	#return(exp(s))
	exp(my.trace(A))
	#exp(s)
}

• 確かめてみよう
  • 2x2の複素行列ならdeterminantの計算は簡単

my.det.2x2 <- function(M){
	M[1,1]*M[2,2] - M[1,2]*M[2,1]
}
n <- 2
A <- matrix(rnorm(n^2),n,n)
my.det.2x2(A)
my.det.complex(A)
A <- matrix(my.random.complex(n^2),n,n)
my.det.2x2(A)
my.det.complex(A)

• • 大きな正方行列になると直に計算するのは大変
  • たとえば、ユニタリ行列のdeterminantの絶対値は1だというので確かめてみる

n <- 4
n.mat <- 10
B.list <- list()
for(i in 1:n.mat){
  B <- matrix(my.random.complex(n^2),n,n)
  B. <- B + my.Herm.t(B)
  B.list[[i]] <- eigen(B.)[[2]]
}
lapply(B.list,function(M){Mod(my.det.complex(M))})

• • また、たとえば、を任意の複素行列で確かめる

n.mat <- 100
B.list <- list()
out.det <- out.tr <- rep(0,n.mat)
for(i in 1:n.mat){
  B <- matrix(my.random.complex(n^2),n,n)
  out.det[i] <- my.det.complex(exp.m(B))
  out.tr[i] <- exp(my.trace(B))
}
range(Mod(out.det-out.tr))