使用 uBLAS 进行实对称正定矩阵的 Cholesky 分解

最新推荐文章于 2022-12-22 09:58:33 发布

五道口纳什

最新推荐文章于 2022-12-22 09:58:33 发布

阅读量2.5k

点赞数

分类专栏： C/C++ 矩阵分析

本文链接：https://blog.csdn.net/lanchunhui/article/details/50897654

版权

C/C++ 同时被 2 个专栏收录

212 篇文章 20 订阅

订阅专栏

矩阵分析

57 篇文章 12 订阅

订阅专栏

Cholesky 分解理论

矩阵分解——三角分解（Cholesky 分解）
矩阵分解——三角分解（二）

注：只有实对称矩阵才有 Cholesky 分解理论。

已知实对称正定矩阵 $A$ ，其 Cholesky 分解形式为： $A=LL^T$ ， $L$ 为下三角矩阵，计算 $L=(\ell_{ij})$ 的递推公式为：

ℓ i j = ⎧ ⎩ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎛ ⎝ a i j - \sum k = 1 j - 1 ℓ 2 i k ⎞ ⎠ 1 / 2, 1 ℓ j j ⎛ ⎝ a i j - \sum k = 1 j - 1 ℓ i k ℓ j k ⎞ ⎠, 0, i = j i > j i < j

$\ell_{ij}=\left \{ \begin{array}{ll} \left (a_{ij}-\sum_{k=1}^{j-1}\ell_{ik}^2\right )^{1/2},&i=j\\ \frac1{\ell_{jj}}\left ( a_{ij}-\sum_{k=1}^{j-1}\ell_{ik}\ell_{jk}\right ),&i>j\\ 0,&i<j \end{array} \right.$

uBLAS 相关线性代数 API

uBLAS——Boost 线性代数基础程序库
 uBLAS——Boost 线性代数基础程序库（二）

（1）内积

ublas::inner_prod()
（2）矩阵切片——取出矩阵的行

ublas::row()
（3）向量切片——取出向量的部分元素

ublas::range()

C++ 实现

void chol(const ublas::symmetric_matrix<double>& A, ublas::triangular_matrix<double>& L)
{
    for (unsigned i = 0; i < A.size1(); ++i)
    {
        for (unsigned j = 0; j <= i; ++j)
        {
            if (i == j)
                L(i, i) = std::pow(A(i, i) - ublas::inner_prod(ublas::project(ublas::row(L, i), ublas::range(0, i)),
                ublas::project(ublas::row(L, i), ublas::range(0, i))), 1. / 2);
            L(i, j) = (A(i, j) - ublas::inner_prod(ublas::project(ublas::row(L, i), ublas::range(0, j)),
                ublas::project(ublas::row(L, j), ublas::range(0, j)))) / L(j, j);
        }
    }
}