首先分析下矩陣是否存在特殊性質,是否稀疏、是否對稱

其次尋求對應的求解方案

個人認為稀疏對稱正定矩陣最容易求解(共軛梯度法,預調節器用Incomplete Cholesky),

其次是稀疏矩陣,LLT LDLT QR LU都有對應的稀疏分解方案;大規模矩陣下用迭代法效率(最小二乘LSCG、BiCGSTAB等)更佳,不過需要選擇合適的預調節器

最糟糕的情況是稠密矩陣,儘管稠密對稱矩陣也有對應的求解方案(如Cholesky分解),但是求解速度依舊感人。

Matlab對矩陣求逆的支持已經很不錯了,我的建議是能否通過某些轉換將待求逆矩陣轉換成具有稀疏對稱性質的新矩陣,同時,盡量不要再去求逆矩陣,而是求解一個線性方程組(例如牛頓法里需要求H^{-1}f,可以轉換成求解線性方程組Hx=f),這樣可以減少存儲空間。

可以嘗試Eigen,通過閱讀裡面提及的benchmark來選擇合適的求解工具(第三方庫)。

Eigen: Main Page?

eigen.tuxfamily.org

作為參考,我之前求解一個10^5,稀疏度(非零元佔比)為1%的稀疏對稱矩陣,求解平均耗時9秒左右,思路是用Incomplete Cholesky的共軛梯度法,Eigen有提供相關的介面(Eigen的速度會慢一些,我是採用了他的思路,並且根據求解矩陣的特有性質和應用場景進行了一些優化,最後是自己用C++重新造了小小的輪子)。

Eigen::ConjugateGradient&< _MatrixType, _UpLo, _Preconditioner &> Class Template Reference?

eigen.tuxfamily.org圖標

希望可以幫到你


如果是一般的矩陣求逆,沒有什麼太好的辦法,如果是正定矩陣,可以用cholesky decompisition。如果矩陣的rank又比較低,會有一些row rank的方法。


迭代大法


matlab的inv也是變成線性方程組解的,你這個問題相當於81E8規模的計算量,明顯工作量太大。你只能根據你的矩陣的特點自己看看能不能利用稀疏性什麼的減工作量。


每10000*10000格視為 1區, 每一區取其內含格子的平均數,做為代表數。 於是一個 90000*90000的矩陣化為9*9的矩陣, 機器可解之。

再改成5000*5000格視為 1區,機解之。

再改成2500*2500格視為 1區,機解之。

觀看解之演變趨勢,猜測終極解答。

--- --- --- ---

或是視300*300格為一區, 機解出300*300的小逆矩陣,

再把小逆矩陣里的aij數字改寫為相同的300*300個aij數字,於是得90000*90000的概略逆矩陣,

再把概略逆矩陣里的a11原佔地盤改為300*300個未知數,讓概略逆矩陣與原矩陣相乘,企圖得90000*90000的單位矩陣,但我們只專心看最左上方300*300個數,它要等於300*300的單位矩陣,於是叫機器解開此300*300個未知數。

接著把概略逆矩陣里的a22原佔地盤改為300*300個未知數,機解之。

接著把概略逆矩陣里的a33原佔地盤改為300*300個未知數,機解之。


gpu加速不知道行不行?


要考慮矩陣是否有可利用的結構特徵,比如稀疏矩陣,如果有,可以採用一些快速演算法。


用迭代法,還得看這矩陣有沒有特殊性質,如果是稀疏矩陣就可以用spdiags的命令,如果是正定矩陣可以用楚列斯基分解或共軛梯度法。


可以考慮矩陣近似方法,通過對小矩陣進行svd分解,來近似原始矩陣的分解形式


改成並行演算法


存下這個矩陣大約要 64GB 的內存。


看你是什麼矩陣了,稠密矩陣一般沒有特別好的辦法。

如果是稀疏矩陣,那就好辦多了,有很多現成的求解器可以用。直接求解器比如Pardiso,如果是對稱正定的稀疏矩陣,開了OMP並行後,估計1s內算出來,如果在多核伺服器上還能更快,不知道能不能滿足你的要求。如果對稱不定,速度稍微慢點,但也不會超過2s。

迭代求解器考慮PESTc,不過我覺的你這個規模的稀疏矩陣沒必要用迭代求解器,規模太小了。直接求解器足夠了。


用intel庫求解或者調用pytorch,TensorFlow之類的支持GPU的軟體,另外Julia可能也可以幫到你


推薦閱讀:
相关文章