Bullet_General.gif (3885 bytes) Solve the matrix eigenvalue problem by using Hotelling's deflation method. (線上執行) (原始碼下載)

一般的 matrix eigenvalue problem 大概以如下的型式出現

HX =λX

其中 H 是一個已知的方陣 (square matrix)，X 和λ分別是未知的 eigenvector  和 eigenvalue。如果 H 是一個 n 階方陣，則可以找出 n 組λ和 X 的 解，使其滿足 HX=λX。以下的程式碼使用 power method 配合 Hotelling's deflation 即可求出所有 n 組λ和 X 的解，為了避免遇到 complex eigenvalues 的問題，以下的程式碼只適用在 H 為對稱矩陣(symmetric matrix)的情況，若 H 為非對稱矩陣，則可能出現 complex eigenvalue，但是這種問題比較複雜，並不在本篇討論的範圍。關於 matrix eigenvalue problem 更進一步的資料，請參閱市面上的數值分析書籍。

　

副程式：

Public Sub HotellingDeflation(H() As Single, EigenValue() As Single, _
       EigenVector() As Single, Optional tol)
     Dim i As Long, j As Long, k As Long, row As Long, lambda(1 To 2) As Single, X() As Single
     Dim RHS() As Single, TX() As Single, RHSU() As Single, RHSD() As Single
     If IsMissing(tol) Then tol = 0.0001
     row = UBound(H, 1)
     ReDim X(1 To row, 1 To 1), EigenValue(1 To row), EigenVector(1 To row, 1 To row)
     For j = 1 To row: X(j, 1) = 1: Next
     For i = 1 To row
          Do
               lambda(1) = lambda(2)
               Call MatrixMultiply(H, X, RHS)
               lambda(2) = RHS(1, 1)
               For j = 1 To row: X(j, 1) = RHS(j, 1) / lambda(2): Next
          Loop Until Abs(lambda(2) - lambda(1)) < tol
          EigenValue(i) = lambda(2)
          For j = 1 To row: EigenVector(j, i) = X(j, 1): Next
          If i <> row Then
               Call MatrixTranspose(X, TX)
               Call MatrixMultiply(RHS, TX, RHSU)
               Call MatrixMultiply(TX, X, RHSD)
               For j = 1 To row
                    For k = 1 To row
                         H(j, k) = H(j, k) - RHSU(j, k) / RHSD(1, 1)
                    Next
               Next
          End If
     Next
End Sub

其中 MatrixTranspose 和 MatrixMultiply 這兩個副程式，請參閱本站的 矩陣的基本運算。

　

參數說明：

  • H：一個已知的對稱矩陣。
  • EigenValue：一維向量，傳回 H 的 n 個 eigenvalue。
  • EigenVector：二維向量，傳回 H 的 n 個 eigenvector，第一行(column)代表第一組 eigenvector，第二行代表第二組eigenvector，依此類推。
  • tol：以 power method 求 dominant eigenvalue 時的容忍誤差。這是個選擇性引數，如果忽略這個引數的話，tol 將自動被設成 0.0001。

　

範例：

1.由鍵盤輸入矩陣元素：

Private Sub Command1_Click()
     Dim i As Long, j As Long, Keyin As String, H() As Single
     Dim EigenValue() As Single, EigenVector() As Single
     Keyin = InputBox("請輸入對稱矩陣，以分號"";""來分隔列元素，以空白來分隔欄元素。" _
               & vbNewLine & "例如有一3×3的對稱矩陣:" & vbNewLine & "1 2 3 " & vbNewLine & _
               "2 3 4 " & vbNewLine & "3 4 5 " & vbNewLine & "則輸入: 1 2 3;2 3 4;3 4 5 ")
     If SepStrToMatrix(Keyin, ";", " ", H) Then
          If Not IsSymmetric(H) Then MsgBox "這不是對稱矩陣": Exit Sub
          Call HotellingDeflation(H, EigenValue, EigenVector)
          For i = 1 To UBound(H, 1)
               Debug.Print EigenValue(i)
               For j = 1 To UBound(H, 1)
                    Debug.Print EigenVector(j, i);
               Next
               Debug.Print
          Next
     Else
          MsgBox "矩陣輸入有誤，請重新輸入。"
     End If
End Sub

其中 SepStrToMatrix 這個副程式，請參閱本站的 將字串拆解成陣列。IsSymmetric這個副程式，請參閱本站的 矩陣的基本運算。

若有一對稱矩陣如下：

2 3 8
3 9 4
8 4 1

執行 Command1，並在 Inputbox 中輸入 2 3 8;3 9 4;8 4 1，在 VB 的即時運算視窗可得此對稱矩陣之 3 組 eigenvalue (奇數列) 和 eigenvector (偶數列)如下。

14.20678 
 1  1.357338  1.016852 
-6.568397 
 1  9.127734E-02 -1.105273 
 4.361548 
 1 -1.332108  0.794727

　

2.由檔案讀入矩陣元素：

Private Sub Command2_Click()
     Dim i As Long, j As Long, Keyin As String, H() As Single, FileNo As Long
     Dim EigenValue() As Single, EigenVector() As Single, temp As String
     FileNo = FreeFile
     Open "c:\eigen.txt" For Input As #FileNo
     Do While Not EOF(FileNo)
          Line Input #FileNo, temp
          If Trim(temp) <> "" Then Keyin = Keyin & temp & ";"
     Loop
     Close #FileNo
     Keyin = Mid(Keyin, 1, Len(Keyin) - 1)
     If SepStrToMatrix(Keyin, ";", " ", H) Then
          If Not IsSymmetric(H) Then MsgBox "這不是對稱矩陣": Exit Sub
          Call HotellingDeflation(H, EigenValue, EigenVector)
          For i = 1 To UBound(H, 1)
               Debug.Print EigenValue(i)
               For j = 1 To UBound(H, 1)
                    Debug.Print EigenVector(j, i);
               Next
               Debug.Print
          Next
     Else
          MsgBox "矩陣輸入有誤，請重新輸入。"
     End If
End Sub

其中 SepStrToMatrix 這個副程式，請參閱本站的 將字串拆解成陣列。IsSymmetric這個副程式，請參閱本站的 矩陣的基本運算。

若有一矩陣儲存在"C:\eigen.txt"如下所示

11 2 3 1 4 2
2 9 3 5 2 1
3 3 15 4 3 2
1 5 4 12 4 3
4 2 3 4 17 5
2 1 2 3 5 8

執行 Command2 之後，在 VB 的即時運算視窗可得此對稱矩陣之 6 組 eigenvalue (奇數列) 和 eigenvector (偶數列) 如下。

27.99144 
 1  0.9414517  1.448762  1.367093  1.866725  0.9641068 
 13.86531 
 1 -2.472871 -4.860624 -2.366373  5.374027  1.632028 
 10.9771 
 1 -0.581289  0.9389649 -1.139224 -6.212605E-02 -0.1448846 
 8.712057 
 1  0.6025767 -0.6092691  0.1080833 -0.3641766 -0.1582333 
 5.988571 
 1 -2.508974 -0.3602135  1.353895 -2.622937  5.11286 
 4.463624 
 1 -2.618987 -0.1057338  2.37532  0.1748679 -2.027661


以 power method 配合 Hotelling's deflation 並不是求 eigenvalue 最有效率的方法，只是程式碼容易寫而已，如果遇到比較大的矩陣，可能必須考慮使用 Householder's Method 配合 LR 或是 QR algorithm 才會得到較快的執行速度，這部份以後有機會再介紹。另外，HotellingDeflation 副程式中的 power method 並沒有使用任何 acceleration technique 的技巧，如果覺得收斂速度太慢，可以自行加入 relaxation factor 或 Rayleigh Quotient 等等加快 power method 收斂速度的方法之程式碼。

3/31/1999 補充：

有某種情況會使得 HotellingDeflation 副程式無法收斂，導致形成無窮迴圈。請注意本專案中的浮點變數都是宣告成 Single，因此當 tol 容忍誤差限制的太小時，就會出現無法收斂的情形，解決方法是將 tol 容忍誤差設的大一點；如果一定要很小的容忍誤差的話，則請將整個專案的 Single 變數改宣告成 Double 也可以解決這類問題，請參考以下的範例：

使用相同的範例，若有一對稱矩陣如下：

2 3 8
3 9 4
8 4 1

欲使用的 tol 容忍誤差是 0.0000001 (如紅色部分所示)，則程式碼如下：

Private Sub Command3_Click()
     Dim i As Long, j As Long, Keyin As String, H() As Single
     Dim EigenValue() As Single, EigenVector() As Single
     Keyin = InputBox("請輸入對稱矩陣，以分號"";""來分隔列元素，以空白來分隔欄元素。" _
               & vbNewLine & "例如有一3×3的對稱矩陣:" & vbNewLine & "1 2 3 " & vbNewLine & _
               "2 3 4 " & vbNewLine & "3 4 5 " & vbNewLine & "則輸入: 1 2 3;2 3 4;3 4 5 ")
     If SepStrToMatrix(Keyin, ";", " ", H) Then
          If Not IsSymmetric(H) Then MsgBox "這不是對稱矩陣": Exit Sub
          Call HotellingDeflation(H, EigenValue, EigenVector, 0.0000001)
          For i = 1 To UBound(H, 1)
               Debug.Print EigenValue(i)
               For j = 1 To UBound(H, 1)
                    Debug.Print EigenVector(j, i);
               Next
               Debug.Print
          Next
     Else
          MsgBox "矩陣輸入有誤，請重新輸入。"
     End If
End Sub

執行 Command3，並在 Inputbox 中輸入 2 3 8;3 9 4;8 4 1，即發現程式進入無窮迴圈，停不下來了，此時把容忍誤差改成大於 0.000001 即可；或是把整個專案的 Single 變數改宣告成 Double 也可解決這類問題，請參考以下的範例：

Private Sub Command4_Click()
     Dim i As Long, j As Long, Keyin As String, H() As Double
     Dim EigenValue() As Double, EigenVector() As Double
     Keyin = InputBox("請輸入對稱矩陣，以分號"";""來分隔列元素，以空白來分隔欄元素。" _
               & vbNewLine & "例如有一3×3的對稱矩陣:" & vbNewLine & "1 2 3 " & vbNewLine & _
               "2 3 4 " & vbNewLine & "3 4 5 " & vbNewLine & "則輸入: 1 2 3;2 3 4;3 4 5 ")
     If SepStrToMatrix(Keyin, ";", " ", H) Then
          If Not IsSymmetric(H) Then MsgBox "這不是對稱矩陣": Exit Sub
          Call HotellingDeflation(H, EigenValue, EigenVector, 0.0000001)
          For i = 1 To UBound(H, 1)
               Debug.Print EigenValue(i)
               For j = 1 To UBound(H, 1)
                    Debug.Print EigenVector(j, i);
               Next
               Debug.Print
          Next
     Else
          MsgBox "矩陣輸入有誤，請重新輸入。"
     End If
End Sub

執行 Command4，並在 Inputbox 中輸入 2 3 8;3 9 4;8 4 1，在 VB 的即時運算視窗可得此對稱矩陣之 3 組 eigenvalue (奇數列) 和 eigenvector (偶數列)如下。

14.2068125885759 
 1  1.3573372140458  1.01685011285943 
-6.56836129734352 
 1  9.12814641801065E-02 -1.10527570326121 
 4.3615487104883 
 1 -1.33211382052807  0.794736266149012

請注意是整個專案的 Single 變數都要改宣告成 Double，而不只是HotellingDeflation 副程式，不然會出現型態不符合的錯誤訊息

類似這樣的問題並不只是會發生在 HotellingDeflation 副程式中，只要是程式碼中有牽涉到是否收斂的問題(例如使用牛頓法求根)，都要注意容忍誤差是否小於變數本身的精確度，如果是的話，程式可能會進入無窮迴圈。

[ 上一個 | 首頁 | 數值分析 | 下一個 ]

This page was written by Jaric on Feb. 26, 1999. All rights reserved.

Revised : Mar. 31 , 1999 .

Total pageview since 4/6/1999.