対称行列の固有値・固有ベクトルを求めるマクロ

エクセルで主成分分析をしたくて固有値・固有ベクトルを求めるマクロを作った．

参考
標準モジュール
クラスモジュール Matrix
クラスモジュール SymmetryMatrix
クラスモジュール GivensMatrix

参考

標準モジュール

グローバル変数

Private x As SymmetryMatrix '最初は元のデータ行列。最終的に対角行列になる。
Private u As Matrix  '固有ベクトルを並べた直交行列

Const CONVERGENCE_CRITERIA As Double = 10 ^ -20 '収束判定基準
Const MAX_ITERATION As Long = 1000   '最大反復回数

ワークシート関数

固有値を並べたベクトルを出力する関数EigenValue

Function EigenValue(r As Range) As Double()
    Call EigenValueDecomposition(r)
    EigenValue = x.Diag
End Function

固有ベクトルを並べた直交行列を出力する関数EigenVector

Function EigenVector(r As Range) As Double()
    Call EigenValueDecomposition(r)
    EigenVector = u.ToArray
End Function

メイン関数

Private Sub EigenValueDecomposition(r As Range)
    Set x = New SymmetryMatrix
    x.SetFromRange r

    Set u = New Matrix
    u.SetAsUnitMatrix x.Size
    
    Dim isConverge As Boolean, i As Long
    isConverge = False
    i = 0

    Do Until isConverge Or i >= MAX_ITERATION
        i = i + 1
        
        Dim g As GivensMatrix
        Set g = New GivensMatrix
        g.SetFromMaxOfUpperTriangularComponent x
        
        Set x = x.Rotate(g) 'x -> gxg
        Set u = u.Times(g)
        
        isConverge = (Abs(g.Theta) < CONVERGENCE_CRITERIA)
    Loop

    Call SortEigenvalue(x, u)
End Sub

ヤコビ法は固有値が大きい順にならないのでクイックソートで整頓する．

Private Sub SortEigenvalue(y As SymmetryMatrix, s As Matrix, _
                           Optional firstIndex As Long = 0, Optional lastIndex As Long = -1)
    'xとrを固有値の項順にソート
    
    '// ソート終了位置省略（最初の呼び出し）時は配列要素数を設定
    If lastIndex = -1 Then lastIndex = y.Size - 1
        
    '// 中央値を取得
    Dim medianIndex As Long, median As Double
    medianIndex = CLng((firstIndex + lastIndex) / 2)
    median = y.Value(medianIndex, medianIndex)

    Dim leftIndex As Long, rightIndex As Long
    leftIndex = firstIndex '// 左ループカウンタ
    rightIndex = lastIndex '// 右ループカウンタ

    Do
        '// 中央値の左側をループし，配列の左側から中央値より小さい値を探す
        Do Until y.Value(leftIndex, leftIndex) <= median
            leftIndex = leftIndex + 1
        Loop
        
        '// 中央値の右側をループし，配列の右側から中央値より大きい値を探す
        Do Until y.Value(rightIndex, rightIndex) >= median
             rightIndex = rightIndex - 1
        Loop
        
        If leftIndex >= rightIndex Then
            '// 左側の方が大きければここで処理終了
            Exit Do
        Else
            '// 右側の方が大きい場合は、
            Call 左右を入れ替える(y, s, leftIndex, rightIndex)
            
            leftIndex = leftIndex + 1 '// 左側を１つ右にずらす
            rightIndex = rightIndex - 1 '// 右側を１つ左にずらす
        End If
    Loop
    
    If firstIndex < leftIndex - 1 Then
        Call SortEigenvalue(y, s, firstIndex, leftIndex - 1)  '// 中央値の左側を再帰でクイックソート
    End If
    
    If rightIndex + 1 < lastIndex Then
        Call SortEigenvalue(y, s, rightIndex + 1, lastIndex) '// 中央値の右側を再帰でクイックソート
    End If
End Sub

Private Sub 左右を入れ替える(y As SymmetryMatrix, s As Matrix, leftIndex As Long, rightIndex As Long)
    Dim temp As Double
    
    'Call swap(y.Value(leftIndex, leftIndex), y.Value(rightIndex, rightIndex))
    temp = y.Value(leftIndex, leftIndex)
    y.Value(leftIndex, leftIndex) = y.Value(rightIndex, rightIndex)
    y.Value(rightIndex, rightIndex) = temp
    
    Dim i As Long
    For i = 0 To y.Size - 1
        'Call swap(s.Value(i, leftIndex), s.Value(i, rightIndex))
        temp = s.Value(i, leftIndex)
        s.Value(i, leftIndex) = s.Value(i, rightIndex)
        s.Value(i, rightIndex) = temp
    Next
End Sub

クラスモジュール Matrix

Private n_ As Long
Private r_() As Double

Sub SetAsUnitMatrix(n As Long)
    n_ = n

    Dim result() As Double
    ReDim result(n - 1, n - 1)

    Dim i As Long, j As Long
    For i = 0 To n - 1
    For j = 0 To n - 1
        result(i, j) = IIf(i = j, 1, 0)
    Next j, i

    r_ = result
End Sub

Property Get Size() As Long
    Size = n_
End Property

Property Let Size(n As Long)
    n_ = n
    ReDim r_(n - 1, n - 1)
End Property

Property Get Value(i As Long, j As Long) As Double
    Value = r_(i, j)
End Property

Property Let Value(i As Long, j As Long, v As Double)
    r_(i, j) = v
End Property

Function ToArray() As Double()
    Dim result() As Double
    ReDim result(n_ - 1, n_ - 1)
    
    Dim i As Long, j As Long
    
    For i = 0 To n_ - 1
    For j = 0 To n_ - 1
        result(i, j) = r_(i, j)
    Next j, i
    
    ToArray = result
End Function

Function ToSymmetryMatrix() As SymmetryMatrix
    Dim result As SymmetryMatrix
    Set result = New SymmetryMatrix
    result.Size = n_

    Dim i As Long, j As Long
    For i = 0 To n_ - 1
    For j = i To n_ - 1
        result.Value(i, j) = (r_(i, j) + r_(j, i)) / 2
    Next j, i

    Set ToSymmetryMatrix = result
End Function

Function Times(g As GivensMatrix) As Matrix
    'r -> r*g
    Dim rg As Matrix
    Set rg = Me.Clone

    Dim i As Long
    For i = 0 To n_ - 1
        rg.Value(i, g.Row) = r_(i, g.Row) * g.RowRow + r_(i, g.Column) * g.ColumnRow
        rg.Value(i, g.Column) = r_(i, g.Row) * g.RowColumn + r_(i, g.Column) * g.ColumnColumn
    Next i

    Set Times = rg
End Function

Function Clone() As Matrix
    Dim result As Matrix
    Set result = New Matrix
    result.Size = n_

    Dim i As Long, j As Long
    For i = 0 To n_ - 1
    For j = 0 To n_ - 1
        result.Value(i, j) = r_(i, j)
    Next j, i

    Set Clone = result
End Function

クラスモジュール SymmetryMatrix

仕様

入力が正方行列でない場合は行数・列数の小さいほうを大きさとする
x_(i,j)はi<=jのときのみ定義することで強制的に対称化する
入力が対称行列でない場合は平均をとることで強制的に対称化する

Private n_ As Long
Private x_() As Double

Sub SetFromRange(r As Range)
    n_ = IIf(r.Rows.Count < r.Columns.Count, r.Rows.Count, r.Columns.Count)

    Dim result() As Double
    ReDim result(n_ - 1, n_ - 1)

    Dim i As Long, j As Long
    For i = 0 To n_ - 1
    For j = i To n_ - 1
        result(i, j) = (CDbl(r(i + 1, j + 1)) + CDbl(r(j + 1, i + 1))) / 2
    Next j, i

    For i = 1 To n_ - 1
    For j = 0 To i - 1
        result(i, j) = 0
    Next j, i

    x_ = result
End Sub

Property Let Size(n As Long)
    n_ = n
    ReDim x_(n - 1, n - 1)
End Property

Property Get Size() As Long
    Size = n_
End Property

Property Let Value(i As Long, j As Long, v As Double)
    If i <= j Then x_(i, j) = v Else x_(j, i) = v
End Property

Property Get Value(i As Long, j As Long) As Double
    Value = IIf(i <= j, x_(i, j), x_(j, i))
End Property

Property Get Diag() As Double()
    '固有値を並べた行ベクトル
    Dim result() As Double
    ReDim result(n_ - 1)

    Dim i As Long
    For i = 0 To n_ - 1
        result(i) = x_(i, i)
    Next

    Diag = result
End Property

Function ToMatrix() As Matrix
    Dim result As Matrix
    Set result = New Matrix
    result.Size = n_

    Dim i As Long, j As Long
    
    For i = 0 To n_ - 1
    For j = 0 To n_ - 1
        result.Value(i, j) = Me.Value(i, j)
    Next j, i

    Set ToMatrix = result
End Function

Function Times(g As GivensMatrix) As Matrix
    'x -> xg
    Dim xg As Matrix
    Set xg = Me.ToMatrix

    Dim i As Long
    For i = 0 To n_ - 1
        xg.Value(i, g.Row) = Me.Value(i, g.Row) * g.RowRow + Me.Value(i, g.Column) * g.ColumnRow
        xg.Value(i, g.Column) = Me.Value(i, g.Row) * g.RowColumn + Me.Value(i, g.Column) * g.ColumnColumn
    Next

    Set Times = xg
End Function

Function Rotate(g As GivensMatrix) As SymmetryMatrix
    'x -> gxg
    Set Rotate = g.Transpose.Times(Me.Times(g)).ToSymmetryMatrix
End Function

クラスモジュール GivensMatrix

$\begin{bmatrix}1 \\ &1 \\ &&\cos\theta&&-\sin\theta \\ &&&1 \\ &&\sin\theta&&\cos\theta \\ &&&&&1 \\ &&&&&&1\end{bmatrix}$
みたいな形の行列．
あえて1本の式で書くと
　 $g_{i, j} = \delta_{i, j}+ (\cos\theta - 1) (\delta_{i, col} \delta_{j, col} + \delta_{i, row} \delta_{j, row})+ \sin\theta (\delta_{i, col} \delta_{j, row} - \delta_{i, row} \delta_{j, col}).$

Private row_ As Long
Private col_ As Long
Private theta_ As Double

Sub SetFromMaxOfUpperTriangularComponent(x As SymmetryMatrix)
    Dim maxValue As Double

    row_ = 0
    col_ = 1
    maxValue = x.Value(0, 1)

    Dim i As Long, j As Long

    For i = 0 To x.Size - 2
    For j = i + 1 To x.Size - 1
        If Abs(x.Value(i, j)) > Abs(maxValue) Then
            row_ = i
            col_ = j
            maxValue = x.Value(i, j)
        End If
    Next j, i

    If x.Value(row_, row_) - x.Value(col_, col_) <> 0 Then
        Dim tan2t As Double
        tan2t = 2 * maxValue / (x.Value(row_, row_) - x.Value(col_, col_))
        theta_ = Atn(tan2t) / 2 '確認：Atnの値域は-pi/2～pi/2
    Else
        theta_ = Atn(1) * Sgn(maxValue)
    End If
End Sub

Property Get Row() As Long
    Row = row_
End Property

Property Let Row(i As Long)
    row_ = i
End Property

Property Get Column() As Long
    Column = col_
End Property

Property Let Column(j As Long)
    col_ = j
End Property

Property Get Theta() As Double
    Theta = theta_
End Property

Property Let Theta(t As Double)
    theta_ = t
End Property

Property Get RowRow() As Double
    RowRow = cos(theta_)
End Property

Property Get ColumnColumn() As Double
    ColumnColumn = cos(theta_)
End Property

Property Get RowColumn() As Double
    RowColumn = -sin(theta_)
End Property

Property Get ColumnRow() As Double
    ColumnRow = sin(theta_)
End Property

Function Times(x As Matrix) As Matrix
    'gx
    Dim gx As Matrix
    Set gx = x.Clone

    Dim i As Long, j As Long

    For j = 0 To x.Size - 1
        gx.Value(row_, j) = Me.RowRow * x.Value(row_, j) + Me.RowColumn * x.Value(col_, j)
        gx.Value(col_, j) = Me.ColumnRow * x.Value(row_, j) + Me.ColumnColumn * x.Value(col_, j)
    Next j

    Set Times = gx
End Function

Function Clone() As GivensMatrix
    Dim result As GivensMatrix
    Set result = New GivensMatrix
    result.Row = row_
    result.Column = col_
    result.Theta = theta_
    Set Clone = result
End Function

Function Transpose() As GivensMatrix
    Dim result As GivensMatrix
    Set result = Me.Clone
    result.Theta = -theta_
    Set Transpose = result
End Function

wetchのブログ

他人に見られることを想定していない書き散らかし独習ノート．物理学とかVBAとか．