固體結(jié)構(gòu)大規(guī)模矩陣正定性判定的快速算法1)

徐 輝 劉 彬

(清華大學(xué)航天航空學(xué)院工程力學(xué)系應(yīng)用力學(xué)教育部重點(diǎn)實(shí)驗(yàn)室,北京 100084)

固體結(jié)構(gòu)大規(guī)模矩陣正定性判定的快速算法1)

徐 輝 劉 彬2)

(清華大學(xué)航天航空學(xué)院工程力學(xué)系應(yīng)用力學(xué)教育部重點(diǎn)實(shí)驗(yàn)室,北京 100084)

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

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

引言

在力學(xué)問題中,對系統(tǒng)穩(wěn)定性的判定非常重要,例如在固體力學(xué)問題中,通常需要判定結(jié)構(gòu)的穩(wěn)定性,而系統(tǒng)的穩(wěn)定性和系統(tǒng)剛度矩陣的正定性直接相關(guān).對系統(tǒng)的剛度矩陣進(jìn)行正定性判定,可以直接從在理論上來推測結(jié)構(gòu)穩(wěn)定性.

對矩陣的正定性進(jìn)行判定也是一個(gè)經(jīng)典的數(shù)學(xué)問題.Johnson[1]在1970年提出實(shí)正定矩陣的概念及一些判定準(zhǔn)則.很多學(xué)者對正定矩陣的性質(zhì)展開了研究工作和總結(jié)[2].

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

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

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

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

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

1.1 算法思想的提出

矩陣K的正定性的數(shù)學(xué)表示是若對于任意的非零列向量x,xTKx始終大于零,則矩陣K是一個(gè)正定矩陣.若需判定一個(gè)不對稱矩陣的正定性,可以轉(zhuǎn)化為判定一個(gè)對稱矩陣K的正定性,其中

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

因此對矩陣正定性進(jìn)行判定可以轉(zhuǎn)化為一個(gè)優(yōu)化問題.借鑒上述力學(xué)系統(tǒng)中能量下降的思想,可以基于優(yōu)化理論對矩陣正定問題進(jìn)行判定,即

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

下文將基于優(yōu)化理論提出迭代算法來判定矩陣的正定性.給定一個(gè)初始擾動(dòng),即迭代初值x1之后開始迭代.隨著迭代過程的進(jìn)行,目標(biāo)函數(shù)的數(shù)值不斷下降.在將要進(jìn)行第i次迭代時(shí),優(yōu)化變量為xi,迭代結(jié)束后,優(yōu)化變量為xi+1.迭代過程中優(yōu)化變量滿足如下關(guān)系

其中,hi為第i次迭代的前進(jìn)方向,α為前進(jìn)步長.

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

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

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

迭代過程中優(yōu)化變量的前進(jìn)方向即是目標(biāo)函數(shù)梯度的反方向,此時(shí)目標(biāo)函數(shù)的梯度方向?yàn)?/p>

其中,gi為優(yōu)化變量所處位置目標(biāo)函數(shù)值梯度的負(fù)方向.

迭代算法的計(jì)算量和優(yōu)化變量前進(jìn)步長的選取直接相關(guān).若迭代步長選取的過小,每一次下降的會(huì)準(zhǔn)確,但是會(huì)導(dǎo)致計(jì)算量的大幅度增加;若迭代步長選取的過大,在一定程度上,可能導(dǎo)致更快的收斂,但是也可能出現(xiàn)越過極值點(diǎn)而導(dǎo)致的重復(fù)計(jì)算,同樣造成計(jì)算量的增加.建議按如下方式選取α,使得在沿著當(dāng)前的前進(jìn)方向hi上,目標(biāo)函數(shù)

達(dá)到極小值.此時(shí)目標(biāo)函數(shù)值是一個(gè)關(guān)于α的二次函數(shù).當(dāng)時(shí),拋物線開口向上,取

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

下面討論算法迭代結(jié)束的終止條件,推薦使用如下3個(gè)終止準(zhǔn)則.

終止準(zhǔn)則一:若能夠計(jì)算得到一個(gè)非零列向量x,滿足

則此時(shí)判定矩陣非正定,迭代終止.

終止準(zhǔn)則二:到達(dá)最大的迭代次數(shù)

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

終止準(zhǔn)則三:優(yōu)化變量方向變化不大

算法的實(shí)現(xiàn)過程中主要涉及到矩陣向量乘運(yùn)算,算法復(fù)雜度為O(n2).

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

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

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

前進(jìn)方向落在當(dāng)前梯度方向和前一次前進(jìn)方向張開的空間中,滿足

其中β為權(quán)系數(shù),將通過共軛方向來確定:前后兩次的前進(jìn)方向滿足

此時(shí)解有

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

可以證明,由共軛方向確定的前進(jìn)方向和前進(jìn)步長,在當(dāng)前梯度方向和上一次迭代的前進(jìn)方向張開的空間中最優(yōu)的.即通過共軛梯度法確定的α,β滿足

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

1.4 基于截平面的算法三

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

這種情況下,提出使用截平面的思想來進(jìn)行矩陣正定性的判定:在一個(gè)n維空間截一個(gè)超平面,在這個(gè)n–1維的超平面上尋找是否存在一個(gè)使得xTKx<0 的點(diǎn),若存在,則該矩陣非正定;若不存在,則矩陣正定.若系統(tǒng)穩(wěn)定性剛剛發(fā)生轉(zhuǎn)化,轉(zhuǎn)變?yōu)橐粋€(gè)不穩(wěn)定的系統(tǒng),此時(shí)其剛度矩陣只有一個(gè)非常小的負(fù)特征值.此時(shí)只有沿著某一個(gè)特定的方向,目標(biāo)函數(shù)xTKx才小于0,如圖3空間中細(xì)橢圓錐所示.所截的這個(gè)超平面內(nèi)只有一小塊區(qū)域(即圖3中藍(lán)色區(qū)域)xTKx<0,此時(shí)其中必包含一個(gè)極小值點(diǎn),當(dāng)然也是超平面上函數(shù)xTKx的駐值點(diǎn).

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

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

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

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

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

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

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

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

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

1.5.1 優(yōu)化問題的提出

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

使用罰函數(shù)法對該問題進(jìn)行求解,把有約束優(yōu)化問題轉(zhuǎn)化為無約束優(yōu)化問題

1.5.2 迭代初值的選取

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

1.5.3 優(yōu)化問題求解

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

2 算法性能測試

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

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

使用三個(gè)指標(biāo)來衡量測試所用的對稱矩陣:矩陣規(guī)模n、負(fù)特征值比例p和負(fù)特征值大小指標(biāo)q.測試所用的對稱矩陣通過下式生成

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

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

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

2.2 強(qiáng)非正定矩陣測試結(jié)果

下討論3種算法對非正定矩陣的測試結(jié)果,并和經(jīng)典的分解算法LDLT進(jìn)行對比.

選取負(fù)特征值比例p=0.5,負(fù)特征值大小指標(biāo)q=-1,隨機(jī)生成一系列強(qiáng)非正定矩陣進(jìn)行測試.測試矩陣的規(guī)模為 256 階,重復(fù)試驗(yàn) 1 000 次,選擇算法參數(shù)幾種算法的測試結(jié)果如表1所示.由于截平面算法CP主要針對弱非正定矩陣,此時(shí)不再對該算法進(jìn)行測試.

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

表1 不同算法對強(qiáng)非正定矩陣的判定結(jié)果對比Table 1 Algorithm performance for strong non-positive matrix

對于強(qiáng)非正定矩陣而言,本文提出的基于優(yōu)化理論的矩陣正定判定算法的耗時(shí)遠(yuǎn)小于經(jīng)典的LDLT分解算法,我們可以作如下解釋:在分解算法的計(jì)算過程中,需要把矩陣先進(jìn)行LU分解,再進(jìn)行LDLT分解.等到整個(gè)分解過程完全結(jié)束之后才能夠?qū)蔷仃囍械脑灰挥^察,從而判斷矩陣是否正定,這整個(gè)分解過程相當(dāng)于把整個(gè)矩陣“完全”地分析了,與之不同的是,優(yōu)化算法通過演化的形式,只需要找到“一個(gè)”x,使得xTKx<0,即可給出判定結(jié)果.與分解算法的“完全”分析相比,優(yōu)化算法的“局部”分析顯然更加的快速,固在大規(guī)模情況下,優(yōu)化算法的計(jì)算效率將會(huì)更高.

2.3 弱非正定矩陣測試結(jié)果

下面討論這些算法對弱非正定矩陣的測試結(jié)果,并和經(jīng)典的分解算法LDLT進(jìn)行對比.

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

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

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

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

2.4 測試結(jié)果討論

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

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

3 結(jié)論

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

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

(2)考慮到力學(xué)系統(tǒng)中,結(jié)構(gòu)剛剛失穩(wěn)狀態(tài)時(shí)的剛度矩陣的稀疏性和弱非正定性,提出截平面法來判定矩陣的正定性,即先截超平面后求解線性方程組找駐值點(diǎn).

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

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

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 屠伯塤,徐誠浩,王芬.高等代數(shù).上海:上?？茖W(xué)技術(shù)出版社,1984(Tu Boxun,Xu Chenghao,Wang Fen.Advanced Algebra.Shanghai:Shanghai Science and Technology Press,1984 (in Chinese))

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

5 李偉.關(guān)于正定矩陣判別定理的改進(jìn).高等數(shù)學(xué)研究,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 趙賢淑,楊莉軍.正定矩陣的降階判別法及其應(yīng)用.北京印刷學(xué)院學(xué)報(bào),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 鄒黎敏,胡興凱,伍俊良.正定矩陣的性質(zhì)及判別法.中山大學(xué)學(xué)報(bào)自然科學(xué)版,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 岑燕斌,韋煜,羅會(huì)亮.快速判斷一類實(shí)對稱矩陣正定的極大極小元方法.北京交通大學(xué)學(xué)報(bào),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 網(wǎng)絡(luò)版發(fā)表.

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

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

徐輝,劉彬.固體結(jié)構(gòu)大規(guī)模矩陣正定性判定的快速算法.力學(xué)學(xué)報(bào),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