cholesky分解java代码_Cholesky 变形LDL分解

最新推荐文章于 2024-12-17 20:00:00 发布

weixin_39633493

最新推荐文章于 2024-12-17 20:00:00 发布

阅读量280

点赞数

文章标签： cholesky分解java代码

本文链接：https://blog.youkuaiyun.com/weixin_39633493/article/details/114490889

版权

概念

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

$math?formula=%7B%5Cdisplaystyle%20%5Cmathbf%20%7BA%7D%20%3D%5Cmathbf%20%7BLDL%7D%20%5E%7B*%7D%3D%5Cmathbf%20%7BLD%7D%20%5E%7B%5Cfrac%20%7B1%7D%7B2%7D%7D(%5Cmathbf%20%7BD%7D%20%5E%7B%5Cfrac%20%7B1%7D%20%7B2%7D%7D)%5E%7B*%7D%5Cmathbf%20%7BL%7D%20%5E%7B*%7D%3D%5Cmathbf%20%7BLD%7D%20%5E%7B%5Cfrac%20%7B1%7D%7B2%7D%7D(%5Cmathbf%20%7BLD%7D%20%5E%7B%5Cfrac%20%7B1%7D%7B2%7D%7D)%5E%7B*%7D%7D%7B%5Cdisplaystyle%20%5Cmathbf%20%7BA%7D%20%3D%5Cmathbf%20%7BLDL%7D%20%5E%7B*%7D%3D%5Cmathbf%20%7BLD%7D%20%5E%7B%5Cfrac%20%7B1%7D%7B2%7D%7D(%5Cmathbf%20%7BD%7D%20%5E%7B%5Cfrac%20%7B1%7D%7B2%7D%7D)%5E%7B*%7D%5Cmathbf%20%7BL%7D%20%5E%7B*%7D%3D%5Cmathbf%20%7BLD%7D%20%5E%7B%5Cfrac%20%7B1%7D%7B2%7D%7D(%5Cmathbf%20%7BLD%7D%20%5E%7B%5Cfrac%20%7B1%7D%7B2%7D%7D)%5E%7B*%7D%7D$

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

D为对角矩阵。

转换

math?formula=%7B%5Cdisplaystyle%20%5Cmathbf%20%7BD%7D%20%3D%5Cmathbf%20%7BS%7D%20%5E%7B2%7D%7D%20%5C%5C%20%7B%5Cdisplaystyle%20%5Cmathbf%20%7BL%7D%20%3D%5Cmathbf%20%7BL%7D%20%5E%7BCholesky%7D%5Cmathbf%20%7BS%7D%20%5E%7B-1%7D%7D%7B%5Cdisplaystyle%20%5Cmathbf%20%7BL%7D%20%3D%5Cmathbf%20%7BL%7D%20%5E%7BCholesky%7D%5Cmathbf%20%7BS%7D%20%5E%7B-1%7D%7D

特征

LDL变形如果得以有效运行，构造及使用时所需求的空间及计算的复杂性与经典Cholesky分解是相同的，但是可避免提取平方根。

某些不存在Cholesky分解的不定矩阵，也可以运行LDL分解，此时矩阵

math?formula=%7B%5Cdisplaystyle%20%5Cmathbf%20%7BD%7D%20%7D%5Cmathbf%20%7BD%7D 中会出现负数元素。

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

math?formula=%7B%5Cdisplaystyle%20%5Cmathbf%20%7BA%7D%20%3D%5Cmathbf%20%7BLDL%7D%20%5E%7B%5Cmathbf%20%7BT%7D%20%7D%7D

计算

$math?formula=d_%7Bi%2Ci%7D%3Dx_%7Bi%2Ci%7D-%5Csum_%7Bk%3D0%7D%5E%7Bi-1%7D%20x_%7Bi%2Ck%7D%5E2*d_%7Bk%2Ck%7D%20(1)%5C%5C%20x_%7Bi%2Cj%7D%3D%5Cfrac%20%7Bx_%7Bi%2Cj%7D-%5Csum_%7Bk%3D0%7D%5E%7Bi-1%7D%20x_%7Bi%2Ck%7D*%20x_%7Bj%2Ck%7D*d_%7Bk%2Ck%7D%7D%20%7B%20d_%7Bi%2Ci%7D%7D%20(2)$

减少乘法运算，引入辅助量

math?formula=u_%7Bik%7D%3Dx_%7Bi%2Ck%7D*d_%7Bk%2Ck%7D ,则公式为：

$math?formula=d_%7Bi%2Ci%7D%3Dx_%7Bi%2Ci%7D-%5Csum_%7Bk%3D0%7D%5E%7Bi-1%7D%20x_%7Bi%2Ck%7D*%20u_%7Bi%2Ck%7D%20(1)%5C%5C%20u_%7Bi%2Cj%7D%3Dx_%7Bi%2Cj%7D-%5Csum_%7Bk%3D0%7D%5E%7Bi-1%7D%20u_%7Bi%2Ck%7D*%20x_%7Bj%2Ck%7D%20(2)%5C%5C%20l%7Bi%2Cj%7D%3D%20%5Cfrac%20%7Bu_%7Bi%2Cj%7D%7D%7Bd_%7Bi%2Ci%7D%7D$

代码

#include

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 Solve(Type L[][N], vector 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 B, int n)

{

for(int i = 0; i < n; i++)

{

for(int j = 0; j < n; j++)

cout<

}

cout<

vector X = Solve(L, B, n);

vector::iterator it;

for(it = X.begin(); it != X.end(); it++)

cout<

}

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 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 -2 4 2

-2 10 -2 -7

4 -2 8 4

2 -7 4 7

8 2 16 6

reference