| 葉明養 |
|
有關特徵值的計算，我們常常經由數學或工程方面的教科書，得知可使用行列式 | Α － λΙ| = 0 的方式求得方程式的根，其缺點是這種方法不易求得特徵向量，而且方程式的幕次僅能至4次方，超過5次方時並無公式可使用，這一點已由挪威數學家 Niels Henrik Abel證明了，一般的計算都是使用英國計算機科學家John G. F. Francis所提出的QR演算法，QR演算法也是二十世紀十大演算法之一，此法能幫助我們改善求解矩陣的特徵值和特徵向量，但因矩陣元素之間的關係和結構特性複雜度，尤其有些特殊矩陣至今仍然存在諸多尚待解決的問題。所謂QR演算法就是將已知矩陣QR分解，包含三種分解法，Givens Rotator(旋轉器)、Householder Reflector(翻轉器)和Gram-Schmidt Process(過程)，分解之後再使用左右迭代運算(Left and Right Iterative Computation)。筆者已將上述各種類別分別使用CSharp程式語言實作並在微軟Visual Studio開發環境中作整體性的測試，其結果尚稱滿意。 |
|
一、前言 以一個營造廠專任工程人員兼具土木結構工程師的角色，應用資料結構學(Data Structure)的知識、SQL(Structural Query Language)和CSharp程式語言撰寫經驗，再結合線性代數(Linear Algebra)與矩陣計算(Matrix Computation)等計算數學(Computational Mathematics)方面的領域知識(Domain Knowledge)，撰寫本篇文章，無論是土木結構工程、資訊和電機 . . . 等等各方面的理論和應用都是相通的，雖然自己並非專業，仍然去學習這些新的知識，感嘆人生短暫，時間在不知不覺中悄悄地過去了。演算法(Algorithms)是電資學院重要的課程，並與資料結構密切結合，單就經典的教科書”Introduction To Algorithms”乙書就厚達五公分，要了解都不知要花多少時間？演算法本身雖然與各種程式語言無關，但使用程式語言是與電腦硬體和作業系統相溝通的媒介，掌握熟練先進的程式語言，相關的編譯器和軟體開發工具，這些都可由網路上隨手下載，並提供豐富的軟體程式開發資源，令人享用不盡。 矩陣是一種解決方案(Solution)的工具，應用在機率與統計、微分方程、數學分析、以及種種工程方面的問題，譬如結構動力和靜力分析、系統識別、控制系統 . . . 等等，把矩陣看成柱向量(Column Vector)有次序的集合，由第0個開始至n-1個為止，每個柱向量在座標軸上的投影數據，由第0個開始至m-1個為止，那矩陣就以FmXn表示，F是英文的Field，一般數學的書翻譯成“體”、“域”、“場”，在資料庫系統則稱為欄位上的值，筆者稱為“Frame”或是座標軸上的數據。F可能複數(Complex)也可能是實數(Real)，為了配合C#程式語言，本文僅討論實數矩陣Rmxn的輸入，但矩陣計算後的輸出結果，可以是實數也可以是複數。 一般線性代數的教科書，所討論的主要部分包含矩陣的基礎理論和矩陣的秩等等，而矩陣計算專書討論的內容，包含矩陣的分解、直接解法、迭代解法、求固有值、奇異值等等 。我們常常在線性代數的書看到A = Q R和A Q = Q D這樣的表示式，以一個程式設計者而言，或許應該寫成A = Q * R和A * Q = Q * D較合適，A是已知的矩陣，Q是待求解的正交( Orthonormal ) 矩陣，若Q是複數元素時稱為Unitary Matrix，但筆者撰寫的程式中Q為實數，很多數學家或是物理學家稱Q為量綱，筆者想Q或許稱量綱矩陣或是座標軸矩陣更合適，另外尚待求解的R是上三角矩陣( Upper Triangle Matrix )及D是對角線矩陣( Diagonal Matrix )，事實上D應稱為Schur矩陣更合適，矩陣相乘的運算子( Operator * )與一般實數相乘的符號類似，C#程式語言稱為Operator Overloading，如何將數學符號和C#程式語言結合在一起，以及如何撰寫程式碼也是非常重要，而且環環相扣。 二、矩陣之QR分解 將已知A矩陣分解成座標軸Q( Orthonormal )矩陣，和上三角R矩陣，也就是所謂QR分解，勿將QR分解與QR演算法混為一談，QR分解可一步達成，而QR演算法是以QR分解為基本，重複迭代運算至矩陣成為對角線矩陣，更精確講應為Schur矩陣為止，但如何達到Schur矩陣也是另外一回事。故必須撰寫程式碼檢查，在運算過程中若是Schur矩陣時就停止，由Schur矩陣再去求得固有值和固有向量，但並非一定是Schur矩陣，有時甚至是對角線矩陣，如何拿捏？可能連書本都不會提到這些事情，必須自己撰寫程式碼或許才會體會到。故在此特別強調自己撰寫程式碼的重要性，而不是完全一味地使用商業應用軟體。線性代數的書中 ，一般都是使用Gram-Schmidt Process的方法闡釋QR分解，基本上是由已知矩陣的柱向量為基本單位，產生彼此相互垂直的向量，一般稱為Classical Gram-Schmidt Process，另外很多矩陣計算的書中提供另一種改良版，叫作Modified Gram-Schmidt Process，若已知A是奇異矩陣，雖然產生的Q矩陣不會是奇異矩陣，但不是真正的座標軸矩陣，故兩者並無多大的差別。筆者以C#程式語言實作，分別稱為GS2和GS類別，這種現象也是自己撰寫程式碼後才發現的，兩者都已包含在整體的程式類別庫中。 另外兩種QR分解，分別為美國數學及計算機科學家 Alston Scott Householder所建構的 Reflection(翻轉計算法)，另一種為美國數學家Wallace Givens 所建構的 Rotation(旋轉計算法)，基本上此兩種計算方法都是以矩陣的對角線元素為Pivot點，依向量的順序由左上往右下來處理，在矩陣計算教科書也都有談到，筆者以C#程式語言實作分別稱為Reflector(翻轉器)和Rotator(旋轉器)，兩種類別。 上述三種QR分解的優點和缺點為何？對於End User而言並不在乎，主要能快速跑出結果來即可，那是程式開發者要去處理的事情。甚至End User不懂C#程式語言也無所謂，只要輸入類別名稱和已知的引數矩陣物件(Argument Matrix Object)，就可得到所要的矩陣計算結果，就像一般人使用手機APP，或是使用數學計算軟體Matlab一樣容易，在微軟 .Net的環境下，矩陣本身就是是物件，在C#程式碼中可以任意使用和傳送，換言之就是矩陣物件的傳送(Passed by Matrix Object)，對於End User非常方便，就以筆者個人而言，筆者僅會使用幾個Matlab的函數作必要的驗證，對於Matlab的語法筆者不懂，筆者對C#程式語言較有興趣，所要必須自行開發線性代數和矩陣計算的C#類別庫，使用這樣的類別庫就像使用Matlab商業數學軟體一樣容易。 三、QR左右迭代運算 已知矩陣A經由QR分解後，即以A = Q * R表示 ，Q是可逆的座標軸矩陣，若以i的superscript表示逆(inverse)矩陣，則R = Qi * A，再將方程式等號兩邊右側各乘Q矩陣，即 R * Q = Qi * A * Q，如此恰好是一個循環，經由k次循環後，即Rk * Qk = Qki * Ak * Qk 。再設定 Rk * Qk = Ak + 1 ，A = A0 則 : Ak+1 = Qki * Qk–1i * . . . * Q0i * A0 * Q0 * Q1 * Q2 . . . * Qk-1 * Qk . . . ( 1 ) 設定 : Ak+1 = D， Qi = Qki * Qk–1i * . . . * Q0i ，Q = Q0 * Q1 * Q2 . . . * Qk-1 * Qk 結果如下 : D = Qi * A * Q 或 A * Q = Q * D . . . . . . . . . . . . . . . . . ( 2 ) 公式( 2 )是QR演算法導出的結果，A是已知正方形矩陣，Q是的座標軸矩軸(Orthonormal Matrix)，D是對角線矩陣(嚴格說是Schur矩陣，不一定是對角線矩陣)，雖然公式非常簡單，但是每一步都需要計算機重複計算，重複檢查運算中的矩陣是否為Schur矩陣？不可能再使用手算，而且已知A矩陣元素(Entries)間的結構關係，使得對角線元素為零，就無法成為Schur矩陣，計算機計算的結果失敗是常見的事情，譬如中空(Hollow)矩陣，Toeplitz矩陣等等，往往不見得都能找到固有值和固有向量。 依據筆者個人測試的結論，已知A矩陣是對稱時，很容易求得特徵值，亦即對角線矩陣D，若A矩陣是非對稱矩陣時，特徵值可能是實數(Real Number)也可能是複數(Complex Number) ，甚至求不出特徵值。矩陣計算的教科書提到改善的方法，譬如實數矩陣使用Single-Shift的策略，若是複數矩陣使用Double-Shift的策略，此外尚有使用Precondition的方法來改善Ill-Posed或是Ill-Conditioned的矩陣，矩陣計算理論的闡述和矩陣計算的專書都是使用Matlab數學軟體。Matlab是評價極高的數學軟體，但筆者僅會使用Matlab函數，譬如eig() . . . 等等，筆者不會Matlab的語法，目前也不會為了學線性代數與矩陣計算而勉強自己學Matlab的語法，CSharp才是筆者喜歡的程式語言，故必須自己再撰寫矩陣計算的類別庫。筆者曾經讀到中國大陸有人使用C#與Matlab數學軟件混合編程，是以Matlab的m檔案在 Visual Studio的環境下執行，這是個混合體的應用，筆者沒有試過，可能會遭遇諸多問題，而且可能效率也不是很好？正統的方法應該是C#程式語言在 .Net Visual Studio的環境下執行。 四、撰寫CSharp程式類別庫 所謂類別庫(Class Library)，並非僅單一類別，而是由多個類別組織而成，就以QR演算法所撰寫的這個類別而言，不是只有這個類別就可以一步到位，尚需多個輔助類別，譬如需要利用兩個和三個矩陣相乘的類別、矩陣的行列式(Determinant)類別、矩陣數據顯示在螢幕上的類別、矩陣是否已經是對角線矩陣的類別 . . . 等等，所以需要開發很多個輔助類別(Help Classes)， 並協助共同運算，C#程式開發者，必需自行撰寫不同的類別程式，所以稱為類別庫，這就是所謂軟體IC的概念，無需撰寫龐大單一複雜的程式碼，而是小型獨立的程式碼，無論是維護或是修改都比較容易，現代程式碼的撰寫都是這種趨勢，這也是物件導向程式(Object-Oriented Programming Language)的優點之一。在此也特別提到微軟共同語言(Common Language)的開發環境，包含.Net Framework豐富的軟體資源和Visual Studio(VS)應用程式開發環境，VS是蠻先進的編譯器，譬如C#語法有錯誤，會立即顯示在螢幕上，Intellisense(智慧型顯示)對於程式的撰寫有很大的幫助，譬如顯示類別的元素(Property，Method，. . .)以及相對應的參數，物件瀏覽器(Object Brawser) . . . 等等，對於程式開發者是一大福音，筆者撰寫的類別庫，目前約有一百三十個類別，總共計一萬二千行左右，VS甚至把類別按字母排列，能快速找到所要的類別位置，非常方便。其他如Just-In-Time編譯器，能快速將C#程式碼轉譯成中間(Intermedia)語言，再轉為Native機械碼，動態連接程式庫(Dynamic Link Library)，在Dynamic Time(Runtime)而不是Compiling Time時再把需要的類別連進來，其他如垃圾收集器 . . . 等等，先進的開發工具(Develop Tools)對於程式設計師撰寫程式碼有很大的幫助。 談到筆者個人撰寫的三種不同QR分解法，即Reflector、Rotator和GS，接著需要再使用左右迭代運算(Left-Right Iteration Computation)，並配合Shift和Preconditioned的技術，建構LRReflector、LRRotator、和LRGS三個QR演算法類別，計算特徵值和特徵向量(筆者自稱之為座標軸矩陣或量綱矩陣)。另奇異值的分解(Singular Value Decompositon)與特徵值的分解可說是兩個軟體的孿生兄弟，為了配合SVD再建構LRReflector2和LRGS2兩個類別，共計五個QR演算法類別，筆者個人自行測試的結果尚稱滿意。 在本文最後段列出三個實例測試，並將執行的結果顯示出來，即使不懂CSharp程式碼，仍然可以感覺，看程式碼就像看文稿(manual script)一樣容易。物件導向程式語言的優點之一，就是把較複雜的程式碼隱藏在類別中，而使用者只要引用類別即可，就像使用Matlab數學軟體的函數一樣容易。 五、結論 理論的驗證往往需要借由撰寫程式碼來實現，軟體模擬(Software Simulation)結構分析和設計 . . . 等等，就以結構分析而言，一般都是使用商業結構應用軟體，譬如SAP2000、ETAB、MIDAS、STAAD PRO . . . 等等，稱為結構專業(Structural Domain Specified Application)軟體，而其他專業領域方面的軟體，譬如流體力學、訊號系統、系統識別(System Identification) . . .等等，常常使用MATLAB數學軟體，MATLAB提供很多矩陣計算方面的函數，譬如矩陣相加、相乘、求逆矩陣、求特徵值和奇異值 . . . 等等，這些軟體統稱為專業領域程式語言(Domain Specificed Programming Language)，而一般目的的程式語言(General Purpose Programming Language)，包含Fortran、C、C++、Java、C# 、Python，R、. . . 等等，筆者嘗試以C#程式語言和線性代數與矩陣計算的知識，撰寫類似MATLAB矩陣計算的應用軟體，筆者個人能力平庸、時間不足，僅作初步的測試，雖然結果尚稱滿意，希望將來能更廣泛的測試與驗證。筆者深深認為，若要完全了解專業知識之後，再開始撰寫程式軟體，那不知要等多久，即使了解三四成的專業知識，都可以開始撰寫程式碼，藉由撰寫程式碼，可以進一步輔助並了解更多的專業知識，兩者完全相輔相乘。 記得多年前參加一個土木結構工程研討會，有位資深的教授談到，以前土木系比電機系有更多的人，使用Fortran撰寫應用程式，時至今天使用個人電腦比過去的超級電腦，更方便、更便宜、而且跑的更快，不同的程式語言及開發介面提供選擇，也借此機會拋磚引玉，無論工程界和學術界在這方面大家共同努力。 參考書目(Bibliography) 1.林聰悟、林佳慧，”數值方法與程式”，圖文技術服務有限公司2.葉明養，”求解矩陣之固有值與奇異值-使用C#程式語言實作“，現代營建 2016年12月 444期 3.周志成，”線代啟示錄“，https://ccjou.wordpress.com 4.Howard Anton, “Elementary Linear Algebra 7 Edition”，John Wiley & Soncs, Inc. 5.Howard Anton & Cris Rorres, “Elementary Linear Algebra with Supplemental Applications 10 Edition”，Wiley Plus 6.林清凉，“從物理學切入的線性代數導論”，五南出版社 7.Stephen H. Friedberg，Arnold J. Insel，Lawrence E. Spence，“Linear Algebra”，Pearson Education International 8.Biswa Nath Datta，“Numbeical Linear Algebra and Applications Second Edition“， Society for Industrial and Applied Mathematics(SIAM) 9.James W. Demmel，”Applied Numerical Linear Algebra“，SIAM 10.David S. Watkins，”Fundamentals of Matrix Computations Third Edition“，A John Wiley & Sons, Inc., Publication 11.Gene H. Golub，Charles F. Van Loan，”Matrix Computations 4th Edition“，The Johns Hopkins University Press Baltimore 12.徐樹方，錢江，”矩陣計算六講“，高等教育出版社 13.J. H. Wilkinson，”The Algebraic Eigenvalue Problem”，Oxford Science Publications 14.Nicholas J. Higham，”Accuracy and Stability of Numerical Algorithms“，SIAM 15.Gilbert Strang，”Linear Algebra and Its Applications Fourth Edition“，Brooks/Cole Cengage Learnin maintainability and safety (RAMS)’, 1999 |
