[ruby-list:39317] Re: 固有値、固有ベクトルの計算

原です。

In [ruby-list:39288]

>佐藤と申します。
>
>行列の固有値、固有ベクトルを求めたいと
>思っています。
>
>webで調べてみた結果、下記の事がわかりましたが、
>調べもれ、事実誤認、その他の方法、
>実際に行列の計算をやっている人からの意見、ヒント
>などがありましたら教えてください。

多分佐藤さんの求めるものとは違うと思いますが，私の作った
ライブラリ algebra を紹介します。

[ruby-list:39310] の例で計算すると，

--------------------
require "algebra"
M = SquareMatrix(Rational, 3)

a = M[[1, 0, -1],
      [1, 2,  1], 
      [2, 2,  3]]

puts "A = "; a.display; puts
    
e = a.diagonalize
puts "特性多項式: #{e.chpoly} => #{e.facts}"; puts
puts "添加元:"
e.roots.each do |po, rs|
  puts "#{rs.join(', ')} : roots of #{po} == 0"
end; puts

puts "固有空間: "
e.evalues.uniq.each do |ev|
  puts "W_{#{ev}} = <#{e.espaces[ev].join(', ')}>"
end; puts

puts "変換行列:"
puts "P ="
e.tmatrix.display; puts

puts "P^-1 * A * P = "; 
(e.tmatrix.inverse * a * e.tmatrix).display; puts
--------------------

で，結果は

--------------------
A =
  1,   0,  -1
  1,   2,   1
  2,   2,   3

特性多項式: t^3 - 6t^2 + 11t - 6 => (t - 1)(t - 2)(t - 3)

添加元:

固有空間:
W_{1} = <[-1, 1, 0]>
W_{2} = <[-1, 1/2, 1]>
W_{3} = <[-1/2, 1/2, 1]>

変換行列:
P =
 -1,  -1, -1/2
  1, 1/2, 1/2
  0,   1,   1

P^-1 * A * P =
  1,   0,   0
  0,   2,   0
  0,   0,   3
--------------------

となります。

この algebra というのは，ピュア Ruby で書かれていて，代数的な
処理を目的としています（浮動小数点小数は念頭に置かれてない）。
従って得られる結果は厳密解で，理論的には何次の行列もＯＫな（対
角化不能の場合は例外が挙がる）のですが，著しく処理に時間がかか
る場合があります。特性多項式が因数分解できないなら，３次行列ぐ
らいが限界でしょうか。

例えば行列を

  1,   0,  -1
  1,   1,   1
  2,   2,   3

にすると，

--------------------
A =
  1,   0,  -1
  1,   1,   1
  2,   2,   3

特性多項式: t^3 - 5t^2 + 7t - 1 => t^3 - 5t^2 + 7t - 1

添加元:
a, b, -b - a + 5 : roots of t^3 - 5t^2 + 7t - 1 == 0

固有空間:
W_{a} = <[1/2a^2 - 2a + 3/2, -1/2a^2 + 5/2a - 3, 1]>
W_{b} = <[(-1/2a + 1/2)b - 1/2a^2 + 5/2a - 2, 1/2ab + 1/2a^2 - 5/2a + 1/2, 1]>
W_{-b - a + 5} = <[(1/2a - 1/2)b - 1/2a + 1/2, -1/2ab + 1/2, 1]>

変換行列:
P =
1/2a^2 - 2a + 3/2, (-1/2a + 1/2)b - 1/2a^2 + 5/2a - 2, (1/2a - 1/2)b - 1/2a + 1/2
-1/2a^2 + 5/2a - 3, 1/2ab + 1/2a^2 - 5/2a + 1/2, -1/2ab + 1/2
  1,   1,   1

P^-1 * A * P =
  a,   0,   0
  0,   b,   0
  0,   0, -b - a + 5
------------------------------

となって，無理数 a, b が現れます。そして a, b の大きさの
計量は全くしません。

というわけで，このライブラリは今のところ線形代数の試験問題を
自動生成させてやろうとか，よこしまなな用途にしかつかえません。
(^^;


In [ruby-list:39303]

> # 学校の授業では、学生にプログラムを書かせるより(車輪の再生産)
> # lapackなどのライブラリの存在と使い方を教えたほうが
> # 有意義だとつねづね思っています。

線形代数の授業としては，プログラムを書かせることは，自分がア
ルゴリズムを理解しているか確かめるために非常に有益だと思いま
すね。（実際にそんな時間はないけど。）教えている先生だって掃
き出し法が何重のループになるか答えられる人は少ないでしょう。
技術的なテーマの授業だとまた別でしょうが。


ちょっと大きめの行列を数値的に対角化したいなんて場合，手軽に
使えるライブラリは私もほしいですね。田中さんのを試してみよう。