球面調和関数係数も一般化逆行列で推定

  • 球面に場があったとする
  • それを球面調和関数で分解する
  • 球面調和関数は直交基底関数なので、個々の関数が場を説明する成分を積分して計算してもよいが、SPHARMというツール(こちら)では、一般化逆行列を使っている(この論文のメソッドセクション)
  • どういうことかというと、球面に多数の点があって、その点にある値v(\theta,\phi)があるとする。\theta,\phiは球面上の点の位置を角度で表したもの
  • 他方、この点には、球面調和関数ごとに値がある。もちろんそれらは\theta,\phiの関数としても表せる
  • さて。今、たくさんの点があったとき、個々の球面調和関数の重みがw_l^mだとしよう。l,mは球面調和関数のID番号のようなものである
  • この問題を次のような推定問題として考えることにする
  • ある点V_iには、従属変数の値y_iがあるが、この従属変数の値が、球面場の値とする
  • また、この点V_iには、説明変数の列X_iがある。この説明変数の値を、このV_iという球面上の点が持つ、球面調和関数の値とする
  • このようにすることで、標本数=球面上の点の数、説明変数の数=用いる球面調和関数の数、という多変量解析問題とみなせる
  • さて。説明変数(球面場)を、説明変数の一次線形和で表すというモデルにすることにする(球面調和関数同士は独立だから、このモデルで全く問題ない)
  • 『ある意味での最もよい推定値』は線形代数として求めても良い
  • Y = Xbを解くだけ
  • b=X^{-1}Yとすることにして、X^{-1}にムーアペンローズ一般化逆行列を使いましょう、と。
  • しかも、球面場の値(Y)を観測する点を固定しておくことにすれば、ムーアペンローズ逆行列は、1度だけ計算しておけば、後は、その逆行列を観察値ベクトルにかけるだけで、bは推定される
  • 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>
  • ちなみに、(Z^t Z)^{-1} Z^tとして、正方行列の逆行列を計算するプロセスを使っても求められる。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