VBScriptで固有値・固有ベクトルを求める (反復法)

n × n の正方行列 A と
n次元のベクトル x について
Ax = λx (ただし x ≠ 0) が成り立つとき
λを固有値, x を固有ベクトルという.
最初に適当なベクトルx₀から始めて x_k+1 = Ax_k を反復すると
x_k は行列 A の最大固有値に対応する固有ベクトルに収束する.
固有値はレイリー(Rayleigh)商
λ = (x_k^Tx_k+1) / (x_k^Tx_k)
により求める.
べき乗法、累乗法とも言う.

Option Explicit

Private Const N = 3

'ベキ乗法で最大固有値を求める
Call Main

Private Sub Main()
    Dim a: a = Array(_
               Array(5.0, 4.0, 1.0, 1.0), _
               Array(4.0, 5.0, 1.0, 1.0), _
               Array(1.0, 1.0, 4.0, 2.0), _
               Array(1.0, 1.0, 2.0, 4.0))

    Dim x: x = Array(1.0 ,0.0 ,0.0 ,0.0)

    'ベキ乗法
    Dim lambda: lambda = power(a, x)

    WScript.StdOut.WriteLine
    WScript.StdOut.WriteLine "eigenvalue"
    WScript.StdOut.WriteLine Right(Space(14) & FormatNumber(lambda, 10, -1, 0, 0), 14)

    WScript.StdOut.WriteLine "eigenvector"
    Call disp_vector(x)
End Sub

'ベキ乗法
Private Function power(ByVal a, ByRef x0)
    Dim i, j, k
    
    Dim lambda: lambda = 0.0

    '正規化 (ベクトル x0 の長さを１にする)
    Call normarize(x0)
    Dim e0: e0 = 0.0
    For i = 0 To N
        e0 = e0 + x0(i)
    Next

    For k = 1 To 200
        '１次元配列を表示
        WScript.StdOut.Write Right(Space(3) & k, 3) & vbTab
        Call disp_vector(x0)

        '行列の積 x1 = A × x0 
        Dim x1: x1 = Array(0.0, 0.0, 0.0, 0.0)
        For i = 0 To N
            For j = 0 To N
                x1(i) = x1(i) + a(i)(j) * x0(j)
            Next
        Next

        '内積
        Dim p0: p0 = 0.0
        Dim p1: p1 = 0.0
        For i = 0 To N
            p0 = p0 + x1(i) * x1(i)
            p1 = p1 + x1(i) * x0(i)
        Next
        '固有値
        lambda = p0 / p1

        '正規化 (ベクトル x1 の長さを１にする)
        Call normarize(x1)
        '収束判定
        Dim e1: e1 = 0.0
        For i = 0 To N
            e1 = e1 + x1(i)
        Next
        If Abs(e1 - e0) < 0.00000000001 Then Exit For

        x0 = x1
        e0 = e1
    Next

    power = lambda
End Function

'１次元配列を表示
Private Sub disp_vector(ByVal row)
    Dim col
    For Each col In row
        WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab
    Next
    WScript.StdOut.WriteLine
End Sub

'正規化 (ベクトルの長さを１にする)
Private Sub normarize(ByRef x())
    Dim s: s = 0.0

    Dim i
    For i = 0 To N
        s = s + x(i) * x(i)
    Next
    s = Sqr(s)
    
    For i = 0 To N
        x(i) = x(i) / s
    Next
End Sub

z:\vbs>cscript //nologo 1101.vbs
  1       1.0000000000    0.0000000000    0.0000000000    0.0000000000
  2       0.7624928517    0.6099942813    0.1524985703    0.1524985703
  3       0.6745979924    0.6589096670    0.2353248811    0.2353248811
  4       0.6517382447    0.6501601860    0.2761602732    0.2761602732
  5       0.6421032522    0.6419452154    0.2963188771    0.2963188771
  6       0.6373267026    0.6373108932    0.3063082594    0.3063082594
  7       0.6349074765    0.6349058954    0.3112772079    0.3112772079
  8       0.6336860412    0.6336858831    0.3137548428    0.3137548428
  9       0.6330719648    0.6330719490    0.3149919005    0.3149919005
 10       0.6327640473    0.6327640457    0.3156099832    0.3156099832
 11       0.6326098648    0.6326098646    0.3159189122    0.3159189122
 12       0.6325327172    0.6325327172    0.3160733485    0.3160733485
 13       0.6324941293    0.6324941293    0.3161505596    0.3161505596
 14       0.6324748319    0.6324748319    0.3161891634    0.3161891634
 15       0.6324651822    0.6324651822    0.3162084649    0.3162084649
 16       0.6324603572    0.6324603572    0.3162181155    0.3162181155
 17       0.6324579446    0.6324579446    0.3162229408    0.3162229408
 18       0.6324567383    0.6324567383    0.3162253534    0.3162253534
 19       0.6324561352    0.6324561352    0.3162265597    0.3162265597
 20       0.6324558336    0.6324558336    0.3162271629    0.3162271629
 21       0.6324556828    0.6324556828    0.3162274644    0.3162274644
 22       0.6324556074    0.6324556074    0.3162276152    0.3162276152
 23       0.6324555697    0.6324555697    0.3162276906    0.3162276906
 24       0.6324555509    0.6324555509    0.3162277283    0.3162277283
 25       0.6324555415    0.6324555415    0.3162277472    0.3162277472
 26       0.6324555367    0.6324555367    0.3162277566    0.3162277566
 27       0.6324555344    0.6324555344    0.3162277613    0.3162277613
 28       0.6324555332    0.6324555332    0.3162277637    0.3162277637
 29       0.6324555326    0.6324555326    0.3162277648    0.3162277648
 30       0.6324555323    0.6324555323    0.3162277654    0.3162277654
 31       0.6324555322    0.6324555322    0.3162277657    0.3162277657
 32       0.6324555321    0.6324555321    0.3162277659    0.3162277659
 33       0.6324555321    0.6324555321    0.3162277659    0.3162277659
 34       0.6324555321    0.6324555321    0.3162277660    0.3162277660
 35       0.6324555320    0.6324555320    0.3162277660    0.3162277660

eigenvalue
 10.0000000000
eigenvector
  0.6324555320    0.6324555320    0.3162277660    0.3162277660

参考文献