フォン・ノイマン環と射影行列 - ryamadaのコンピュータ・数学メモ

フォン・ノイマン環 Von Neumann algebra - Wikipediaを勉強していると、射影行列というものが出てくる。
その環の要素の説明に射影行列が出てくる。
であると言う。
- これは、量子力学では状態の性質(観測すると0または1が得られる。0は存在しない、1は存在する)に対応する作用素。その固有値は0と1
射影行列は線形回帰の文脈でも登場する。Projection matrix - Wikipedia
最小２乗法により、観測値セットを表す点から、モデル部分空間への垂線の足が推測値セットを表す点になることを言っている。
WikipediaVon Neumann algebraの記事には、von Neumann algebraとProjectionsとの関連として、"they (projections) are exactly the operators which give an orthogonal projection of H onto some closed subspace."と記載されているように、射影行列は空間を亜空間に直交射影する。
行列を考えているから、この射影は線形射影である。
空間のすべての点が(線形)部分空間に射影されるから、射影行列を構成するベクトルは線形独立ではなく、固有値の中に0が含まれる～行列式は0である。
２ｘ２行列でなる行列がどうなるか考えてみる。
- $E=E^*$ ： $*$ は、共役転置のことなので、２ｘ２行列は $\begin{pmatrix} a, b + i c \\ b- ic, d \end{pmatrix},a,b,c,d \in R$ のように、対角成分は実数、２つの非対角成分は相互に共役複素数となる。
- これを書き換えると $\begin{pmatrix} a, r e^{it} \\ r e^{-it}, d \end{pmatrix},a,d,r,t \in R$ となる。
- $E = EE$ という条件から、 $(\begin{pmatrix} a, r e^{it} \\ r e^{-it}, d \end{pmatrix})^2 = \begin{pmatrix} a^2 + r^2, (a+d)r e^{it} \\(a+d)r e^{-it}, d^2 + r^2 \end{pmatrix} = \begin{pmatrix} a, r e^{it} \\ r e^{-it}, d \end{pmatrix}$ となる。
- $r=0$ の場合には、 $a=d=1$ となり、 $E= I$ (単位行列)となる。
- の場合には、から、が得られる。また、から。
  - ただし、 $a(1-a) > 0$ のときは $r = \sqrt{a(1-a)} \in R$ となるが、 $a(1-a) < 0$ のときは、 $r = \sqrt{a(1-a)} \in C$ (純虚数)になっていることに注意
- したがって、 $\begin{pmatrix} a, \sqrt{a(1-a)} e^{it}\\ \sqrt{a(1-a)}e^{-it},1-a \end{pmatrix}, a,t \in R$ となる。
- 固有値は、1と0になっている。
- 固有値１に対応する固有ベクトルで、かつ、二乗ノルムが１なるベクトルは、 $e^{i \phi}(\sqrt{a} e^{it},\sqrt{1-a})$ を二乗ノルムが１になるように補正したベクトル、固有値0に対応する固有ベクトルは $e^{i \phi}(\sqrt{1-a}e^{it},\sqrt{a})$ を二乗ノルムが１になるように補正したベクトル
- 0なる固有値があることが、２次元平面を直線につぶすことを意味し、1なる固有値があることが、固有値1に対応する固有ベクトルが表す直線上の点は動かないことを表しており、その直線と交叉し0なる固有値に対応する固有ベクトル方向の直線は、その交点につぶされることを意味する。
この射影行列が次のように分解されるという話がある。。このとき、ととはequivalentという話がある。これを２ｘ２の場合について、a,tでパラメタライズした場合との関係を考えてみる。
- ２ｘ２複素行列を $u = (r_{i,j} e^{i \theta_{i,j}} ), i,j \in \{1,2\}$ とする。
- $u u^* = \begin{pmatrix} r_{1,1}^2 + r_{1,2}^2, r_{1,1}r_{2,1} e^{i (\theta_{1,1}-\theta_{2,1})} + r_{1,2} r_{2,2} e^{i (\theta_{1,2} -\theta_{2,2})} \\ r_{1,1} r_{2,1} e^{-i (\theta_{1,1}-\theta_{2,1})} + r_{1,2} r_{2,2} e^{-i (\theta_{1,2} - \theta_{2,2})}, r_{2,1}^2 + r_{2,2}^2 \end{pmatrix}$
- $u u^* = E = \begin{pmatrix} a, (a(1-a))^{0.5} e^{it} \\ (a(1-a))^{0.5} e^{-i t}, 1-a \end{pmatrix}$ であるから
- $r_{1,1}^2 + r_{1,2}^2 = a, r_{2,1}^2 + r_{2,2}^2 = 1-a$
- - $a(1-a) =(r_{1,1}^2 + r_{1,2}^2 )(r_{2,1}^2 + r_{2,2}^2)$ を右辺に代入して、式変形すると。。。
  - $r_{1,1}^2 r_{2,2}^2 + r_{1,2}^2 r_{2,1}^2 - 2 r_{1,1} r_{1,2} r_{2,1} r_{2,2} \cos{\theta_{1,1} - \theta_{1,2}- \theta_{2,1} + \theta_{2,2}} = 0$ が得られる。さらに変形して
  - $(r_{1,1}r_{2,2} - r_{1,2} r_{2,1})^2 = 2 r_{1,1}r_{1,2}r_{2,1}r_{2,2} (\cos{\theta_{1,1} - \theta_{1,2}- \theta_{2,1} + \theta_{2,2}} - 1)$ となる。
  - ここから、 $r_{1,1}r_{2,2} - r_{1,2} r_{2,1} = 0$ , $\theta_{1,1} - \theta_{1,2}- \theta_{2,1} + \theta_{2,2} = 0$ が得られる。
  - $r$ と $\theta$ との制約は分離されていて次のような制約にまとめられる。
  - $r_{1,2} = \pm (a- r_{1,1})^{0.5}$ , $r_{2,1} = \pm ((1-a)/a)^{0.5} r_{1,1}$ , $r_{2,2} = \frac{r_{1,2} r_{2,1}}{r_{1,1}}$
  - $\theta_{1,1} - \theta_{1,2}- \theta_{2,1} + \theta_{2,2} = 0$
- 射影行列のことを別の観点から考えてみる
  - 今、ある状態確率ベクトル $\psi$ （複素ベクトルであって、その要素のModの二乗和が１）が存在しているとする。その状態の存否を0,1で判定したくなる
  - $| \psi > < \psi$ なる行列は射影行列になっている。そして、その行列を状態ベクトル $\psi$ に作用すると $\psi$ になる（そして $\psi$ と作用してできた $\psi$ の積の二乗ノルムが１になるという意味で「存否＝１」)。また、同じことだが、この射影行列の固有値は1と0で、固有値1に対する固有ベクトルは $\psi$ と（位相違いはあるが）同じベクトルになっている(こちらに、量子計算におけるProjectorオペレータの定義が出ている)

v <- runif(3) + 1i * runif(3)
v. <- v/sqrt(sum(Mod(v)^2))
sum(Mod(v.)^2) # 1
v.m <- matrix(v.,ncol=1)
p <- v.m %*% t(Conj(v.m))
eigen(p) -> eigen.out
v.m / eigen.out[[2]][,1]
v.m / eigen.out[[2]][,1] -> a
Mod(a[1])
p %*% v.m -> q
q/v.m

a <- 0.3
t <- 0.2
r <- sqrt(a*(1-a) + 0 * 1i)
M <- matrix(c(a, r * exp(1i * t), r * exp(-1i * t), 1-a),byrow=TRUE,2,2)
M

M %*% M
round(M%*% M - M) # 0

Conj(t(M))
round(Conj(t(M)) - M) # 0
eout <- eigen(M)
round(eout[[1]]) # 固有値は1と0

ev1 <- c((a+0*1i)^0.5 * exp(1i *t),(1-a+0*1i)^0.5)
ev2 <- c(-(1-a+0*1i)^0.5 * exp(1i * t), (a+0*1i)^0.5)

ev1 <- ev1 / sqrt(sum(Mod(ev1)^2))
ev2 <- ev2 / sqrt(sum(Mod(ev2)^2))

sum(Mod(ev1)^2)
sum(Mod(ev2)^2)
sum(Mod(eout[[2]][,1])^2)
sum(Mod(eout[[2]][,2])^2)
eout[[2]][,1]/ev1
Mod(eout[[2]][,1]/ev1)
cbind(ev1,ev2) / eout[[2]]
Mod(cbind(ev1,ev2) /eout[[2]])

> round(M%*% M - M) # 0
     [,1] [,2]
[1,] 0+0i 0+0i
[2,] 0+0i 0+0i
> round(Conj(t(M)) - M) # 0
     [,1] [,2]
[1,] 0+0i 0+0i
[2,] 0+0i 0+0i
> round(eout[[1]]) # 固有値は1と0
[1] 1 0
> sum(Mod(ev1)^2)
[1] 1
> sum(Mod(ev2)^2)
[1] 1
> sum(Mod(eout[[2]][,1])^2)
[1] 1
> sum(Mod(eout[[2]][,2])^2)
[1] 1
> eout[[2]][,1]/ev1
[1] -0.9800666+0.1986693i -0.9800666+0.1986693i
> Mod(eout[[2]][,1]/ev1)
[1] 1 1
> cbind(ev1,ev2) / eout[[2]]
                       ev1                  ev2
[1,] -0.9800666-0.1986693i 0.9800666+0.1986693i
[2,] -0.9800666-0.1986693i 0.9800666+0.1986693i
> Mod(cbind(ev1,ev2) /eout[[2]])
     ev1 ev2
[1,]   1   1
[2,]   1   1

a <- -0.3
t <- 0.2
r <- sqrt(a*(1-a) + 0 * 1i)
M <- matrix(c(a, r * exp(1i * t), r * exp(-1i * t), 1-a),byrow=TRUE,2,2)
M

M %*% M
round(M%*% M - M) # 0

Conj(t(M))
round(Conj(t(M)) - M) # 0
eout <- eigen(M)
round(eout[[1]]) # 固有値は1と0



ev1 <- c((a+0*1i)^0.5 * exp(1i *t),(1-a+0*1i)^0.5)
ev2 <- c(-(1-a+0*1i)^0.5 * exp(1i * t), (a+0*1i)^0.5)

ev1 <- ev1 / sqrt(sum(Mod(ev1)^2))
ev2 <- ev2 / sqrt(sum(Mod(ev2)^2))

sum(Mod(ev1)^2)
sum(Mod(ev2)^2)
sum(Mod(eout[[2]][,1])^2)
sum(Mod(eout[[2]][,2])^2)
eout[[2]][,1]/ev1
Mod(eout[[2]][,1]/ev1)
cbind(ev1,ev2) / eout[[2]]
Mod(cbind(ev1,ev2) /eout[[2]])

> round(M%*% M - M) # 0
     [,1] [,2]
[1,] 0+0i 0+0i
[2,] 0+0i 0+0i
> round(Conj(t(M)) - M) # 0
     [,1] [,2]
[1,] 0+0i 0-1i
[2,] 0-1i 0+0i
> round(eout[[1]]) # 固有値は1と0
[1]  1.000000e+00+0i -2.220446e-16+0i
> sum(Mod(ev1)^2)
[1] 1
> sum(Mod(ev2)^2)
[1] 1
> sum(Mod(eout[[2]][,1])^2)
[1] 1
> sum(Mod(eout[[2]][,2])^2)
[1] 1
> eout[[2]][,1]/ev1
[1] 1+0i 1+0i
> Mod(eout[[2]][,1]/ev1)
[1] 1 1
> cbind(ev1,ev2) / eout[[2]]
      ev1                   ev2
[1,] 1-0i -0.9800666-0.1986693i
[2,] 1+0i -0.9800666-0.1986693i
> Mod(cbind(ev1,ev2) /eout[[2]])
     ev1 ev2
[1,]   1   1
[2,]   1   1

$p = u u^*$ の確認

a <- 0.25 + 0 * 1i
r1 <- rep(0.1 + 0 * 1i,4)

r2 <- c((a-r1^2)^0.5,(a-r1^2)^0.5,-(a-r1^2)^0.5,-(a-r1^2)^0.5)
r3 <- c(((1-a)/a)^0.5 * r1,-((1-a)/a)^0.5 * r1,((1-a)/a)^0.5 * r1,-((1-a)/a)^0.5 * r1)
r4 <- r2 * r3 / r1

theta1 <- runif(1)
theta2 <- runif(1)
theta3 <- runif(1)
theta4 <- theta2 + theta3 - theta1

for(i in 1:4){
	u <- matrix(c(r1[i] * exp(1i * theta1), r2[i] * exp(1i * theta2), r3[i] * exp(1i * theta3), r4[i] * exp(1i * theta4)), byrow=TRUE,2,2)

	ustar <- ustar <- t(Conj(u))

	p <- u %*% ustar

	q <- ustar %*% u

	print(p - t(Conj(p))) # 0 : p = p*
	print(round(p %*% p - p)) # 0 : p^2 = p
	print(q - t(Conj(q))) # 0 : q = q*
	print(round(q %*% q - q)) # 0 : q^2 = q

}