Cholesky 变形LDL分解

2024-07-03 16:46| 来源: 网络整理| 查看: 265

概念

LDL分解是经典Cholesky分解的一个变形： $\mathbf {A} =\mathbf {LDL} ^{*}=\mathbf {LD} ^{\frac {1}{2}}(\mathbf {D} ^{\frac {1} {2}})^{*}\mathbf {L} ^{*}=\mathbf {LD} ^{\frac {1}{2}}(\mathbf {LD} ^{\frac {1}{2}})^{*}}{\displaystyle \mathbf {A} =\mathbf {LDL} ^{*}=\mathbf {LD} ^{\frac {1}{2}}(\mathbf {D} ^{\frac {1}{2}})^{*}\mathbf {L} ^{*}=\mathbf {LD} ^{\frac {1}{2}}(\mathbf {LD} ^{\frac {1}{2}})^{*}$

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

$\mathbf {D} =\mathbf {S} ^{2}} \\ {\displaystyle \mathbf {L} =\mathbf {L} ^{Cholesky}\mathbf {S} ^{-1}}{\displaystyle \mathbf {L} =\mathbf {L} ^{Cholesky}\mathbf {S} ^{-1}$

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

因此人们更倾向于使用LDL分解。对于实数矩阵，该种分解的形式可被改写成 $\mathbf {A} =\mathbf {LDL} ^{\mathbf {T} }$

计算

$d_{i,i}=x_{i,i}-\sum_{k=0}^{i-1} x_{i,k}^2*d_{k,k} (1)\\ x_{i,j}=\frac {x_{i,j}-\sum_{k=0}^{i-1} x_{i,k}* x_{j,k}*d_{k,k}} { d_{i,i}} (2)$

减少乘法运算，引入辅助量 $u_{ik}=x_{i,k}*d_{k,k}$ ,则公式为： $d_{i,i}=x_{i,i}-\sum_{k=0}^{i-1} x_{i,k}* u_{i,k} (1)\\ u_{i,j}=x_{i,j}-\sum_{k=0}^{i-1} u_{i,k}* x_{j,k} (2)\\ l{i,j}= \frac {u_{i,j}}{d_{i,i}}$

代码 #include #include #include #include #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

【本文地址】

公司简介

联系我们