行列の固有値を求めたい。
なにをいまさら感たっぷりなんだけど、理系の出版社でテクニカルレビューっぽいことをしていると、特に統計の本なんかで行列の固有値を検算したい場合がなきにしもあらず。たいていは Maxima に突っ込んじゃえばいいんだけど、Maxima では特性多項式を作ってそれを解いているらしく、うまく解が求まりやがらないことがある。実際、maxima/share/matrix/eigen.macを見ると、charpolyという関数で特定多項式を作り、それをsolveという関数で解いている。で、この solve という多項式を解く関数がうまくないといったことがMaximaのマニュアルにも書いてある。

eigenvalues calls the function solve to find the roots of the characteristic polynomial of the matrix. Sometimes solve may not be able to find the roots of the polynomial

幸い、固有値を求めたい行列は対称行列ばっかりなので、自分でJacobi法をナイーブに実装してみた。対称行列だとうれしい理由やJacobi法については『プログラミングのための線形代数』の第5章で。

jacobi.scm

gosh> 
(define A
  ((make-matrix 6 6)
   0.68 0.65 0.8 0.11 0.01 0.14 0.65 0.88 0.89 0.02 0.19 0.01 0.8 0.89 0.91 0.02 0.04 0.1 0.11
   0.02 0.02 0.81 0.82 0.77 0.01 0.19 0.04 0.82 0.81 0.64 0.14 0.01 0.1 0.77 0.64 0.66))

gosh> A
((0.68 0.65 0.8 0.11 0.01 0.14) (0.65 0.88 0.89 0.02 0.19 0.01) (0.8 0.89 0.91 0.02 0.04 0.1) 
(0.11 0.02 0.02 0.81 0.82 0.77) (0.01 0.19 0.04 0.82 0.81 0.64) (0.14 0.01 0.1 0.77 0.64 0.66))

gosh> (eigenvalues-jacobi A)
(-0.010575430595377565 -0.09400171282374802 2.109504119397093 -0.08138275206455295 0.279213041557744 2.5472427345288424)

行列のデータをどういうふうに表現すべきなのか、ひどく悩んだ。Dybvig 本のようにベクトルを使うのがセオリーなんだろう。ベクトルなら行列のij要素を取り出したり変更したりするのも簡単だ。でも、行と列を順繰りになめて各要素をちこちこ掛けたり足したりする処理ばっかり書くことになりそうで、ちょとブルー。そんなわけで今回はあえてリストのリストとして定義した（行ベクトルが1つのリスト。それを要素に持つ縦ベクトルのリストが行列）。リストなので、行列の積や和なんかも高階関数だけで定義できる。

(define (matrix-transform m)
  (apply zip m))
(define (matrix-product m1 m2)
  (if (not (= (line-length m1) (row-length m2)))
      (error "Matrix-Product Not Defined" (list m1 m2))
      (apply (make-matrix (line-length m1) (row-length m2))
             (map (cut fold + 0 <>)
                  (map (cut apply map * <>)
                       (cartesian-product (list m1 (matrix-transform m2))))))))
(define (matrix-sum m1 . ms)
  (map (cut apply map + <>)
       (apply zip m1 ms)))

ちょっと満足。ただし単位行列や回転行列を定義するのに要素を頭からなめるはめになってるんだけどね。満足しちゃだめ！