概念

LDL分解是经典Cholesky分解的一个变形：

• L 为下三角矩阵，且对角元素必须为1；
• D为对角矩阵。

转换

特征

• LDL变形如果得以有效运行，构造及使用时所需求的空间及计算的复杂性与经典Cholesky分解是相同的，但是可避免提取平方根。
• 某些不存在Cholesky分解的不定矩阵，也可以运行LDL分解，此时矩阵中会出现负数元素。

因此人们更倾向于使用LDL分解。对于实数矩阵，该种分解的形式可被改写成

计算

减少乘法运算，引入辅助量,则公式为：

代码

#include <iostream>
#include <string.h>
#include <stdio.h>
#include <vector>
#include <math.h>
 
using namespace std;
const int N = 1005;
typedef double Type;
 
Type A[N][N], L[N][N];
 
/** 分解A得到A = L * L^T */
void Cholesky(Type A[][N], Type L[][N], int n)
{
    for(int k = 0; k < n; k++)
    {
        Type sum = 0;
        for(int i = 0; i < k; i++)
            sum += L[k][i] * L[k][i];
        sum = A[k][k] - sum;
        L[k][k] = sqrt(sum > 0 ? sum : 0);
        for(int i = k + 1; i < n; i++)
        {
            sum = 0;
            for(int j = 0; j < k; j++)
                sum += L[i][j] * L[k][j];
            L[i][k] = (A[i][k] - sum) / L[k][k];
        }
        for(int j = 0; j < k; j++)
            L[j][k] = 0;
    }
}
 
/** 回带过程 */
vector<Type> Solve(Type L[][N], vector<Type> X, int n)
{
    /** LY = B  => Y */
    for(int k = 0; k < n; k++)
    {
        for(int i = 0; i < k; i++)
            X[k] -= X[i] * L[k][i];
        X[k] /= L[k][k];
    }
    /** L^TX = Y => X */
    for(int k = n - 1; k >= 0; k--)
    {
        for(int i = k + 1; i < n; i++)
            X[k] -= X[i] * L[i][k];
        X[k] /= L[k][k];
    }
    return X;
}
 
void Print(Type L[][N], const vector<Type> B, int n)
{
    for(int i = 0; i < n; i++)
    {
        for(int j = 0; j < n; j++)
            cout<<L[i][j]<<" ";
        cout<<endl;
    }
    cout<<endl;
    vector<Type> X = Solve(L, B, n);
    vector<Type>::iterator it;
    for(it = X.begin(); it != X.end(); it++)
        cout<<*it<<" ";
    cout<<endl;
}
 
int main()
{
    int n;
    cin>>n;
    memset(L, 0, sizeof(L));
    for(int i = 0; i < n; i++)
    {
        for(int j = 0; j < n; j++)
            cin>>A[i][j];
    }
    vector<Type> B;
    for(int i = 0; i < n; i++)
    {
        Type y;
        cin>>y;
        B.push_back(y);
    }
    Cholesky(A, L, n);
    Print(L, B, n);
    return 0;
}
 
/**data**
4
4 -2 4 2
-2 10 -2 -7
4 -2 8 4
2 -7 4 7
8 2 16 6
*/

reference

1.矩阵分解 LDL^T分解
2.矩阵分解 (乘法篇)