固體結構大規模矩陣正定性判定的快速算法1)

徐 輝 劉 彬

(清華大學航天航空學院工程力學系應用力學教育部重點實驗室,北京 100084)

固體結構大規模矩陣正定性判定的快速算法1)

徐 輝 劉 彬2)

(清華大學航天航空學院工程力學系應用力學教育部重點實驗室,北京 100084)

對于結構穩定性分析中超大規模矩陣正定性判定,必須采用并行計算方法,傳統方法如計算特征值、主子式行列式及LDLT等直接方法難以實現.本文提出了一些可適用于并行的迭代判定算法.借鑒力學系統中能量下降的思想,發展了一種判定矩陣正定性的新思路,即將矩陣的正定性判定問題轉化為一個優化問題,并基于優化算法來判定矩陣的正定性.提出了基于最速下降法和共軛梯度法來進行矩陣正定性判定的算法.然后考慮到力學系統剛度矩陣的稀疏性和結構剛失穩狀態的弱非正定性,提出可以先截超平面后解方程求駐值點的方法來判定弱非正定矩陣的正定性.為了保證對強非正定矩陣判定的準確性,本文提出可以高效混雜使用截平面法和共軛梯度法.數值實驗結果表明,本文提出的算法具有準確性和高效性.對于強非正定矩陣而言,共軛梯度算法更加高效;而對于弱非正定矩陣,則是截平面法和混雜算法更加高效.這些算法都容易在機群上實現并行計算,能夠快速判定大規模矩陣的正定性.

矩陣,正定性,共軛梯度法,截平面

引言

在力學問題中,對系統穩定性的判定非常重要,例如在固體力學問題中,通常需要判定結構的穩定性,而系統的穩定性和系統剛度矩陣的正定性直接相關.對系統的剛度矩陣進行正定性判定,可以直接從在理論上來推測結構穩定性.

對矩陣的正定性進行判定也是一個經典的數學問題.Johnson[1]在1970年提出實正定矩陣的概念及一些判定準則.很多學者對正定矩陣的性質展開了研究工作和總結[2].

經典的判定矩陣的正定性方法主要通過兩種途徑來實現[3-4].第一是從標準型上出發,有若干等同的表述:矩陣K的特征值全大于零,則矩陣K正定;矩陣K與n階單位矩陣En合同,則矩陣K正定;存在可逆矩陣S,使得K=STS,則矩陣K正定等等.第二是從主子式上出發,也有若干等同的表述:矩陣K的所有主子式全大于零,則矩陣K正定;矩陣K的所有順序主子式全大于零,則矩陣K正定等等.

近年來,學者對矩陣正定性的判定展開了一些研究工作.有一些工作基于矩陣分解和變換來判定矩陣的正定性:三角及對角陣LDLT分解和Cholesky分解算法常用于矩陣正定性的判定.李偉[5]通過對矩陣進行初等保號變換,將矩陣轉化為三角矩陣并判斷其中元素的正負來判定原矩陣的正定性.趙淑賢等[6]通過合同變換,證明n階矩陣的正定性可以由n– 1階矩陣的正定性來進行確定,提出使用降階的方法來判定矩陣的正定性.

另有一些部分工作從矩陣的特征值分步的角度來展開.Gerschgorin圓盤定理可用于估計矩陣的特征值分布[7].基于圓盤定理,鄒黎敏等[8]通過比較包含所有特征值的圓盤半徑和圓心到原點的距離,基于優化理論來判定矩陣的正定性.岑燕斌等[9]基于圓盤定理、半正定矩陣的滿秩分解和主對角線嚴格占優判別法,提出極大極小元矩陣正定性判定算法.Sorensen[10]提出使用Arnoldi法來進行多項式濾波,該方法可以用于直接求解矩陣的最小特征值,對矩陣的正定性進行判定.

在有限元問題中,系統的剛度矩陣通常具有很大的規模.當矩陣規模變大之后,使用經典的判定算法所需計算量快速加大,并且這類算法都基于單機開發,難以實現并行化.這時將難以再使用經典的判定方法來判定大規模矩陣的正定性.本文將提出和探討幾種更適合于大規模并行計算框架下矩陣的正定性判定算法.

1 矩陣正定性判定的快速算法

1.1 算法思想的提出

矩陣K的正定性的數學表示是若對于任意的非零列向量x,xTKx始終大于零,則矩陣K是一個正定矩陣.若需判定一個不對稱矩陣的正定性,可以轉化為判定一個對稱矩陣K的正定性,其中

因此下文所討論的矩陣都是對稱矩陣.如前所述,固體力學結構或系統的穩定性與其剛度矩陣的正定性是等價的.因此提出一類算法,為了判定剛度矩陣K的正定性,可以來研究其對應結構的穩定性.具體來說,對于一個n階矩陣K,可以假定“系統的平衡位置”為x0=0,在此基礎上,外界施加一個“擾動”x,其中x是一個n維列向量并代表偏離平衡位置量.之后模擬“系統能量”的下降過程,即使xTKx的數值不斷減小,來模擬力學結構受到擾動后的演化情況,若經過一段時間之后,“系統回到平衡位置”,即x=x0=0,則系統是穩定的,矩陣是正定的.若經過一段時間之后,“系統偏離原來的平衡位置”,即xTKx<0,則系統是不穩定的,即矩陣是非正定的.

因此對矩陣正定性進行判定可以轉化為一個優化問題.借鑒上述力學系統中能量下降的思想,可以基于優化理論對矩陣正定問題進行判定,即

其中,?(x)為目標函數,x為優化變量.求解上述優化問題,若經過一系列的計算之后,能夠找到一個x,使得xTKx<0,則矩陣是非正定的;若能夠找到一個x≠ 0,但xTKx=0,則矩陣是半正定的;若經過大量的計算,始終找不到一個x,使得xTKx<0,則矩陣是正定的.

下文將基于優化理論提出迭代算法來判定矩陣的正定性.給定一個初始擾動,即迭代初值x1之后開始迭代.隨著迭代過程的進行,目標函數的數值不斷下降.在將要進行第i次迭代時,優化變量為xi,迭代結束后,優化變量為xi+1.迭代過程中優化變量滿足如下關系

其中,hi為第i次迭代的前進方向,α為前進步長.

1.2 基于最速下降法的算法一

首先,基于優化理論[11]中的最速下降思想(steepest descent method,SD),提出矩陣正定判定算法一.使用最速下降法進行迭代的示意圖如圖1所示.使用最速下降法進行求解優化問題(2),每次優化變量的前進方向是當前的梯度下降方向,在力學系統中,這相當于對一個無慣性效應的系統進行動力學顯式計算,每次位移將沿著局部不平衡力的方向進行演化.隨著演化的進行,系統的能量逐漸降低.這也是將最直接的最具力學意義的想法進行算法實現.

圖1 最速下降法迭代示意圖Fig.1 Schematic of steepest descent method

迭代過程中優化變量的前進方向即是目標函數梯度的反方向,此時目標函數的梯度方向為

其中,gi為優化變量所處位置目標函數值梯度的負方向.

迭代算法的計算量和優化變量前進步長的選取直接相關.若迭代步長選取的過小,每一次下降的會準確,但是會導致計算量的大幅度增加;若迭代步長選取的過大,在一定程度上,可能導致更快的收斂,但是也可能出現越過極值點而導致的重復計算,同樣造成計算量的增加.建議按如下方式選取α,使得在沿著當前的前進方向hi上,目標函數

達到極小值.此時目標函數值是一個關于α的二次函數.當時,拋物線開口向上,取

此時α使目標函數值?(xi+1)關于α達到極小值,若?(xi+1)>0 暫時不能判定矩陣是否正定,需繼續迭代;而當時,此時目標函數的最小值可趨向–∞,可以顯示矩陣K為非正定的.

下面討論算法迭代結束的終止條件,推薦使用如下3個終止準則.

終止準則一:若能夠計算得到一個非零列向量x,滿足

則此時判定矩陣非正定,迭代終止.

終止準則二:到達最大的迭代次數

若隨著迭代的進行,迭代次數i到達算法容許的最大迭代次數imax,仍然沒有找到非正定方向x,使得xTKx<0,則此時認為矩陣是正定的.算法容許的最大迭代次數imax通常和矩陣規模有關,矩陣規模越大,所需的迭代次數越多.在后文的測試算例中,選取imax=2n.

終止準則三:優化變量方向變化不大

算法的實現過程中主要涉及到矩陣向量乘運算,算法復雜度為O(n2).

1.3 基于共軛梯度法的算法二

需要指出的是,直接基于最速下降法來進行迭代計算,在靠近目標函數極值點的附近,總是曲折前進,優化的效果往往不好,算法收斂速度較慢[10].為了避免這個問題,可以基于前后兩次迭代的前進方向相互共軛的思想,對前進的方向進行修正,使用共軛梯度法 (conjugate gradient method,CG)來判定矩陣的正定性,該方法考慮了上一次的前進方向對此次前進方向的影響,算法的收斂速度將大幅度的加快.在此基于共軛梯度法,提出矩陣正定性判定算法二.使用共軛梯度法進行迭代的示意圖如圖2所示.

圖2 共軛梯度法迭代示意圖Fig.2 Schematic of conjugate gradient method

前進方向落在當前梯度方向和前一次前進方向張開的空間中,滿足

其中β為權系數,將通過共軛方向來確定:前后兩次的前進方向滿足

此時解有

基于共軛梯度算法確定的前進步長和終止準則與最速下降法一致,將由式(7)～式(10)決定.

可以證明,由共軛方向確定的前進方向和前進步長,在當前梯度方向和上一次迭代的前進方向張開的空間中最優的.即通過共軛梯度法確定的α,β滿足

對于規模為n二次優化問題,在理論上使用共軛梯度法迭代n次可以即找到最優解.為了保險起見,選取最大迭代次數imax=Sn,其中S為大于 1 的安全因數.算法的實現過程中涉及到矩陣向量乘運算,算法復雜度為O(n2).

1.4 基于截平面的算法三

另一方面,在實際問題中,經常需要對弱非正定矩陣進行分析.例如在固體力學的問題中,在一個結構是穩定的時候,其剛度是正定的.隨著載荷的增大,結構出現失穩.當結構剛剛出現失穩時,其剛度矩陣一般只有唯一且非常小的負特征值.此時矩陣非正定的非常不明顯.如果通過最速下降算法或者共軛梯度法來進行求解的時候,很可能需要迭代非常多次,才能夠判定矩陣是非正定的,隨著負特征值的絕對值越小,通過若干次搜索找到相應的非正定方向x越困難.因此需要開發一種對弱非正定矩陣有著較好的性能的算法.

這種情況下,提出使用截平面的思想來進行矩陣正定性的判定:在一個n維空間截一個超平面,在這個n–1維的超平面上尋找是否存在一個使得xTKx<0 的點,若存在,則該矩陣非正定;若不存在,則矩陣正定.若系統穩定性剛剛發生轉化,轉變為一個不穩定的系統,此時其剛度矩陣只有一個非常小的負特征值.此時只有沿著某一個特定的方向,目標函數xTKx才小于0,如圖3空間中細橢圓錐所示.所截的這個超平面內只有一小塊區域(即圖3中藍色區域)xTKx<0,此時其中必包含一個極小值點,當然也是超平面上函數xTKx的駐值點.

綜上所述,對于弱非正定矩陣的判定,提出可以先截取一個超平面然后求其駐值點.并稱這第三種算法為截平面法 (cutting plane method,CP).具體做法是,選擇一個不通過原點的超平面,其過原點的法向量記為d,則該超平面可以表示為對于弱非正定矩陣而言,認為所截的超平面上的點的極小值點就是駐值點.因此構建拉格朗日函數

圖3 針對弱非正定矩陣的截平面法示意圖Fig.3 Schematic of cutting plane method for weak non-positive definite matrix

其中λ為拉格朗日乘子,對上式取駐值有

求上述超平面的駐值問題轉化為一個n+1維的線性方程組求解問題,其解x就對應于駐值點.由于該線性方程組的系數矩陣滿秩,因此有唯一的解或駐值點存在.若剛度矩陣只有一個小負特征值的時候,如圖3所示,此時該駐值點x將極有可能是超平面中目標函數的極小值點,所以用這個方法一次求解就可以發現非正定方向.

需要指出的是,固體結構求解中的線性方程組矩陣K一般是稀疏的,而易知式(17)中的系數矩陣依然稀疏,充分利用矩陣的稀疏性,可以使算法復雜度接近O(n).在求解式(17)的線性方程組的過程中,如果是超大規模的問題則必須采用并行算法來求解,但線性方程組直接解法(如LU分解或LDLT分解)的并行求解規模難以擴大而不能被采用,因此這些直接算法無法用于確定矩陣的正定性,需采用更適用于并行的迭代算法,如多級平衡的并行求解器[12].

需要說明的是,截平面法對于一些弱非正定的矩陣一般比較有效.倘若超平面(或法向量d)選取不當,或者系統不穩定方向的圓錐開口較大時,如圖 4所示.此時超平面內xTKx<0的區域如圖 4中藍色區域所示,xTKx>0 的區域如圖 4 中淺黃色區域所示,此時超平面上函數xTKx的駐值點將不再是極小值點,此時僅僅通過截平面求駐值點來判定矩陣的正定性可能會失效.在實際的應用中,為了判定一個弱非正定矩陣的正定性,為了保險起見,可以截多個超平面來進行測試,這樣就更有可能超平面上的駐值點就是極小值點,以保證截平面法求解的準確性.

圖4 針對強非正定矩陣的截平面法示意圖Fig.4 Schematic of cutting plane method for strong non-positive definite matrix

1.5 算法四: 高效混雜使用上述算法

如前所述,前面各種方法都有優缺點和較適用的情形.如對于一些強非正定矩陣而言采用截平面法時,可能需使用多個法向量來測試才能使得超平面上的駐值點為極小值點.為了避免這種弊端,可以把截平面法和共軛梯度法相結合:固定優化變量只能在所截的超平面附近運動,把對矩陣正定性的判定問題轉化為一個優化問題,首先使用截平面法來確定迭代的初值,再使用共軛梯度法來迭代求解.這樣可以同時保證對弱非正定矩陣進行判定的高效性,及對強非正定矩陣判定的準確性.在此基于截平面法 (cutting plane method,CP) 和共軛梯度法(conjugate gradient method,CG)提出混雜算法四(CP&CG),按以下三部分論述.

1.5.1 優化問題的提出

假設存在一個向量x,其末端落在所截超平面中,則向量x–d是這個平面中的向量,故原來的式(1)對應的無約束優化問題轉化為一個有約束的優化問題

使用罰函數法對該問題進行求解,把有約束優化問題轉化為無約束優化問題

1.5.2 迭代初值的選取

一個好的迭代初值可以使目標函數值快速下降,使用更少的迭代次數就能達到較小的目標函數值.迭代的初值將選取在所截超平面的駐值點上,即由截平面法通過求解式(17)得到.

1.5.3 優化問題求解

經過求解式(17)對應的線性方程組,找到了一個較好的迭代初值,這將使后續的優化計算有更加快速的收斂速度.倘若此時超平面上的駐值點是極小值點,迭代的初始值處xTKx<0,此時就不再需要使用共軛梯度法進行后續的迭代,可以直接判定矩陣是非正定的,此時CP&CG算法即退化為CP算法.倘若此時超平面上的駐值點不是極小值點,如圖4所示,則使用共軛梯度法進行后續的迭代,尋找在超平面附近是否存在使得xTKx<0的點.

2 算法性能測試

前文提出了4種判定矩陣正定性的快速算法:基于最速下降法的算法一(SD)、基于共軛梯度法的算法二(CG)、基于直接截平面的算法三(CP)和混雜算法四(CP&CG).在本節中,將生成不同類型的矩陣,來測試這些算法的性能.考慮到力學系統中總體剛度矩陣稀疏的特點,分別討論算法對稠密矩陣和稀疏矩陣的性能.

2.1 測試所用矩陣的生成方法

使用三個指標來衡量測試所用的對稱矩陣:矩陣規模n、負特征值比例p和負特征值大小指標q.測試所用的對稱矩陣通過下式生成

其中W是大小為n×n隨機矩陣,其中的每一個非零元素在[–0.5,0.5]之間服從均勻分布.D是大小為n×n的對角矩陣,D中元素分成兩類,一類為正數,大小為 1,一類為負數,大小為q.

當負特征值比例p=0時,對角陣D中所有元素都是正數,此時生成的測試矩陣K是正定的.當負特征值比例p=1時,對角陣中的所有元素都是負數,此時矩陣所有特征值都是負數,此時生成的測試矩陣是負定的.在負特征值比例p∈(0,1)時,測試矩陣K是非正定的,此時負特征值比例p越小,負特征值大小指標q的絕對值越小,矩陣非正定的越不明顯.

若隨機矩陣W是一個稠密矩陣,經過式(21)得到測試矩陣也是稠密的.在力學問題中,系統的剛度矩陣往往是稀疏的,選取稀疏的矩陣作為隨機矩陣W,經過式(21)得到的測試矩陣也是稀疏的.本文測試所用的稀疏矩陣每一行中非零元素占比約為10%.

2.2 強非正定矩陣測試結果

下討論3種算法對非正定矩陣的測試結果,并和經典的分解算法LDLT進行對比.

選取負特征值比例p=0.5,負特征值大小指標q=-1,隨機生成一系列強非正定矩陣進行測試.測試矩陣的規模為 256 階,重復試驗 1 000 次,選擇算法參數幾種算法的測試結果如表1所示.由于截平面算法CP主要針對弱非正定矩陣,此時不再對該算法進行測試.

從強非正定矩陣測試結果可以看出,SD算法耗時最短,CG 算法其次,CP&CG 算法隨后.SD 與CG算法用時相近,明顯短于CP&CG算法.這是對強非正定矩陣進行判定時,優化算法迭代幾次之后很快就能夠找到一個優化變量x,使得xTKx<0,從而判定矩陣是非正定的;而CP&CG算法為了選取一個較好的迭代初始值,需要對一個n+1階的線性方程組進行求解,這導致了額外的計算量.此時如果直接使用CP算法,超平面上的駐值點一般都不是極小值點,CP算法將會失效,而把共軛梯度法和截平面算法相結合的CP&CG混雜算法保證了對強非正定矩陣判定的準確性.

表1 不同算法對強非正定矩陣的判定結果對比Table 1 Algorithm performance for strong non-positive matrix

對于強非正定矩陣而言,本文提出的基于優化理論的矩陣正定判定算法的耗時遠小于經典的LDLT分解算法,我們可以作如下解釋:在分解算法的計算過程中,需要把矩陣先進行LU分解,再進行LDLT分解.等到整個分解過程完全結束之后才能夠對對角矩陣中的元一一觀察,從而判斷矩陣是否正定,這整個分解過程相當于把整個矩陣“完全”地分析了,與之不同的是,優化算法通過演化的形式,只需要找到“一個”x,使得xTKx<0,即可給出判定結果.與分解算法的“完全”分析相比,優化算法的“局部”分析顯然更加的快速,固在大規模情況下,優化算法的計算效率將會更高.

2.3 弱非正定矩陣測試結果

下面討論這些算法對弱非正定矩陣的測試結果,并和經典的分解算法LDLT進行對比.

隨機生成一系列弱非正定矩陣來進行測試,對角陣D中只有一個元素是負數,負特征值大小指標q= -10–5,這使得測試的矩陣只有唯一特別小的特征值.測試矩陣的規模為256階,重復試驗1000次,選擇算法參數S=2,ε=1×10–12,幾種算法的測試結果如下表2所示.對弱非正定矩陣而言,由于CP&CG混雜算法直接求解線性方程組作為迭代初值之后即可以判定矩陣的正定性,不再需要進行后續的優化迭代,此時CP&CG算法退化為CP算法,因此不再需要對單獨CP&CG算法進行測試.

表2 不同算法對弱非正定矩陣的判定結果對比Table 2 Algorithm performance for weak non-positive matrix

從表2中可以看出,對于弱非正定矩陣而言,LDLT,CG和CP算法仍然能夠準確的判定其正定性,而SD算法由于其較慢的收斂速度而失效.此時最優的是CP算法,其通過求解n+ 1階的線性方程組,通過計算超平面上的駐值點來尋找超平面上極小值點,直接判定出矩陣的正定性.而CG算法則需要迭代很多次之后才能找到一個優化變量x,使得xTKx<0,從而判定矩陣的正定性,并且矩陣的負特征值越小,使用CG算法來尋找非正定方向就更加的困難.

與稠密矩陣相比,稀疏矩陣由于其進行矩陣向量乘法和線性方程組求解的時候所需的計算量更少,此時優化算法的計算效率加快.

2.4 測試結果討論

綜上所述,最速下降算法SD由于其收斂速度較慢,對于一些弱非正定矩陣不能準確的判定,在實際運用中舍去.截平面法CP對弱非正定矩陣進行判定時有很好的性能,可以和共軛梯度法混雜一起使用,以保證對強非正定矩陣的判定的準確性.不管是稠密矩陣還是稀疏矩陣,共軛梯度算法和混雜算法對矩陣的正定性進行能夠準確的判定.對弱非正定矩陣而言,混雜算法CP&CG有著更好的性能;對強非正定矩陣而言,共軛梯度法CG有著更好的性能.

共軛梯度算法、截平面算法和混雜算法都容易進行程序的編寫.同時算法涉及的主要的操作是矩陣向量乘法和線性方程組的迭代求解計算,對于一些大規模的問題容易實現并行化.需要指出的是,在力學問題中,系統的剛度矩陣往往是非常稀疏的.此時求解大規模的線性方程組的算法復雜度接近O(n).并且當結構將要發生失穩的時候,系統的剛度矩陣只有一個很小的負特征值.此時對于弱非正定系統而言,此時使用截平面算法CP和混雜算法CP&CG來進行判定矩陣的正定性將會更加地高效.

3 結論

本文主要探討了對于固體結構分析中超大規模矩陣的正定性判別算法,由于必須采用并行算法,所以一些如求特征值或LDLT分解算法都無法適用,因此探討了一些可以用于迭代并行的算法.主要結論如下.

(1)借鑒力學系統中能量下降的思想,提出一種判定矩陣正定性的新思路,即將矩陣的正定性判定問題轉化為一個優化問題.并基于優化理論提出二種判定矩陣正定性的快速算法:最速下降法和共軛梯度法.

(2)考慮到力學系統中,結構剛剛失穩狀態時的剛度矩陣的稀疏性和弱非正定性,提出截平面法來判定矩陣的正定性,即先截超平面后求解線性方程組找駐值點.

(3)提出可以混雜使用共軛梯度法和截平面法,在高效判定弱非正定矩陣正定性的同時,保證對強非正定矩陣的判定準確性.

(4)共軛梯度算法能夠更加快速的判定強非正定矩陣的正定性.截平面法和混雜算法能夠更加快速的判定弱非正定矩陣的正定性.

1 Johnson CR.Positive definite matrices.American Mathematical Monthly,1970,77(3):259-264

2 Bhatia R.Positive Definite Matrices.New Jersey:Princeton University Press,2015

3 屠伯塤,徐誠浩,王芬.高等代數.上海:上海科學技術出版社,1984(Tu Boxun,Xu Chenghao,Wang Fen.Advanced Algebra.Shanghai:Shanghai Science and Technology Press,1984 (in Chinese))

4 張禾瑞,郝炳新.高等代數.第4版.北京:高等教育出版社,2001(Zhang Herui,Hao Bingxin.Advanced Algebra.4th ed.Beijing:Higher Education Press,2011(in Chinese))

5 李偉.關于正定矩陣判別定理的改進.高等數學研究,2007,10(6):42-43(Li Wei.Improvement of positive definite matrix discriminant theorem.Studies in College Mathematics,2007,10(6):42-43(in Chinese))

6 趙賢淑,楊莉軍.正定矩陣的降階判別法及其應用.北京印刷學院學報,2000(2):40-43(Zhao Xianshu,Yang Lijun.Dicrimination method of reduction of order of positive definite matrix and its application.Journal of Beijing Institute of Graphic Communication,2000(2):40-43(in Chinese))

7 Varga RS.Ger?gorin and His Circle.Berlin Heidelberg:Springer 2004

8 鄒黎敏,胡興凱,伍俊良.正定矩陣的性質及判別法.中山大學學報自然科學版,2009,48(5):16-23(Zou Limin,Hu Xingkai ,Wu Junliang.The properties and discrimination of the positive definite matrices.Acta Scientiarum Naturalium Universitatis Sunyatseni,2009,48(5):16-23(in Chinese))

9 岑燕斌,韋煜,羅會亮.快速判斷一類實對稱矩陣正定的極大極小元方法.北京交通大學學報,2011,35(6):140-143(Cen Yanbin,Wei Yu,Luo Huiliang.Max and mini-element method:A rapid determination of one class rear symmetric positive definite matrices.Journal of Beijing Jiaotong University,2011,35(6):140-143(in Chinese))

10 Sorensen DC.Implicit application of polynomial filters in a k-step Arnoldi method.Society for Industrial and Applied Mathematics,1992

11 Boyd,Vandenberghe,Faybusovich.Convex optimization.IEEE Transactions on Automatic Control,2006,51(11):1859

12 Xu R,Liu B,Dong Y.Scalable hierarchical parallel algorithm for the solution of super large-scale sparse linear equations.Journal of Applied Mechanics,2012,80(2):a111-121

FAST ALGORITHMS FOR DETERMINING POSITIVE DEFINITENESS OF LARGE SCALE MATRICES IN SOLID STRUCTURE ANALYSIS1)

Xu Hui Liu Bin2)
(AML,Department of Engineering Mechanics,Tsinghua University,Beijing100084,China)

In order to determine the positive definiteness of the super-large-scale matrix in the structural stability analysis,the parallel computing method must be adopted.However,the traditional direct methods such as the eigenvalue’s analysis,the master subordinate determinant’s computation and LDLT decomposition are difficult to realize in parallel computing.In this paper,some iterative algorithms which are suitable for parallel computing are proposed.A new approach is developed that determining the positive definiteness of a matrix can be transformed into an optimization problem,which is solved by various optimization algorithms.The algorithms based on the steepest descent method and the conjugate gradient method are proposed.Considering the sparseness of the stiffness matrix of the mechanical system and the weakly non-positive definite property of at the critical buckling state,we propose a method via calculating the stationary point on a cutting plane to determine the non-positive definite matrices.In order to ensure the accuracy of the determination of strong non-positive definite matrices,a hybrid method combing the cutting plane method and the conjugate gradient method is developed.The numerical results show that the proposed algorithms are accurate and efficient.The conjugate gradient algorithm is more efficient for strongly non-positive definite matrices while the hybrid method is more efficient for the weak non-positive definite matrices.These algorithms are easy to realize in parallel computing on the cluster and can quickly determine the positive definiteness of large-scale matrices.

matrix,positive definite,conjugate gradient method,cutting plane

O346.1

A

10.6052/0459-1879-17-292

2017–08–26 收稿,2017–11–15 錄用,2017–11–15 網絡版發表.

1)國家自然科學基金資助項目 (51232004,11372158,11425208).

2)劉彬,教授,主要研究方向:大規模多尺度多物理場計算、復合材料力學和斷裂力學.E-mail:liubin@tsinghua.edu.cn

徐輝,劉彬.固體結構大規模矩陣正定性判定的快速算法.力學學報,2017,49(6):1223-1230

Xu Hui,Liu Bin.Fast algorithms for determining positive definiteness of large scale matrices in solid structure analysis.Chinese Journal of Theoretical and Applied Mechanics,2017,49(6):1223-1230