“數值分析”課程中稀疏矩陣理論的教學策略

鄭 華，羅 亮

（韶關學院 數學與統計學院， 廣東 韶關 512005）

隨著信息技術的不斷革新和發展，大數據技術已經逐步滲透到實際應用的各個行業和業務職能領域，逐漸成為重要的生產要素，如何更好地運用海量數據，決定了新一輪生產率增長和消費者激增浪潮的到來.大數據背景下的信息技術需要良好的科學計算能力作為基礎，這對高校相關專業的人才培養提出了新的要求，特別是對與之關聯緊密的信息與計算科學專業核心課程進行改革創新.“數值分析”是信息與計算科學專業最核心的專業必修課，該課程關注怎樣用計算機求解各類數學模型轉換得到的數學問題， 并作相關的算法分析.

1 “數值分析”教學概況

“數值分析”課程分為數值代數、數值逼近和微分方程數值解3 個模塊， 其中數值代數部分主要關注實際數學模型離散化得到的線性方程組和特征值計算的數值算法. 很多大數據應用問題，比如微分方程的有限差分離散化［1］、網頁排序問題［2-3］、復雜關系網絡問題［4］等，往往會轉化為大規模稀疏矩陣的處理，在目前的一些常見的《數值分析》教材［5-6］中，對稀疏矩陣理論包括與之匹配的上機實驗的討論較少，如果在教學過程中能設計有針對性的教學環節，就能更好地幫助學生理解和掌握相關的數值算法理論.

針對“數值分析”課程中與大規模稀疏矩陣相關的線性方程組求解和特征值計算兩部分內容，本文給出了稀疏矩陣理論的教學策略，完成相關章節教學目標的同時，使相關內容的教學更具有針對性，以適應大數據時代的需要.

2 稀疏矩陣教學策略

如果一個矩陣的非零元較少，并且排列沒有規律，則稱它為稀疏矩陣.在各類數值算法的計算機實現中，對于稀疏矩陣的存儲，尤其是矩陣階數比較大的情形，如果采用常規矩陣的存儲方式，那將會占用大量的計算機內存，導致上機實驗的結果出現不準確，甚至計算機的運行出現內存溢出.因此，對待稀疏矩陣一般采用的是三元數組存儲的方法，只存儲非零元及其位置信息.另一方面，為了使矩陣在運算過程中保持其稀疏結構不被破壞，還要盡量避免矩陣乘矩陣和矩陣求逆的運算.這些算法實現的基本準則，在“數值分析”的算法理論教學中一般較少涉及.接下來按數值代數中線性方程組求解和特征值計算兩部分內容，分別給出針對性的稀疏矩陣教學策略.

2.1 求解線性方程組的直接法

對于求解線性方程組直接法部分的教學內容，是基于Gauss 消去法展開的，包括考慮減少小主元影響所得的主元素Gauss 消去法、對稱正定矩陣情形的平方根法和三對角矩陣情形的追趕法.由于直接法的本質是基于矩陣初等變換進行的，理論上等價于相應的矩陣分解（如LU 分解、Cholesky 分解等），因此，各類直接法在計算機實現中不可避免地會出現矩陣乘矩陣的運算，對于本科階段的數值分析學習要求來說，直接法破壞了矩陣的稀疏結構，是不適用于大規模稀疏矩陣的.因此，在教學過程中，先按照“數值分析”課程教學要求完成常規的教學內容，上機實驗所涉及的矩陣都采用小規模矩陣進行，讓學生先掌握各類直接法的算法思想和算法流程.在階段總結過后，再引入稀疏矩陣的概念，并以隨機生成的稀疏矩陣為例，向學生展示直接法應用過程中遇到的內存溢出、計算時間過長等麻煩，展示直接法的缺點，為后續求解線性方程組的迭代法的引入做好鋪墊.

在線性方程組求解的實驗教學過程中，如果采用MATLAB 進行，學生可能會使用“”運算符［7］.由于MATLAB 的“”運算符所調用的mldivide 函數的內置算法涉及到了SuperLU 軟件包［8］，對于大規模稀疏矩陣（尤其是具有一定結構的稀疏矩陣）同樣適用，這容易讓學生對前期的理論講解感到困惑.

表1 展示了用“”運算符求解不同類別系數矩陣例子的效果對比.對于矩陣A（稠密），隨著矩陣階數的提升，計算時間和內存明顯增加；對于不具有任何結構的矩陣B（稀疏），雖然存成稀疏方式后可以明顯節省內存，但由于矩陣階數更大，所耗費的計算時間比例1 更多；對于矩陣C（三對角），使用“”運算符的效果就得到了明顯的提升，可求解的矩陣階數可達到千萬級.

表1 “”運算符的效果對比

結合這樣計算實例，可結合本章主要知識點給學生進行說明.和矩陣A 相比較，對矩陣B 的計算時間沒有減少，說明“”運算符本質上是直接法，而對矩陣C 的計算效率明顯，說明了“”運算符的內置代碼有其特殊性，進而對SuperLU 軟件包做個簡單介紹，強調這是在后續數值代數研究中會遇到的問題，本科教學階段暫時不去涉及，一方面可以打消學生可能產生的疑惑，另一方面也為部分學生今后可能從事的相關研究培養興趣.

2.2 求解線性方程組的迭代法

求解線性方程組的迭代法，是專門針對大規模稀疏矩陣設計的，這部分算法迭代格式的構造主要包括兩類方式：一類是基于矩陣分裂進行的，包括經典的Jacobi 迭代、Gauss-Seidel 迭代和SOR 迭代；另一類是基于子空間迭代理論進行的，比如共軛梯度法.這類算法的特點是每一步迭代只涉及矩陣乘向量的運算，不會破壞矩陣的稀疏結構，適用于大規模稀疏矩陣.對于該部分內容的教學，在講解完各算法的基本思想和流程后，上機實驗的實例都充分利用大型稀疏矩陣，具體的實例生成可以采用隨機生成和特殊生成的方式進行.算法的實現過程不但關注各類迭代法的實現，同時注意跟直接法進行對比，展示算法運行的計算時間，使學生更好地理解迭代法的優勢所在.

在實驗教學中，由于MATLAB 內置了大量的稀疏矩陣操作函數，這為迭代法的算法實現提供了極大便利.為了便于實驗教學的展開，在迭代法實現的上機實驗之前，需設置MATLAB 稀疏矩陣相關函數的調用任務，主要包括sparse、full、nzz、spy、sprand、spdiags 等.同時，考慮到教學課時的限制，這些函數的學習和使用只要求學生掌握最基本的調用方式，如表2 所示.

表2 MATLAB 稀疏矩陣常見函數教學要點

2.3 特征值計算

特征值計算包含兩個數學問題：一是求大規模稀疏矩陣的某個特征值和對應特征向量，二是求解小型稠密矩陣的全部特征值.

對于第一個數學問題，涉及冪法和反冪法.根據大規模稀疏矩陣運算的需要，冪法格式的構造需要避免矩陣乘冪，以矩陣乘向量的迭代格式進行，《數值分析》教材對此的講解較少，如果學生直接認同教材最終給出的冪法迭代格式，就容易忽略冪法中矩陣乘冪的本質.因此，在冪法格式構造的課堂教學中，務必向學生演示直接做矩陣乘冪和以迭代格式進行的區別，加深學生對“理論上等價，數值上未必等價”這個課程理念的理解.

在實驗教學過程中，借助來源于實際應用的超大規模稀疏矩陣，比如從University of Florida Sparse Matrix Collection （https：//sparse.tamu.edu/）中獲取的Web 矩陣，緊密結合實際應用的實驗教學任務，能更好地讓學生明確所學知識的應用價值.對于反冪法的教學，由于該算法的本質是對矩陣的逆應用冪法，因此迭代格式的引入是簡單并且自然的.到了反冪法的算法實現環節，明顯出現了矩陣求逆運算，此處正好符合針對稀疏矩陣應盡量避免矩陣求逆的準則應用，所以在反冪法的最終格式構造上，也是稀疏矩陣理論引入的一個切入點.

對于第二個數學問題，涉及的主要算法是QR 方法，算法細節包括利用Householder 變換對矩陣進行QR 分解以及把矩陣通過正交相似變換化為Hessenberg 矩陣.由于教材上對QR 方法的介紹明確指出該方法只適用于小型稠密矩陣，容易讓學生認為上述相關理論都不適用于大規模稀疏矩陣.由于細節操作中的Householder 變換是可以通過數學公式變形避免矩陣乘矩陣的運算，因此在QR 方法準備知識的講解中，要給學生做大規模稀疏矩陣的課堂實驗演示，讓學生掌握這個經典正交變換算法的優勢.例如，設H ＝I-2wwT是初等反射陣，在進行Householder 變換時，算法過程的基本操作是矩陣向量乘，即Hx.這種簡單的表達式很容易讓學生在進行算法實現時把MATLAB 代碼寫成：

上述代碼對于小規模矩陣并沒有很明顯的缺陷，但一旦涉及大規模矩陣的情形，w*w’的操作顯然形成了大規模的稠密矩陣，會造成內存溢出.因此，需要根據稀疏矩陣的運算規則，修改代碼來避免存儲初等反射陣，即：

顯然，新的代碼避免了大規模稠密矩陣的生成.

另一方面，對于最終的QR 方法迭代格式，由于其中的步驟涉及了把矩陣進行QR 分解后再反轉相乘的操作，導致了該算法不適用于大規模稀疏矩陣，這是需要對學生強調的要點.除了課程知識的講解，還可以結合實際應用給學生做簡單的介紹，對于大規模稀疏矩陣，一般在應用中也不需要去求解全部特征值，從應用的角度幫助學生肯定QR 方法的價值.

3 小結

基于“數值分析”的課程標準，在不影響課程本身教學目標的基礎上，給出了稀疏矩陣理論在關聯知識中的教學策略.所提出的教學策略可以總結為：以鞏固算法理論為基礎、以應用為導向和以實驗教學為輔助.該教學策略可以推廣到其他應用型數學專業課的教學中.