2014-12-12

1 Quaternion Algebra ぱらぱらめくる『Quaternion Fourier Transforms for Signal and Image Processing』

Rのonionパッケージを使って、実際の処理などをなぞることとする(onionパッケージのごく基礎)
四元数
- 基本

library(onion)
a <- runif(4)
q <- a[1] * H1 + a[2] * Hi + a[3] * Hj + a[4] * Hk
q. <- a[1] * 1 + a[2] * 1i + a[3] *Hj + a[4] * Hk # これでもＯＫ
q == q.
qs <- c(H1,Hi,Hj,Hk)
lab <- c("1","i","j","k")
for(i in 1:4){
 for(j in 1:4){
  tmp <- paste(lab[i],lab[j],sep="*")
  print(tmp)
  print(qs[i]*qs[j])
 }
}

- 実部(S(q)、虚ベクトル部(V(q))の分離

q[[1]] # 実
c(q[[2]],q[[3]],q[[4]]) # 虚ベクトル
get.comp(q,1) # これでもＯＫ
get.comp(q,2)
Re(q) # これでも
Im(q)
i(q)
j(q)
k(q)

- 四元数の積の分解
  - - commutative partとnon-commutative part(cross-product部分)とに分けてある
  - さらに右辺第１項は[tex:\frac{qp+pq}{2}=S(q)S(p) - + S(q)V(p)+S(p)V(q)]となる
  - 第２項は、虚ベクトル部分のみで計算しても同じなので[tex:qp=S(q)S(p) - + S(q)V(p)+S(p)V(q)+V(q) \times V(p)]
  - これをRで書くと

p <- rquat(1) # 整数係数のランダム四減数を１個、発生
q <- rquat(1)
q
p
# <V(q),V(p)>は
-g.even(Im(q),Im(p))
# V(q) x V(p)は
g.odd(Im(q),Im(p))
q * p
g.even(q,p) + g.odd(q,p)
# 細分
Re(q) * Re(p) + g.even(Im(q),Im(p)) + Re(q) * Im(p) + Re(p) * Im(q) + g.odd(Im(q),Im(p))

- ノルムとNormed algebraであること

Norm(q)
Norm(p)
Norm(q*p)
Norm(q)*Norm(p)

- 共役四元数
  - 実部・虚部の取り出しと共役四元数
  - 共役と積の順序交換

q
Conj(q)
(q+Conj(q))/2
Re(q)
(q-Conj(q))/2
Im(q)
Conj(q*p)
Conj(p)*Conj(q)

- Involutions( を満足する関数)
  - Hi,Hj,Hkごとに定まるInvolutionのセットがCanonical involutions

q
-Hi * q * Hi
-(Hi * (-Hi * q * Hi) * Hi)
my.involution <- function(q,p){
 p. <- Im(p)
 p. <- p./sqrt(Norm(p.))
 -p. * q * p.
}
my.involution(my.involution(q,Hi),Hi)
my.involution(my.involution(q,Hj),Hj)
my.involution(my.involution(q,Hk),Hk)

- - 任意のノルム１の純虚四元数でInvolution

p
p. <- Im(p) # 純虚
p. <- p./sqrt(Norm(p.)) # ノルム 1
Norm(p.) # |p.|^2
q
my.involution(my.involution(q,p.),p.)

- - Involution いろいろ

p <- rquat(1) # 整数係数のランダム四減数を１個、発生
q <- rquat(1)
r <- rquat(1)
my.involution(q*p,r) - my.involution(q,r) * my.involution(p,r)
my.involution(my.involution(q,r),r) - q
my.involution(my.involution(q,Hi),Hj) - my.involution(q,Hk)

- 逆元

q^(-1)
Conj(q)/sqrt(Norm(q))
q * q^(-1)

- 除算
  - p/q はp*(q^(-1))か(q^(-1)) * pかの区別がつかないので注意
- 角座標的にunit quaternionとスカラーの積で表す

q
q. <- q/sqrt(Norm(q))
sqrt(Norm(q))
q.
Norm(q.)

- 正規直交基底(1,p,q,r=pq=qp)

n <- 3
library(GPArotation)
R <- Random.Start(n) # 3x3 正規直交基底
p <- q <- quaternion(1)
Im(p) <- R[,1]
Im(q) <- R[,2]
r <- p*q # 第３の基底ベクトルは他二つの積
r. <- q*p
g.even(p,q)
# g.even(p,q)=0なので、r = p*q = q*p
r
R[,3]

- 指数・対数
  - 複素数を長さと角とに分けたとき、角の方の指数は単位円周をぐるぐる回り、長さが実数の指数関数になったように、四元数でも単位四元数とそのスカラー倍に分けて扱う
  - $e^q = e^{S(q)} \times (\cos(|V(q)|)+ V_{st}(q) \sin(|V(q)|))$ 、ただし $V_{st}(q)$ はV(q)の単位純虚四元数
  - 対数は
    - 複素数の対数は、複素数で、実部と虚部とに分けられる。四元数のそれも、実部と虚部とに分ける。
    - $q = |q|e^{\mu \phi}=|q|(\cos\phi + \mu \sin\phi)$
    - $|q| = \sqrt{a^2+b^2+c^2+d^2}$
    - $\mu = \frac{b i + c j + d k}{\sqrt{b^2+c^2+d^2}}$
    - $\phi = arctan(\frac{\sqrt{b^2+c^2+d^2}}{a})$

# 指数
q <- rquat(1)
Sq <- Re(q)
Vq <- Im(q)
Vq.len <- sqrt(Norm(Vq))
V.st <- Vq/Vq.len
exp(Sq) * (cos(Vq.len) + V.st*sin(Vq.len))
# Rには実装されている
exp(q)
# 対数
log(q)
log(sqrt(Norm(q))) + log(q/sqrt(Norm(q)))
q
mu <- Im(q)/sqrt(Norm(Im(q)))
phi <- atan(sqrt(Norm(Im(q)))/Re(q))
sqrt(Norm(q)) * exp(phi * mu)
sqrt(Norm(q)) * (cos(phi) + mu * sin(phi))

- - 複素数の場合に、純虚数の指数関数は $e^{i\theta} = \cos(\theta) + i \sin{\theta}$ であるが、この虚数記号iに対応して、四元数の単位四元数 $\mu$ を持ち込んでいる形
- 四元数の表記法はいろいろある
  - オイラー
    - $q = |q|e^{\mu \phi}=|q|(\cos\phi + \mu \sin\phi)$
    - $|q| = \sqrt{a^2+b^2+c^2+d^2}$
    - $\mu = \frac{b i + c j + d k}{\sqrt{b^2+c^2+d^2}}$
    - $\phi = arctan(\frac{\sqrt{b^2+c^2+d^2}}{a})$
  - オイラー角
    - $q = |q|e^{\eta i}e^{\kappa j}e^{\phi k}$
  - Cayley-Dickson。二つの複素数に分解する
    - $q = z_1 + z_2 Hj$
    - $z_1 = a + ib,z_2 = c + id$ のようにiのみの虚部を持つ二つの複素数に分ける
  - Ortho-split/symplectic
    - $q=(a'+b' \mu) + (c' + d'\mu) \nu$
    - もしくは $q=(a'+b'\mu) + (c' - d' \mu) \nu$ のように、「共役複素数」的な分解もあって、これがCayley-Dicksonの「一般形」
    - i,j,kのような特定の軸に依存させずに $\mu$ 方向をとった上で、 $\mu$ と垂直な $\nu$ と併せてCayley-Dickson的に表す
    - ２つの直交純虚単位四元数 $\mu,\nu$ を取り、それによって分解する
    - $q = q_+ + q_-$
    - $q_+ = \frac{1}{2}(q+\mu q \nu)$
    - $q_- = \frac{1}{2}(q-\mu q \nu)$

q <- rquat(1)
# 二つの直交純虚単位四元数
R <- Random.Start(n) # 3x3 正規直交基底
m <- n <- quaternion(1)
Im(m) <- R[,1]
Im(n) <- R[,2]
q.p <- 1/2*(q+m*q*n)
q.n <- 1/2*(q-m*q*n)
# もちろん以下は0
q.p + q.n - q

- - Polar Cayley-Dicksonも
- Cj ペア記法
  - Ortho-split記法の軸をHjにして、四元数を複素数のペアで表すことにすると、四元数の演算が簡略化される。それは、Hi,Hj,Hkの間にあるぐるぐる回って、その３つ以外は現れない、という性質からくる
  - また、non-commutative であったはずの四元数積がcommutativeに扱えるようになる
  - $a + b i + c j + d k = (a+ c j) + (b + d j) i$ と分解して $(a+cj, b+dj)$ のペアを考えるということ
  - いくつかの計算例を以下に示す

my.Cj <- function(H){
 ret <- rep(0,2)
 ret[1] <- Re(H)+1i*j(H)
 ret[2] <- i(H)+1i*k(H)
 ret
}
q <- rquat(1)
q
q.Cj <- my.Cj(q)
q.Cj
my.Cj(Conj(q))
q.conj.Cj <- c(Conj(q.Cj[1]),-q.Cj[2])
q.conj.Cj
my.Cj(q^(-1))
c(Conj(q.Cj[1])/Norm(q),-q.Cj[2]/Norm(q))

- 行列表現
  - 4x4実行列表現

M1 <- diag(rep(1,4))
Mi <- Mj <- Mk <- matrix(0,4,4)
Mi[1,2] <- Mi[4,3] <- 1
Mi[2,1] <- Mi[3,4] <- -1
Mj[1,3] <- Mj[2,4] <- 1
Mj[3,1] <- Mj[4,2] <- -1
Mk[1,4] <- Mk[3,2] <- 1
Mk[2,3] <- Mk[4,1] <- -1

q <- rquat(1)
p <- rquat(1)

p.mat <- Re(p)*M1 + i(p)*Mi + j(p)*Mj + k(p)*Mk
q.mat <- Re(q)*M1 + i(q)*Mi + j(q)*Mj + k(q)*Mk

p*q
p.mat %*% q.mat

- - 2x2複素行列表現

Mc1 <- diag(rep(1,2))
Mci <- Mcj <- Mck <- matrix(0,2,2)
Mci[1,1] <- 1i
Mci[2,2] <- -1i
Mcj[1,2] <- -1
Mcj[2,1] <- 1
Mck[1,2] <- -1i
Mck[2,1] <- -1i

p.matc <- Re(p)*Mc1 + i(p)*Mci + j(p)*Mcj + k(p)*Mck
q.matc <- Re(q)*Mc1 + i(q)*Mci + j(q)*Mcj + k(q)*Mck

p*q
p.matc %*% q.matc

- べき乗
  - $q^n = |q|^n e^{n\mu \phi} = |q|(\cos(n\phi) + \mu \sin(n\phi))$