• 岩波書店さんの『行列と変換群』をRでなぞりながら読むための副読本
  • こちらに現れる予定
  • Rmdファイル
---
title: "行列と変換群をRで"
author: "ryamada"
date: "Tuesday, December 09, 2014"
output: html_document
---

## 実二次元ベクトルを扱う

### 基本演算
ベクトルは値の組。
ベクトルの和、ベクトルの定数倍。

```{r}
x <- c(1,2)
y <- c(3,4)
x+y
a <- 3
a *x
```
### 内積
ベクトルに内積を定義する。
要素同士の積の和として定義してもよいし、2つのベクトルを行列とみなして、その片方の転置と積とみなしてもよい。
```{r}
my.ip.v <- function(x,y){
  #sum(x*y)
  t(x)%*%y
}
my.ip.v(x,y)
```

ノルムはベクトルとそれ自身の内積の平方根。
```{r}
my.norm <- function(x){
  sqrt(my.ip.v(x,x))
}
my.norm(x)
```

### 回転

回転行列を作る。
```{r}
my.2d.rot <- function(theta){
  matrix(c(cos(theta),sin(theta),-sin(theta),cos(theta)),2,2)
}
theta <- runif(1)*2*pi
rot <- my.2d.rot(theta)
x.rot <- rot %*% x
y.rot <- rot %*% y
```

### 内積は回転で不変
```{r}
my.ip.v(x,y)
my.ip.v(x.rot,y.rot)
```

### 外積

2次元ベクトルの外積はdetereminant
```{r}
det(cbind(x,y))
```
lこれも回転で不変。
```{r}
det(cbind(x.rot,y.rot))
```
## 回転行列

### 回転行列の性質
回転行列のdetermintは1。
回転行列の転置は逆行列。
```{r}
thetas <- runif(3)
rots <- list()
for(i in 1:length(thetas)){
  rots[[i]] <- my.2d.rot(thetas[i])
}
lapply(rots,det)
ts <- lapply(rots,t)
for(i in 1:length(thetas)){
  print(rots[[i]] %*% ts[[i]])
}
```

### 回転行列の積は回転行列〜回転群〜

群には、結合則、単位元、逆元がある。

結合則。
```{r}
(rots[[1]] %*% rots[[2]]) %*% rots[[3]]
rots[[1]] %*% (rots[[2]]) %*% rots[[3]]
```
単位元。
```{r}
rots.e <- diag(rep(1,2))
rots.e %*% rots[[1]]
rots[[1]] %*% rots.e
```
逆元は、逆行列。
```{r}
lapply(rots,solve)
# または、転置行列
lapply(rots,t)
```
## 2次元ベクトルと複素数

### 実部と虚部
2次元ベクトルの2つの要素を複素数の実部と虚部とみることにする。

```{r}
zx <- x[1] + 1i * x[2]
zy <- y[1] + 1i * y[2]
```

和と実定数倍はベクトルの場合と同じ。
```{r}
x + y
zx + zy
x * y[1]
zx * Re(zy)
```

### 内積と共役複素数
ベクトルの内積はちょっと違う。
ベクトルの内積を計算するには、共役複素数を持ち出す必要がある。

```{r}
zx. <- Conj(zx)
```
ベクトルの内積は
```{r}
my.complex.ip <- function(zx,zy){
  1/2*(zx * Conj(zy) + Conj(zx) * zy)
}
my.complex.ip(zx,zy)
my.ip.v(x,y)
```

### 回転とは複素数をかけること

回転行列は
```{r}
thetas <- runif(3)
rots <- list()
for(i in 1:length(thetas)){
  rots[[i]] <- my.2d.rot(thetas[i])
}
```
これに対応する複素数がある。
```{r}
rots.z <- list()
for(i in 1:length(thetas)){
  rots.z[[i]] <- exp(1i*thetas[i])
}
for(i in 1:length(thetas)){
  print(rots[[i]] %*% x)
  print(rots.z[[i]] * zx)
}
```

### ベクトルとしての複素数、行列としての複素数

2次元ベクトルの2要素を複素数の実部と虚部に対応付けるのが、「実2次元ベクトルの複素数化」。

他方、回転行列の作用を、複素数の積に対応付けるときは$e^{i\theta}$$\theta$回転を表す2x2行列に対応付けるわけなので、$1+i\times 0$と、純虚数$i$に対応するのは、
```{r}
matrix(c(1,0,0,1),2,2)
matrix(c(0,-1,1,0),2,2)
```

### 2次元回転はSO(2)でありU(1)でもある

2次元の回転を2x2行列で表せば、determinantが1である直交行列であった。それをSO(2) : determinantが1であるという特殊性(special)を持つ直交(Orthogonal)な次元2の群、と言う。
他方、複素数を用いると、$e^{i\theta}$と表せて、これを複素数を要素とする1x1ユニタリ行列とみることで、この群をU(1)と書く。ちなみに、ユニタリ行列とは『共役転置行列と逆行列が等しい行列』のこと。

言い換えると、2次元平面にある単位円に、2x2実行列と、1x1複素行列との二つの見方がある、ということ。

これを、3次元空間の単位球面に一般化していきましょう、というのが、この文書の狙い。

## 行列から群に進む前に、行列の基本性質の確認

### トレース

対角成分は特別だ、という性質

```{r}
# ランダウ複素数発生関数
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])
}
my.trace <- function(M){
  sum(diag(M))
}
n <- 4
A <- matrix(my.random.complex(n^2),n,n)
B <- matrix(my.random.complex(n^2),n,n)
my.trace(A+B) - (my.trace(A) + my.trace(B))
a <- runif(1)
my.trace(a*A) - a*(my.trace(A))
my.trace(A%*%B) - my.trace(B%*%A)
M <- matrix(my.random.complex(n^2),n,n)
my.trace(M%*% A %*% solve(M))-my.trace(A)
```

### 転置
```{r}
t(A+B)-(t(A)+t(B))
t(a*A) - a*t(A)
t(A%*%B) - t(B)%*%t(A)
my.trace(t(A)) - my.trace(A)
t(t(A)) - A
```
### 共役転置

ただ、転置するだけでなく、複素共役にした上で転置する。
```{r}
my.Herm.t <- function(M){
  t(Conj(M))
}
A
my.Herm.t(A)
my.Herm.t(A+B) - (my.Herm.t(A)+my.Herm.t(B))
my.Herm.t(a*A) - Conj(a)*my.Herm.t(A)
my.Herm.t(A%*%B) - my.Herm.t(B)%*%my.Herm.t(A)
my.trace(my.Herm.t(A))-Conj(my.trace(A))
my.Herm.t(my.Herm.t(A))-A
```
### 交換子積

行列の積は、順番を入れ替えると変わるので、その差を表すものとして交換子積がある。

```{r}
my.commutator <- function(A,B){
  A%*%B - B %*%A
}
my.commutator(A,solve(A))
my.commutator(A,diag(rep(1,length(A[1,]))))

my.commutator(A,B)-my.commutator(B,A)
my.commutator(A+B,M)-(my.commutator(A,M)+my.commutator(B,M))
b <- my.random.complex(1)
my.commutator(b*A,B)-b*my.commutator(A,B)
t(my.commutator(A,B)) - (-my.commutator(t(A),t(B)))
my.Herm.t(my.commutator(A,B))-(-my.commutator(my.Herm.t(A),my.Herm.t(B)))
```
### 特殊な行列クラスと群

#### 対称行列と反対称行列
対称行列は、行列とその転置の和で表せるし、反対称行列は、行列とその転置の差。
```{r}
A. <- A + t(A)
A.
t(A.)-A.
A.. <- A - t(A)
A..
t(A..) + A..
```
対称行列と反対称行列とは、交換子積に特徴を持つ。
```{r}
A. <- A + t(A)
B. <- B + t(B)
A.B. <- my.commutator(A.,B.)
t(A.B.) + A.B.
A.. <- A - t(A)
B.. <- B - t(B)
A..B.. <- my.commutator(A..,B..)
t(A..B..) + A..B..
```

エルミート行列は、実成分が対称行列、虚成分が反対称行列。
```{r}
A. <- Re(A + t(A)) + 1i * Im(A-t(A))
A.
my.Herm.t(A.)-A.

B. <- Re(A+t(A)) + 1i*Im(B-t(B))
A.+B.
my.Herm.t(A.+B.)-(A.+B.)
A.%*%B.+B.%*%A.
my.Herm.t(A.%*%B.+B.%*%A.)-(A.%*%B.+B.%*%A.)
1i * my.commutator(A.,B.)
my.Herm.t(1i*my.commutator(A.,B.))-1i*my.commutator(A.,B.)
```

反エルミート行列は、実部が反対称行列、虚部が対称行列。
```{r}
A.. <- Re(A-t(A)) + 1i*Im(A+t(A))
my.Herm.t(A..)+A..
```

### 直交行列・ユニタリ行列と群

直交行列のdeterminantは$\pm1$。

直交行列は直交群をなす。

determinantが1の場合は直交行列の部分集合で回転。これは回転群。

ユニタリ行列はユニタリ群をなす。
ユニタリ行列のdeterminantの絶対値が1。
これは群をなす。これはユニタリ群。

ユニタリ行列のdeterminantを1に限るとこれも群をなし、これは特殊ユニタリ群。

直交行列はユニタリ行列の実部版。

直交行列とユニタリ行列を作っておく。

対称行列の固有ベクトルが作る行列は直交行列である。
エルミート行列の固有ベクトルが作る行列はユニタリ行列である。

```{r}
n <- 4
A <- matrix(rnorm(n^2),n,n)
A. <- A+t(A)
t(A.)-A
Ax <- eigen(A.)[[2]]
Ax %*% t(Ax)
det(Ax)

B <- matrix(my.random.complex(n^2),n,n)
B. <- B + my.Herm.t(B)
my.Herm.t(B.) - B.
Bx <- eigen(B.)[[2]]
Bx %*% my.Herm.t(Bx)
```
複素行列であるユニタリ行列のdeterminantの絶対値が1であることの確かめは、いくつかの補助関数を作らないといけないけれど、先でも使うので、ここで作っておく。

```{r}
# 行列の対数関数
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))
	X <- V%*%B%*%U
	if(eigen){
		return(list(matrix = X, eigen.vs = eigen.out[[1]]))
	}else{
		return(X)
	}
}
# 行列の指数関数(後で章を立ててやる)
# 行列の指数関数と対数関数は同じような関数
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)
  }
	
}
# 複素行列のdeterminant関数
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)
}
my.det.complex(Bx)
Mod(my.det.complex(Bx))
```

ちなみに、determinantが1であるユニタリ行列を以下のように作ると、$\begin{pmatrix}\alpha & \beta \\ -\bar{\beta} & \bar{\alpha}\end{pmatrix}$となっていることもわかる。
```{r}
Bx/sqrt(my.det.complex(Bx))
```

ユニタリ行列の積がユニタリ行列であることも確かめておく。

```{r}
n <- 4
B.list <- list()
for(i in 1:3){
  B <- matrix(my.random.complex(n^2),n,n)
  B. <- B + my.Herm.t(B)
  B.list[[i]] <- eigen(B.)[[2]]
}
my.Herm.t(B.list[[1]]%*%B.list[[2]]) %*% (B.list[[1]]%*%B.list[[2]])
```


## 次元を上げる〜実ベクトルの高次元化・ベクトルの複素数化〜

### 高次元ベクトルの内積・複素ベクトルの内積
```{r}
my.ip.v <- function(x,y){
  #sum(x*y)
  t(x)%*%y
}
my.ip.v(x,y)

my.ip.complex.v <- function(zx,zy){
  #t(Conj(zx)) %*% zy
  my.Herm.t(zx) %*% zy
}
my.ip.complex.v(x,y)

zx <- my.random.complex(2)
zy <- my.random.complex(2)
zx
my.ip.complex.v(zx,zy)
```
複素ベクトル内積の特徴。
```{r}
my.ip.complex.v(zx,zy) - Conj(my.ip.complex.v(zy,zx)) # 順序を入れ替えて共役をとる
zz <- my.random.complex(2)
my.ip.complex.v(zx+zy,zz) - (my.ip.complex.v(zx,zz)+my.ip.complex.v(zy,zz))
b <- my.random.complex(1)
my.ip.complex.v(b*zx,zy)-Conj(b)*my.ip.complex.v(zx,zy)
my.ip.complex.v(zx,b*zy)-b*my.ip.complex.v(zx,zy)
my.ip.complex.v(zx,zx) - sum(Mod(zx)^2)
```

### 内積を不変に保つ変換

実ベクトル空間の場合、内積は直交変換で不変。
複素ベクトル空間の場合、内積は、ユニタリ変換で不変。。

```{r}
n <- 4
A <- matrix(rnorm(n^2),n,n)
A. <- A+t(A)
t(A.)-A
Ax <- eigen(A.)[[2]] # 直交行列

x <- rnorm(n)
y <- rnorm(n)
my.ip.v(Ax%*%x,Ax%*%y)-my.ip.v(x,y)

B <- matrix(my.random.complex(n^2),n,n)
B. <- B + my.Herm.t(B)
my.Herm.t(B.) - B.
Bx <- eigen(B.)[[2]] # ユニタリ行列
zx <- my.random.complex(n)
zy <- my.random.complex(n)
my.ip.complex.v(Bx%*%zx,Bx%*%zy)-my.ip.complex.v(zx,zy)
```

### 行列にして次元を上げる

行列は数を二次元配置していて、そこには演算があるから、行列も次元を上げる手段とみることができる。

今、内積を不変に保つことに着目しているから、行列という形で次元を上げるなら、行列の内積の定義が必要。

行列の内積。
```{r}
n <- 2
A <- matrix(my.random.complex(n^2),n,n)
B <- matrix(my.random.complex(n^2),n,n)
# 二通りの書き方がある。
sum(Conj(A) * B)
my.trace(my.Herm.t(A)%*%B)
# 後者で関数化しておく
my.ip.mat <- function(A,B){
  my.trace(my.Herm.t(A)%*%B)
}
```
内積の性質の確認。
```{r}
my.ip.mat(A,B) - Conj(my.ip.mat(B,A))
M <- matrix(my.random.complex(n^2),n,n)
my.ip.mat(A+B,M) - (my.ip.mat(A,M)+my.ip.mat(B,M))
b <- my.random.complex(1)
my.ip.mat(b*A,B) - Conj(b) * my.ip.mat(A,B)
my.ip.mat(A,b*B) - b * my.ip.mat(A,B)
my.ip.mat(A,A) - sum((Mod(A)^2))
```

行列の基底は、1つの要素が1であとは0であるような行列。
任意の行列は、この基底の線形和で表される。
基底の係数は、その基底と表したい行列の内積。
この関係も普通のベクトルと同じになっている(それは内積の定義のせい)。

### 行列の内積を不変に保つ線形変換

2つのユニタリ行列でサンドイッチにする変換を考える。
それによって、行列の内積が保たれることを示す。
```{r}
n <- 4
A1 <- matrix(my.random.complex(n^2),n,n)
A2 <- matrix(my.random.complex(n^2),n,n)
B.list <- list()
for(i in 1:2){
  B <- matrix(my.random.complex(n^2),n,n)
  B. <- B + my.Herm.t(B)
  B.list[[i]] <- eigen(B.)[[2]]
}
U <- B.list[[1]]
V <- B.list[[2]]

my.ip.mat(U%*%A1%*%V,U%*%A2%*%V) - my.ip.mat(A1,A2)
```

特殊な場合として、あるユニタリ行列とその逆行列(ユニタリ行列の逆行列はエルミート共役だが)とでサンドイッチする場合を考えることもある。
2x2の場合は、後述する3次元回転と関係してくる。
この場合は、内積だけでなく、determinantとトレースが保存される。
```{r}
my.trace(U%*%A1%*%my.Herm.t(U)) - my.trace(A1)
```

determinant とトレースが保存されるのは、Uがユニタリであるだけではなく、正則行列であれば、保存される。このような変換を掃除変換と呼ぶ。
```

### 群と行列
群の元を行列とし、
群の積を行列の積とし、
群の単位元を単位行列とし、
群の逆元を逆行列とする。

こうすることで、群に行列による表現ができる。

行列はベクトル空間における線形変換と見ることができるので、群が、ベクトル空間と結びついたことになる。

ベクトル空間では、線形変換が起きているから、「変換群」と呼ぶ。

一つの群には、複数の行列表現を対応付けることができる。その1例が、2次元の回転群に、2x2実行列の回転行列としての表現と、1x1虚数行列としての表現があることである。

## 行列の指数関数

行列の作用を連続的に表そうとするとき、行列の指数関数を使う。

対角成分をサンドイッチした形に分解し、対角成分についてだけ、パラメタに応じて変化を入れる。
その対角成分のパラメタ変化が指数関数になっている。

回転行列やエルミート行列の場合は、「ぐるぐる回る」ので、指数関数もぐるぐる回ったものになる。

```{r}
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)
  }
	
}
n <- 4
A <- matrix(my.random.complex(n^2),n,n)
exp.A <- exp.m(A)
exp.m(t(A)) - t(exp.m(A))
exp.m(my.Herm.t(A)) - my.Herm.t(exp.m(A))
P <- matrix(my.random.complex(n^2),n,n)
exp.m(P %*% A%*%solve(P)) - P%*%exp.m(A)%*%solve(P)
```
行列の指数関数の逆行列は、行列を負にしたものの指数関数。
```{r}
solve(exp.m(A)) -exp.m(-A)
```
行列の指数関数は行列であるが、2つの行列の指数関数の積を、2つの行列の線形和の指数関数で表すことは、通常、難しい。
それは、行列の指数関数は、行列を繰り返し掛け合わせたものになっているわけだが、行列の積の順序を入れ替えられないからである。

逆にいうと、順序を入れ替えてもよいような行列同士(交換子積が0)の場合には、行列の指数関数の積は、行列の和の指数関数になる。

## 行列の指数関数と微分と行列式

$\frac{de^{tA}}{dt} = A e^{tA}$なる関係がある。
これは、一階の微分方程式。

また、行列の指数関数のdeterminantとトレースとには次の関係がある。

```{r}
A <- matrix(rnorm(n^2),n,n)
det(Re(exp.m(A)))-exp(my.trace(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)
my.det.complex(exp.m(A))-exp(my.trace(A))
```

### 特徴的な行列とその指数関数

対称行列・反対称行列・エルミート行列・反エルミート行列の指数関数の特徴を示す。

```{r}
n <- 4
A <- matrix(my.random.complex(n^2),n,n)
B <- matrix(my.random.complex(n^2),n,n)
# 対称
A1 <- A + t(A)
# 対称行列の指数関数は対称行列
t(exp.m(A1)) - exp.m(A1) 
# 反対称
A2 <- A - t(A)
# 反対称行列の指数関数は直交行列
t(exp.m(A2)) - exp.m(-A2)
t(exp.m(A2)) - solve(exp.m(A2)) 
# エルミート
A3 <- Re(A+t(A)) + 1i*Im(A-t(A))
# エルミート行列の指数関数はエルミート行列
my.Herm.t(exp.m(A3)) - exp.m(A3)
# 反エルミート
A4 <- Re(A-t(A)) + 1i*Im(A+t(A))
my.Herm.t(A4)+A4
#A4 <- Re(A-t(A)) + 1i*Im(A+t(A))
# 反エルミート行列の指数関数逆行列はエルミート転置
my.Herm.t(exp.m(A4)) - solve(exp.m(A4))
# それはユニタリ
Mod(my.det.complex(exp.m(A4)))
```

### 変換群と行列との関係

> 直交群(O(n))の群の元はnxn直交行列。それは反対称行列の指数関数
> 特殊直交群(SO(n))の群の元は行列式1のnxn直交行列。それは反対称行列の指数関数

> ユニタリ群(U(n))の群の元はnxnユニタリ行列。それは反エルミート行列の指数関数

> 特殊ユニタリ群(SU(n))の群の元は行列式1のnxnユニタリ行列。それはトレース0の反エルミート行列の指数関数

## 3次元の回転

3次元回転は非可換。

オイラー角による3次元回転行列。
```{r}
# n次元、aが指定する平面に関する角度thetaの回転行列
my.rot.nd.12 <- function(theta,a=c(1,2),n=3){
  ret <- diag(rep(1,n))
  ret[a,a] <- matrix(c(cos(theta),sin(theta),-sin(theta),cos(theta)),2,2)
  ret
}
# 第3軸についてphi回してから、第1軸についてtheta回し、さらに第3軸についてpsi回す
my.rot.3d <- function(phi,theta,psi){
  rot1 <- my.rot.nd.12(phi,a=c(1,2))
  rot2 <- my.rot.nd.12(theta,a=c(2,3))
  rot3 <- my.rot.nd.12(psi, a=c(1,2))
  rot3 %*% rot2 %*% rot1
}
my.rot.3d(pi/3,pi/4,-pi/3)
```

### 行列の指数関数を使って回転行列
$\cos(\theta),\sin(\theta)$を用いて回転行列を作ったが、$\begin{pmatrix}0 & 1 \\ -1 & 0\end{pmatrix}$という行列の指数関数を使って作ることもできる。
この行列が反対称行列であることにも注意。

```{r}
my.rot.nd.12.e <- function(theta,a=c(1,2),n=3){
  ret <- diag(rep(1,n))
  A <- matrix(c(0,-1,1,0),2,2)
  ret[a,a] <- exp.m(A,theta)
  ret
}
theta <- pi/4
my.rot.nd.12.e(theta)
my.rot.nd.12(theta)
my.rot.3d.e <- function(phi,theta,psi){
  rot1 <- my.rot.nd.12.e(phi,a=c(1,2))
  rot2 <- my.rot.nd.12.e(theta,a=c(2,3))
  rot3 <- my.rot.nd.12.e(psi, a=c(1,2))
  rot3 %*% rot2 %*% rot1
}
my.rot.3d.e(pi/3,pi/4,-pi/3)
```

### 回転生成子の交換関係

$\begin{pmatrix}0 & 1 \\ -1 & 0\end{pmatrix}$のような行列を回転の生成子と言う。

これらには、交換関係がある。
交換子を問題にするのは、「交換子=0」ならば、可換であって、これが0出ないことが非可換を特徴付けているから。


```{r}
Js <- list()
A <- matrix(c(0,-1,1,0),2,2)
Js[[1]] <- matrix(0,3,3)
Js[[1]][c(2,3),c(2,3)] <- A
Js[[2]] <- matrix(0,3,3)
Js[[2]][c(3,1),c(3,1)] <- A
Js[[3]] <- matrix(0,3,3)
Js[[3]][c(1,2),c(1,2)] <- A

my.commutator(Js[[1]],Js[[2]]) + Js[[3]]
my.commutator(Js[[2]],Js[[3]]) + Js[[1]]
my.commutator(Js[[3]],Js[[1]]) + Js[[2]]
```

### 交換関係と置換
1,2,...,nの順序を入れ替えたものは置換。
ペアを入れ替えることを偶数回だけ繰り返すのは偶置換。
奇数回の場合は奇置換。

偶置換に1を、奇置換に0を与える「置換符合」を考える。

さらに進んで、n個の並びが置換ではなくて、一部に重複がある場合に、0となるような、「置換重複符合」を与えることにする。

```{r}
my.perm.sign <- function(s){
  n <- length(s)
  M <- diag(rep(1,n))
  det(M[s,])
}
my.perm.sign(c(1,2,3))
my.perm.sign(c(2,3,1))
my.perm.sign(c(3,1,2))
my.perm.sign(c(2,1,3))
my.perm.sign(c(1,3,2))
my.perm.sign(c(3,2,1))
my.perm.sign(c(1,1,2))
my.perm.sign(c(1,2,1))
my.perm.sign(c(2,1,1))
```
この生成子の交換子積が別の生成子に一致することが何やら重要で、この交換関係が「リー代数」だ、ということだが、ここではちょっとわからないので、後の章でこの件は解決したい。

なお、生成子は反対称行列であるが、この生成子の1,-1を虚数記号に換えると、トレースが0のエルミート行列になることにも注意。

## 3次元回転とユニタリ行列

3x3直交行列で表現できない下部構造を知るために、determinantが1の2x2ユニタリ行列と3次元回転の回転を確認する。

ユニタリ行列は、determinant の絶対値が1であるような行列で、エルミートが逆行列になっているものである。
2x2の場合には、$e^{i\phi/2} \begin{pmatrix}\alpha & \beta \\ -\bar{\beta} & \bar{\alpha} \end{pmatrix}$という形をしている。ここで$|\alpha|^2+|\beta|^2=1$である。

これを利用して、determinant が1のユニタリ行列、一般的な2x2ユニタリ行列を作る。

$\alpha = x_4 + i \times x_3, \beta=x_2+i \times x_1$とすると、4つの実数$x_i;i=1,2,3,4$から成り、$x_1^2+...+x_4^2=1$である。

これを
$$ i x_1\begin{pmatrix}0 & 1 \\ 1 & 0 \end{pmatrix} $$

$$ + i x_2\begin{pmatrix}0 & -i \\ i & 0 \end{pmatrix}$$

$$ + i x_3 \begin{pmatrix}1 & 0 \\ 0 & -1 \end{pmatrix}$$

$$ + x_4 \begin{pmatrix}1 & 0 \\ 0 & 1 \end{pmatrix}$$
と表せば、4つの行列が基礎的要素であることがわかる。
この4つの行列から単位行列を除いた3つは、パウリ行列と呼ばれる。
パウリ行列はトレースが0のエルミート行列。


```{r}
paulis <- list()
paulis[[1]] <- matrix(c(0,1,1,0),2,2)
paulis[[2]] <- matrix(c(0,1i,-1i,0),2,2)
paulis[[3]] <- matrix(c(1,0,0,-1),2,2)
paulis[[4]] <- matrix(c(1,0,0,1),2,2)
paulis

d <- exp(1i*runif(1)*2*pi)
library(MCMCpack)
x <- sqrt(rdirichlet(1,rep(1,4)))
U.s <- 1i * x[1] * paulis[[1]] + 1i * x[2] * paulis[[2]] + 1i * x[3] * paulis[[3]] + x[4] * paulis[[4]]
U <- d*U.s
my.det.complex(U.s)
my.det.complex(U)
d^2
apply(Mod(U.s)^2,1,sum)
apply(Mod(U)^2,1,sum)
```

パウリ行列の性質。
```{r}
paulis[[1]]^2
paulis[[2]]^2
paulis[[3]]^2
paulis[[1]] %*% paulis[[2]] -1i * paulis[[3]]
paulis[[1]] %*% paulis[[2]] %*% paulis[[3]]
lapply(paulis,my.det.complex)
paulis[[1]] %*% paulis[[2]] %*% paulis[[1]] + paulis[[2]]
paulis[[1]]^3 - paulis[[1]]
```

パウリ行列の(半分の)交換子積は、3x3行列の回転生成子と同じ性質がある。

```{r}
my.commutator(paulis[[1]]/2,paulis[[2]]/2) - 1i * my.perm.sign(c(1,2,3)) * paulis[[3]]/2
my.commutator(paulis[[2]]/2,paulis[[3]]/2) - 1i * my.perm.sign(c(2,3,1)) * paulis[[1]]/2
my.commutator(paulis[[3]]/2,paulis[[1]]/2) - 1i * my.perm.sign(c(3,1,2)) * paulis[[2]]/2
```

パウリ行列の$1/\sqrt{2}$は正規直交基底。
```{r}
ips <- matrix(0,3,3)
for(i in 1:3){
  for(j in 1:3){
    ips[i,j] <- my.ip.mat(paulis[[i]]/sqrt(2),paulis[[j]]/sqrt(2))
  }
}
ips
```

これは、3つのパウリ行列は、何かの3次元実ベクトル空間を張っていることを意味する。
それは、トレースが0のエルミート行列全体である。
言い換えると、トレースが0のエルミート行列全体は3次元ベクトル空間全体である。

### パウリ行列の線形和と球

パウリ行列の線形和は$$\begin{pmatrix}x_3 & x_1 - ix_2 \\ x_1+ix_2 & -x_3 \end{pmatrix}$$
と表される。

この行列の行列式は$$-(x_3^2+x_1^2+x_2^2)$$である。

今、パウリ行列はトレース0のエルミート行列であるが、$(x_1,x_2,x_3,x_4=0)$という(行列表記されたベクトル)について、内積が不変な変換では、行列式が不変ということになる。

determinantが1である2x2ユニタリ行列はトレース1のエルミート行列に対して、行列式不変な変換。

```{r}
# トレース0のエルミート行列
n <- 2
A <- matrix(my.random.complex(n^2),n,n)
A. <- Re(A + t(A)) + 1i * Im(A-t(A))
A.[n,n] <- A.[n,n]- my.trace(A.)
my.trace(A.)

# determinant 1 のユニタリ行列
B <- matrix(my.random.complex(n^2),n,n)
B. <- B + my.Herm.t(B)
my.Herm.t(B.) - B.
Bx <- eigen(B.)[[2]]
Bx <- Bx/sqrt(my.det.complex(Bx))
# ユニタリ変換Bs A. Bx'
my.det.complex(A.) - my.det.complex(Bx%*%A.%*%my.Herm.t(Bx))
my.Herm.t(Bx%*%A.%*%my.Herm.t(Bx)) - Bx%*%A.%*%my.Herm.t(Bx)
my.trace(Bx%*%A.%*%my.Herm.t(Bx))
```

したがって、パウリ行列の線形和は3次元空間だが、それを半径ごとに球面に分割することができて、それは、行列の世界では、determinantの値で分割することであることがわかる。

2x2ユニタリ変換は、内積不変な複素回転行列だから、パウリ行列をユニタリ変換で移り変われる者同士で分類すれば、それが半径による空間の球分類となる。

3x3回転行列をかけることと、2x2ユニタリ変換をすること(ユニタリ行列とそのエルミートで挟むこと)が同じで、いくつかの回転行列の積は、いくつかのユニタリ行列の積によるユニタリ変換に対応づく。

実際、このユニタリ変換を担うユニタリ行列は、パウリ行列の指数関数になっている。

### 3次元回転の2つの表現の違い

3次元回転を3x3直交行列で表すか、2次元ユニタリ行列で表すかの違い。

SO(3)では、ある角$\theta$の回転に関して、$\theta$そのものがパラメタとして用いられるのに対して、U(2)の方では$\theta/2$が用いられる。

これは$2\pi$の回転が「作用として何もしないのと同じ」とみなすか、$2\pi$ではなくて$4\pi$の回転(2回転)して初めて「作用として何もしないのと同じ」とみなすか、の違いとなる。

この「$2\pi$の回転は、座標としては同じだが、何かが違う」という情報が「スピン」に関する情報であるとして、スピノールという概念が登場する(量子力学)## ベクトル・テンソル・スピノール

### テンソルの計算
n次元空間での回転を考えるときに気になるのは、スカラー(n^0個の要素からなる0階のテンソル)、n次元ベクトル(n^1個の要素からなる1階のテンソル)、n^2個の要素からなる2階のテンソル、n^3個の要素からなる3階のテンソル…、n^k個の要素からなるk階のテンソル。

k階のテンソルは、k個の添え字を持ち、各添え字は1,2,...,nのいずれかである。

テンソルは、nxn行列(変換行列)によって、次の規則で変換する。

$$
t_{i_1,i_2,...,i_k} = \sum_{j_1,j_2,...,j_k}^n ((\prod_{p=1}^k T_{i_p,j_p})\times t_{j_1,j_2,...,j_k})
$$

この計算をする関数を作る。
```{r}
my.Tensor.calc <- function(t,Z){
  n <- length(Z[1,])
	k <- log(length(t),n)

	if(k==0){
		return(t)
	}
	t <- as.array(t)
	ret <- rep(0,length(t))
	tmp.list <- list()
	dimt <- dim(t)
	for(i in 1:k){
		tmp.list[[i]] <- 1:dimt[i]
	}
	ad <- as.matrix(expand.grid(tmp.list))
	for(i in 1:length(t)){
		for(j in 1:length(t)){
			tmp <- 1
			for(kk in 1:k){
				tmp <- tmp * Z[ad[i,kk],ad[j,kk]]
			}
			ret[i] <- ret[i] + tmp * t[j]
		}
	}
	ret2 <- t
	ret2[1:length(t)] <- ret
	ret2
}
```

この関数を使って、ベクトルを回転してみる。
回転行列を使ったベクトルの回転と比較する。

```{r}
# 回転行列
n <- 3
angles <- rnorm(n) * 2*pi
Ax <- my.rot.3d.e(angles[1],angles[2],angles[3])

v <- my.random.complex(n)
Ax %*% v
my.Tensor.calc(v,Ax)
```
2つのベクトルの内積はスカラー。
ベクトルを回転してから、その内積を計算する方法と、スカラーを計算してから、それをテンソル変換する方法を並列で実施する。

実際、テンソルの計算では、スカラーは不変であることを前提に関数を書いてあるので、不変である。
```{r}
v1 <- my.random.complex(n)
v2 <- my.random.complex(n)
ipv12 <- my.ip.complex.v(v1,v2)
ipv12
my.ip.complex.v(Ax %*% v1, Ax %*% v2)
my.Tensor.calc(ipv12,Ax)
```

2階のテンソルとして、2ベクトルの要素の総当たりペアの積を考えてみる。要素数はn^2個である。

```{r}
prodv12 <- outer(v1,v2,"*")
prodv12
prodv12.rot <- outer(Ax %*% v1,Ax %*% v2, "*")
prodv12.rot
my.Tensor.calc(prodv12,Ax)
```

### 回転行列のサンドイッチでテンソル変換
2階のテンソルの場合は、変換行列とその逆行列でサンドイッチにして計算できる特徴がある。

```{r}
Ax %*% prodv12 %*% solve(Ax)
```

### determinantと置換・重複符号とテンソル
nxn行列のdeterminantはn階のテンソルの要素を、置換・重複符号を掛けて足し合わせたもの。

要素数は$n^n$通りある。重複がある場合は、置換・重複符号が0なので、それを除いて計算するのもよい。


```{r}
n <- 5
#M <- matrix(my.random.complex(n^2),n,n)
M <- matrix(rnorm(n^2),n,n)
library(gtools)
#p.list <- permutations(n,n,repeats.allowed=TRUE) # 全通り
p.list <- permutations(n,n) # 重複を除く
s.list <- apply(p.list,1,my.perm.sign)
det.out <- 0
for(i in 1:length(s.list)){
  det.out <- det.out + s.list[i] * prod(diag(M[p.list[i,],]))
}
det.out
my.det.complex(M)
```

実際、この交換・重複符号自体がn階テンソルになっていて、回転によって不変(なぜなら、回転行列のdeterminantが1であるから)。
それを確かめる。

```{r}
n <- 3
m <- 3
S <- array(0,rep(m,n))
dim(S)
ad <- as.matrix(expand.grid(rep(list(1:m),n)))
s <- apply(ad,1,my.perm.sign)
s
s. <- my.Tensor.calc(array(s,rep(m,n)),Ax)
round(c(s.),1)
```

### テンソル積もSO(n)

2階のテンソルに回転を施すとき、回転行列のクロネッカー積(行列のクロネッカー積はテンソル積の特殊版)2階のテンソルをベクトル化したものに施したものと同じ。

今、回転行列のクロネッカー積はやはり回転行列。
単位行列も逆行列もあり、群をなす。

実際、回転行列のk回のクロネッカー積が作る、(n^k)x(n^k)行列は、n次元空間における、k階テンソルの回転を表しており、SO(n)な群をなす。

nxn回転行列はベクトルに作用する回転行列で、(n^2)x(n^2)クロネッカー積行列は、2階のテンソルに作用する「回転行列」だということ。

k階テンソル積はSO(n)のn^k次元表現。

```{r}
Ax2 <- Ax %x% Ax
Ax2 %*% c(prodv12)
my.Tensor.calc(prodv12,Ax)
my.det.complex(Ax2)
```

### SU(2)の場合

3次元回転はSO(3)3次元回転はU2)3次元回転を表す3x3回転行列はベクトルに作用する。
2階のテンソルに作用するときは、回転行列で挟む。

2x2ユニタリ行列が回転を表すとき、「ベクトルに相当する」のは、2x2行列で、それに作用するときは、ユニタリ行列で挟む。

では、2x2ユニタリ行列が作用する複素ベクトル、というものがあるはず。

それがスピノール。

スピノールへの角$2\pi$の作用は符号反転を起こす。

### 4 = 3 + 1、3次元ベクトルと1個のスカラー

二次元複素ベクトルというスピノールの2階のテンソルが、複素2x2行列の形となって、回転を表すユニタリ行列が、サンドイッチとして作用する。

ただし、2x2複素行列の4要素のうち、3要素は3次元ベクトルで、残りの1要素はスカラーという混成になっている。

## リー代数

SO(3)やU(2)の元は行列であるが、今、ある行列Aの指数関数Bが、これらの群の元の行列であるとき、元の行列Aをその群のリー代数、またはリー環と言う。

行列の指数関数になっているので、「連続群」。

exp^(tA)のように、行列Aについて1パラメタtを用いることで、元が連続に軌跡を描く。この軌跡はリー代数の部分群。1パラメタで表せるので、1パラメタ部分群と呼ぶ。

リー代数の1パラメタ部分群は曲線であるから、それには接ベクトルが定義できる。

これは代数の幾何学的解釈をもたらす。

リー代数の部分群はみな、単位元を通る曲線。

単位元ではたくさんの曲線が交わっており、そこにはたくさんの接ベクトルがある。

この接ベクトルの連続変化に対応して、部分群の代表元も連続的に変化する。

リー代数はベクトルとしての性質を持つ。

リー代数では交換子積もリー代数の元になっている(交換子積についてリー代数は閉じている)。