この記事はQiita 数学 Advent Calendar 2015の9日目の記事です。（後付け）
3行で頼む
この記事では、Rのouter()を外積代数に仕立て上げます。
特に、3次元の有名な外積（クロス積）を作ります。
もっと速い物を知っている方は教えてください。たのしめ！
outerとは何ぞ
Rにおけるouter
"outer"という言葉と、掛け算的な演算を想像しようとすると、いわゆる外積代数を思い浮かべてしまいます。
しかし、Rでいうouter(a,b)=a %o% bは、どうもそういう物では無いようです。
この間（テンソルの記事を書いた時）まで分かっていなかったのですが、これ、outerと言いつつ結局はテンソル的な演算をしているということなんですよね。
> a <- c(1,2,3)
> b <- c(4,5,6)
> a %o% b
[,1] [,2] [,3]
[1,] 4 5 6
[2,] 8 10 12
[3,] 12 15 18
どうも日本語の外積という言葉には、exterior productとouter productが対応しているようで、一般的な(?)外積はexterior productのようです。
ここでは、exterior productを再現する方法を考えます。
一般的な外積代数
もっとも有名な外積は3次元の場合で、3次元ベクトル×3次元ベクトル＝3次元ベクトルの演算だと思いますが、普通は次のような代数（と同型な代数）として定義されます。
定義：$V$を(体$k$上の)ベクトル空間とする。$V$のテンソル代数$T(V)$とし、$v\otimes v(v \in V)$によって生成される両側イデアルを$I$とする。$V$の外積代数とは、$\wedge (V)=T(V)/I$のことと定める。
※この定義と同型の意味で同値な定義が他にもいくつか知られていますが、割愛します。また、ベクトル空間よりも一般的なものでも定義できますが、これも割愛します。素朴な例だと$\mathbb{Z}$とか。
よく知られた3次元ベクトル×3次元ベクトル＝3次元ベクトルの演算は、$V=k^3$に対して、$\wedge V$の部分空間$V$と$V\otimes V/<v\otimes v>$を同一視することにより得られます。
%o%は、この定義の意味での外積代数とは、ちょっと意味合いが違います。
そもそも的な話で、例えば2次元ベクトルと4次元ベクトルみたいなものについても取れるので、Rのouterでいう「外積」とは何なのか、みたいなところもあるのですが。
Rで本当の外積代数を実現する
ナイーブな実装
そこで、本当の外積代数を実現することを考えます。$v \otimes v$で生成されるイデアルで割るというのは、要するに$v \otimes v = 0$という関係式を導入するのと一緒です。
本当のテンソル積に対して$v \otimes v = 0$は一般には正しくないので、$\wedge(V)$において$x \otimes y$に対応する要素を$x\wedge y$と書くことにすれば、$(u+v)\wedge(u+v)=0$から$u\wedge v + v\wedge u = 0$が出ます。いわゆる交代性ですね。
特に$V$の適当な基底を$e_1,...,e_n$とすると、$e_2\wedge e_1 = -e_1\wedge e_2$といった具合に、積の順序の入れ替えによって符号が反転していきます。
まずは、この部分を計算する為の道具を用意します。
# Permutation(size,d)をリスト形式で返す sPd 的な
# combn(c(1,2,3),2,simplify=F) と同様の形式
permlist <- function(size, p = c(), l = list(), d = size){
if(d == 0){
l[[length(l) + 1]] = p
} else {
for(i in 1:size){
if(!(i %in% c(p))){
l <- permlist(size, c(p, i), l, d - 1)
}
}
}
l
}
# 転倒数によって偶奇を定める
odd.perm <- function(ls){
oe <- 0
n = length(ls)
for(i in 1:n){
oe <- oe + sum(ls[i:n] < ls[i])
}
oe %% 2
}
permlistは並び変えたものをリストで返却します。これが積の順序の交換に対応します。
また、このpermlistの各要素に対して、上のような入れ替えを行った時の正負を求める為の関数がodd.permです。
0が正、1が負の場合に対応します。
この道具を用いて、次のように計算をします。
# 行列バージョン
outer.array <- function(size, d = size, arr = array(1,dim=rep(size,d))){
mat <- array(0,dim=rep(size,d))
l <- permlist(size,c(),list(),d)
for(i in 1:length(l)){
n <- Reduce(function(a,b){a + (b-1)*size}, l[[i]], right = T)
mat[n] = (-1)^odd.perm(l[[i]])
}
mat*arr
}
# 適当な基底に対するベクトル表示バージョン
outer.base <- function(size, d = size, basis = combn(1:size,d,simplify=F), arr){
r <- c()
l <- permlist(size,c(),list(),d)
for(j in 1:length(basis)){
o <- (-1)^odd.perm(basis[[j]])
mat <- array(0,dim=rep(size,d))
for(i in 1:length(l)){
n <- Reduce(function(a,b){a + (b-1)*size}, l[[i]], right = T)
mat[n] = (-1)^odd.perm(l[[i]])*o*(Reduce(function(a,b){a&b}, sort(l[[i]])==sort(basis[[j]])))
}
r[j] <- sum(mat*arr)
}
r
}
# いわゆる普通の外積
outer.3x3 <- function(l, r){
b <- list()
b[[1]] = c(2,3)
b[[2]] = c(3,1)
b[[3]] = c(1,2)
arr <- l %o% r
outer.base(3, 2, basis=b, arr=arr)
}
はじめの関数は、配列(array)の形式のままで、例えば$e_i\wedge e_i$のような消えてしまう成分を0として、生き残る成分のみを、正負付で計算するものです。
次の関数は、これらの成分について、適当な基底に対して和をとり、各基底の係数を返す関数です。
これを具体的に使用して、一般的な3次元の外積（クロス積）を定義しているのが、outer.3x3です。
それぞれを試しに使ってみると、次のようになります。
おためし
> # 普通の外積の検算
> outer.3x3(c(1,0,0), c(0,1,0))
[1] 0 0 1
> outer.3x3(c(0,1,0), c(0,0,1))
[1] 1 0 0
> outer.3x3(c(0,0,1), c(1,0,0))
[1] 0 1 0
> outer.3x3(c(1,1,0), c(-1,1,0))
[1] 0 0 2
>
> # 係数を確認
> outer.array(3,2)
[,1] [,2] [,3]
[1,] 0 1 1
[2,] -1 0 1
[3,] -1 -1 0
>
> # x <- array(sample(100,9), c(3,3))
> x <- array(c(27,15,37,83,61,56,76,70,75),c(3,3))
> det(x)
[1] 31588
>
> arr <- x[1,] %o% x[2,] %o% x[3,]
> outer.base(3,arr=arr)
[1] 31588
>
> arr <- x[,1] %o% x[,2] %o% x[,3]
> outer.base(3,arr=arr)
[1] 31588
>
> # x <- array(sample(100,16), c(4,4))
> x <- array(c(45,19,91,11,31,82,16,12,1,10,32,50,84,43,97,74),c(4,4))
> det(x)
[1] 16317636
>
> arr <- x[1,] %o% x[2,] %o% x[3,] %o% x[4,]
> outer.base(4,arr=arr)
[1] 16317636
>
> arr <- x[,1] %o% x[,2] %o% x[,3] %o% x[,4]
> outer.base(4,arr=arr)
[1] 16317636
もっと速くて便利なものはないかしら
ここまでの話は、教育的観点やRの仕様についての理解を深める意味では有意義だと思いますが、繰り返し頻繁に使うコードとしてはたぶん遅いと思います。
（低次元でちょろちょろ試すぐらいのレベルでは十分速いですが）
きっと、もっと速い便利なパッケージがあるのではないでしょうか。（願望）
うんちく
テンソル代数$T(V)$は、積の「親玉」というか、もっとも普遍的なもので、外積代数$\wedge(V)$は交代性を持つ積の中でもっとも普遍的なものです。じつは、これと同じように積が可換な場合、つまり$x\otimes y-y\otimes x=0$を条件にすると、類似の代数を構成できます。
これは対称代数という名前がついていますが、$V$が体$k$上のベクトル空間であれば、要するに多項式環のことです。(同型の意味で)
おわり
たのしめ！