“數(shù)值分析”課程中稀疏矩陣?yán)碚摰慕虒W(xué)策略

鄭 華，羅 亮

（韶關(guān)學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院， 廣東 韶關(guān) 512005）

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

1 “數(shù)值分析”教學(xué)概況

“數(shù)值分析”課程分為數(shù)值代數(shù)、數(shù)值逼近和微分方程數(shù)值解3 個(gè)模塊， 其中數(shù)值代數(shù)部分主要關(guān)注實(shí)際數(shù)學(xué)模型離散化得到的線性方程組和特征值計(jì)算的數(shù)值算法. 很多大數(shù)據(jù)應(yīng)用問(wèn)題，比如微分方程的有限差分離散化［1］、網(wǎng)頁(yè)排序問(wèn)題［2-3］、復(fù)雜關(guān)系網(wǎng)絡(luò)問(wèn)題［4］等，往往會(huì)轉(zhuǎn)化為大規(guī)模稀疏矩陣的處理，在目前的一些常見(jiàn)的《數(shù)值分析》教材［5-6］中，對(duì)稀疏矩陣?yán)碚摪ㄅc之匹配的上機(jī)實(shí)驗(yàn)的討論較少，如果在教學(xué)過(guò)程中能設(shè)計(jì)有針對(duì)性的教學(xué)環(huán)節(jié)，就能更好地幫助學(xué)生理解和掌握相關(guān)的數(shù)值算法理論.

針對(duì)“數(shù)值分析”課程中與大規(guī)模稀疏矩陣相關(guān)的線性方程組求解和特征值計(jì)算兩部分內(nèi)容，本文給出了稀疏矩陣?yán)碚摰慕虒W(xué)策略，完成相關(guān)章節(jié)教學(xué)目標(biāo)的同時(shí)，使相關(guān)內(nèi)容的教學(xué)更具有針對(duì)性，以適應(yīng)大數(shù)據(jù)時(shí)代的需要.

2 稀疏矩陣教學(xué)策略

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

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

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

在線性方程組求解的實(shí)驗(yàn)教學(xué)過(guò)程中，如果采用MATLAB 進(jìn)行，學(xué)生可能會(huì)使用“”運(yùn)算符［7］.由于MATLAB 的“”運(yùn)算符所調(diào)用的mldivide 函數(shù)的內(nèi)置算法涉及到了SuperLU 軟件包［8］，對(duì)于大規(guī)模稀疏矩陣（尤其是具有一定結(jié)構(gòu)的稀疏矩陣）同樣適用，這容易讓學(xué)生對(duì)前期的理論講解感到困惑.

表1 展示了用“”運(yùn)算符求解不同類(lèi)別系數(shù)矩陣?yán)拥男Ч麑?duì)比.對(duì)于矩陣A（稠密），隨著矩陣階數(shù)的提升，計(jì)算時(shí)間和內(nèi)存明顯增加；對(duì)于不具有任何結(jié)構(gòu)的矩陣B（稀疏），雖然存成稀疏方式后可以明顯節(jié)省內(nèi)存，但由于矩陣階數(shù)更大，所耗費(fèi)的計(jì)算時(shí)間比例1 更多；對(duì)于矩陣C（三對(duì)角），使用“”運(yùn)算符的效果就得到了明顯的提升，可求解的矩陣階數(shù)可達(dá)到千萬(wàn)級(jí).

表1 “”運(yùn)算符的效果對(duì)比

結(jié)合這樣計(jì)算實(shí)例，可結(jié)合本章主要知識(shí)點(diǎn)給學(xué)生進(jìn)行說(shuō)明.和矩陣A 相比較，對(duì)矩陣B 的計(jì)算時(shí)間沒(méi)有減少，說(shuō)明“”運(yùn)算符本質(zhì)上是直接法，而對(duì)矩陣C 的計(jì)算效率明顯，說(shuō)明了“”運(yùn)算符的內(nèi)置代碼有其特殊性，進(jìn)而對(duì)SuperLU 軟件包做個(gè)簡(jiǎn)單介紹，強(qiáng)調(diào)這是在后續(xù)數(shù)值代數(shù)研究中會(huì)遇到的問(wèn)題，本科教學(xué)階段暫時(shí)不去涉及，一方面可以打消學(xué)生可能產(chǎn)生的疑惑，另一方面也為部分學(xué)生今后可能從事的相關(guān)研究培養(yǎng)興趣.

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

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

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

表2 MATLAB 稀疏矩陣常見(jiàn)函數(shù)教學(xué)要點(diǎn)

2.3 特征值計(jì)算

特征值計(jì)算包含兩個(gè)數(shù)學(xué)問(wèn)題：一是求大規(guī)模稀疏矩陣的某個(gè)特征值和對(duì)應(yīng)特征向量，二是求解小型稠密矩陣的全部特征值.

對(duì)于第一個(gè)數(shù)學(xué)問(wèn)題，涉及冪法和反冪法.根據(jù)大規(guī)模稀疏矩陣運(yùn)算的需要，冪法格式的構(gòu)造需要避免矩陣乘冪，以矩陣乘向量的迭代格式進(jìn)行，《數(shù)值分析》教材對(duì)此的講解較少，如果學(xué)生直接認(rèn)同教材最終給出的冪法迭代格式，就容易忽略?xún)绶ㄖ芯仃嚦藘绲谋举|(zhì).因此，在冪法格式構(gòu)造的課堂教學(xué)中，務(wù)必向?qū)W生演示直接做矩陣乘冪和以迭代格式進(jìn)行的區(qū)別，加深學(xué)生對(duì)“理論上等價(jià)，數(shù)值上未必等價(jià)”這個(gè)課程理念的理解.

在實(shí)驗(yàn)教學(xué)過(guò)程中，借助來(lái)源于實(shí)際應(yīng)用的超大規(guī)模稀疏矩陣，比如從University of Florida Sparse Matrix Collection （https：//sparse.tamu.edu/）中獲取的Web 矩陣，緊密結(jié)合實(shí)際應(yīng)用的實(shí)驗(yàn)教學(xué)任務(wù)，能更好地讓學(xué)生明確所學(xué)知識(shí)的應(yīng)用價(jià)值.對(duì)于反冪法的教學(xué)，由于該算法的本質(zhì)是對(duì)矩陣的逆應(yīng)用冪法，因此迭代格式的引入是簡(jiǎn)單并且自然的.到了反冪法的算法實(shí)現(xiàn)環(huán)節(jié)，明顯出現(xiàn)了矩陣求逆運(yùn)算，此處正好符合針對(duì)稀疏矩陣應(yīng)盡量避免矩陣求逆的準(zhǔn)則應(yīng)用，所以在反冪法的最終格式構(gòu)造上，也是稀疏矩陣?yán)碚撘氲囊粋€(gè)切入點(diǎn).

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

上述代碼對(duì)于小規(guī)模矩陣并沒(méi)有很明顯的缺陷，但一旦涉及大規(guī)模矩陣的情形，w*w’的操作顯然形成了大規(guī)模的稠密矩陣，會(huì)造成內(nèi)存溢出.因此，需要根據(jù)稀疏矩陣的運(yùn)算規(guī)則，修改代碼來(lái)避免存儲(chǔ)初等反射陣，即：

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

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

3 小結(jié)

基于“數(shù)值分析”的課程標(biāo)準(zhǔn)，在不影響課程本身教學(xué)目標(biāo)的基礎(chǔ)上，給出了稀疏矩陣?yán)碚撛陉P(guān)聯(lián)知識(shí)中的教學(xué)策略.所提出的教學(xué)策略可以總結(jié)為：以鞏固算法理論為基礎(chǔ)、以應(yīng)用為導(dǎo)向和以實(shí)驗(yàn)教學(xué)為輔助.該教學(xué)策略可以推廣到其他應(yīng)用型數(shù)學(xué)專(zhuān)業(yè)課的教學(xué)中.