ちょっとずらして解く - ryamadaのコンピュータ・数学メモ

$M x = 0$ (Mは正方行列、xはベクトル)を解くと、Mのdeterminantが0出なければ、xは0ベクトルになるけれど、それが欲しいわけじゃなくて、 $M x=0$ をだいたい満足するxが欲しいときの話
こちらにあるようにWillmoreフローの計算機解の一環として、次のような問題がある
- $(D-\rho) \lambda=0$ を解きたい。ただし、 $D$ は行列、 $\rho$ は実数列で、ここでは、対角成分に相当する(また、 $\lambda$ が四元数なので、四元数を実＋３虚に引き伸ばして行列・ベクトル化している、などもあり、「実相当部分の対角成分」というのが正確なところだが、ここでは、四元数のことは忘れて、複素数で考えることにする)。 $\lambda$ は四元数のベクトル(だが、ここでは普通に複素数までとする)
このままだと、に0ベクトルという解が出てしまうので、次のようにする
- $(D-\rho) \lambda = \gamma \lambda$ と固有値・固有ベクトル問題にして、 $(D-(\rho+\gamma))\lambda = 0$ の解で代用する
- これは、ある量 $\rho$ を各変数ごとに定めて動かすとして、それにうまく合致する $\lambda$ を求めたかったけれど、それだと解は出ないので、各変数ごとに動かす量として、定めた量に定数を加減してやったもので代用してよいのなら、解が出る、という話
- ここで得られる解は、固有ベクトルなので、その定数倍はすべて「解」になる。今、代用しているわけなので、誤差を含むから、この固有ベクトルを小さく定数倍して変化を定め、また、次の時点で同じことをやり直す、という作戦をとるべし
固有値は複数出るだろうけれど、そのうちの最小のものを取っておけば、結構いい感じの代用結果になるはず
やってみる

d <- 3

D <- matrix(rnorm(d^2),d,d)
rho <- rnorm(d)

D.rho <- D-diag(rho,d)


# (D-rho)lambrda = 0 はゼロベクトルが答えになって困る
solve(D.rho,rep(0,d))

e.out <- eigen(D.rho)
e.val <- e.out[[1]]
e.vec <- e.out[[2]]


D %*% e.out[[2]]


for(i in 1:d){
	D.rho2 <- D-diag(rho+e.val[i],d)
	print(D.rho2 %*% e.vec[,i])
}

              [,1]
[1,]  2.220446e-16
[2,] -1.339207e-15
[3,] -3.330669e-16
              [,1]
[1,]  3.330669e-16
[2,] -3.191891e-16
[3,]  3.885781e-16
              [,1]
[1,] -2.775558e-16
[2,] -4.440892e-16
[3,] -1.665335e-16