- 球面に場があったとする
- それを球面調和関数で分解する
- 球面調和関数は直交基底関数なので、個々の関数が場を説明する成分を積分して計算してもよいが、SPHARMというツール(こちら)では、一般化逆行列を使っている(この論文のメソッドセクション)
- どういうことかというと、球面に多数の点があって、その点にある値
があるとする。
は球面上の点の位置を角度で表したもの
- 他方、この点には、球面調和関数ごとに値がある。もちろんそれらはの関数としても表せる
- さて。今、たくさんの点があったとき、個々の球面調和関数の重みが
だとしよう。
は球面調和関数のID番号のようなものである
- この問題を次のような推定問題として考えることにする
- ある点
には、従属変数の値
があるが、この従属変数の値が、球面場の値とする
- また、この点には、説明変数の列
がある。この説明変数の値を、このという球面上の点が持つ、球面調和関数の値とする
- このようにすることで、標本数＝球面上の点の数、説明変数の数＝用いる球面調和関数の数、という多変量解析問題とみなせる
- さて。説明変数(球面場)を、説明変数の一次線形和で表すというモデルにすることにする(球面調和関数同士は独立だから、このモデルで全く問題ない)
- 『ある意味での最もよい推定値』は線形代数として求めても良い
を解くだけ
とすることにして、
にムーアペンローズ一般化逆行列を使いましょう、と。- しかも、球面場の値(Y)を観測する点を固定しておくことにすれば、ムーアペンローズ逆行列は、１度だけ計算しておけば、後は、その逆行列を観察値ベクトルにかけるだけで、
は推定される
- RのMASSパッケージのginv()関数はsvd()特異値分解をしているし、複素行列にも対応していて、以下に示すように、複素行列のときに、ユニタリをとっていたりする。
- 一般化逆行列(Wiki)
> ginv
function (X, tol = sqrt(.Machine$double.eps))
{
if (length(dim(X)) > 2L || !(is.numeric(X) || is.complex(X)))
stop("'X' must be a numeric or complex matrix")
if (!is.matrix(X))
X <- as.matrix(X)
Xsvd <- svd(X)
if (is.complex(X))
Xsvd$u <- Conj(Xsvd$u)
Positive <- Xsvd$d > max(tol * Xsvd$d[1L], 0)
if (all(Positive))
Xsvd$v %*% (1/Xsvd$d * t(Xsvd$u))
else if (!any(Positive))
array(0, dim(X)[2L:1L])
else Xsvd$v[, Positive, drop = FALSE] %*% ((1/Xsvd$d[Positive]) *
t(Xsvd$u[, Positive, drop = FALSE]))
}
<bytecode: 0x00000000332756f0>
<environment: namespace:MASS>
- ちなみに、
として、正方行列の逆行列を計算するプロセスを使っても求められる。ginv()よりそのほうが速かった
library(MASS)
library(Matrix)
n <- 500
Z <- matrix(rnorm(n^2),n,n) + 1i * matrix(rnorm(n^2),n,n)
Zinv1 <- ginv(Z)
Zinv2 <- solve(t(Z)%*%Z) %*% t(Z)
range(Re(Zinv1-Zinv2))
range(Im(Zinv1-Zinv2))
> system.time(
+ for(i in 1:10){
+ Zinv1 <- ginv(Z)
+ }
+ )
ユーザ システム 経過
21.08 0.16 21.24
>
> system.time(
+ for(i in 1:10){
+ Zinv2 <- solve(t(Z)%*%Z) %*% t(Z)
+ }
+ )
ユーザ システム 経過
13.42 0.09 13.51
