在C++中半定矩陣的Cholesky分解

Question

我要實現一個函數來在C++中完成半定矩陣的Cholesky分解，並且想知道是否存在已經優化過的庫/任何東西。這是因爲工作或類似於在此描述：

http://www.sciencedirect.com/science/article/pii/S0096300310012713

這是正定的一個例子，但它並不適用於半正定工作：http://en.wikipedia.org/wiki/Cholesky_decomposition#The_Cholesky.E2.80.93Banachiewicz_and_Cholesky.E2.80.93Crout_algorithms

程序必須不用C/FORTRAN庫，（想想尖尖的老闆給出指令），這意味着ATLAS，LAPACK等。出局。我通過MTL + Boost看過，但他們只適用於正定矩陣。有沒有我沒有找到的圖書館，甚至沒有寫過的單一功能？

Answer 1

半決定矩陣的Cholesky分解的問題是1）它不是唯一的2）Crout的算法失敗。

分解的存在通常是證明非建設性，通過一個限制參數（如果M_n - > M，並且M_n = U^t_n U_n，然後|| || U_n = || || M_n^1/2其中||。|| = Hilbert-Schmidt範數，U_n是一個有界序列，提取一個子序列以找到一個滿足U^t U = M和U三角形的極限U.）

我發現在情況我感興趣的是，將對角線元素乘以1 +ε，使ε小（獲得機器ε的幾千倍）以給出完全可接受的分解，這是令人滿意的。事實上，如果M爲正半定，那麼對於每個ε> 0，M +εI是確定的正。

由於當epsilon變爲零時該方案收斂，因此可以考慮計算多個epsilon的分解並執行Richardson外推。

至於肯定的情況，你可以自己實現Crout的算法（在Numerical Recipes中有一個示例代碼），但是我強烈建議不要自己編寫它，並建議使用LAPACK來代替。

這可能涉及讓您的老闆支付英特爾MKL的費用，如果他擔心LAPACK潛在的糟糕實施。大多數時候我聽到這樣的講話，其理由是「但是我們無法控制代碼，我們希望自己編寫代碼，以便在出現問題時進行調試」。愚蠢的說法。 LAPACK已有40年曆史並經過全面測試。

要求不使用LAPACK是要求不使用的正弦，餘弦和對數標準庫一樣傻。