矩阵分解 LDL^T分解矩阵分解

文章图片
分解

实际问题中，当求解方程组的系数矩阵是对称矩阵时，则用下面介绍的
文章图片
分解法可以简化程序设计并减少计算量。
从定理可知，当矩阵A的各阶顺序主子式不为零时，A有唯一的Doolittle分解A= LU。矩阵U的对角线元素Uii 不等于0，将矩阵U的每行依次提出，

文章图片

下面将U分解为
文章图片

文章图片

定理：若对称矩阵A的各阶顺序主子式不为零时，则A可以唯一分解为A=
文章图片
，这里

文章图片

L^T为L的转置矩阵。
当A有
文章图片
分解时，利用矩阵运算法则及相等原理易得计算ljk和dk的公式为

文章图片

将对称正定矩阵
文章图片
通过分解成
文章图片
，其中
文章图片
是单位下三角矩阵，
文章图片
是对角均为正数的对角矩阵。把这
一分解叫做
文章图片
分解，是Cholesky分解的变形。对应两边的元素，很容易得到

文章图片

由此可以确定计算
文章图片
和
文章图片
的公式如下

文章图片

在实际计算时，是将
文章图片
的严格下三角元素存储在
文章图片
的对应位置上，而将
文章图片
的对角元存储在
文章图片
的对应的对角位置上。
类似地求解线性方程组
文章图片
的解步骤如下

（1）对矩阵
文章图片
进行分解得到
文章图片

（2）求解
文章图片
，得到
文章图片

（3）求解
文章图片
，得到
文章图片

代码：
[cpp]view plaincopy

#include
#include
#include
#include
#include
using namespace std;
const int N = 1005;
typedef double Type;
Type A[N][N], L[N][N], D[N][N];
/** 分解A得到A = LDL^T */
void Cholesky(Type A[][N], Type L[][N], Type D[][N], int n)
{
for(int k = 0; k < n; k++)
{
for(int i = 0; i < k; i++)
A[k][k] -= A[i][i] * A[k][i] * A[k][i];
for(int j = k + 1; j < n; j++)
{
for(int i = 0; i < k; i++)
A[j][k] -= A[j][i] * A[i][i] * A[k][i];
A[j][k] /= A[k][k];
}
}
memset(L, 0, sizeof(L));
memset(D, 0, sizeof(D));
for(int i = 0; i < n; i++)
{
D[i][i] = A[i][i];
L[i][i] = 1;
}
for(int i = 0; i < n; i++)
{
for(int j = 0; j < i; j++)
L[i][j] = A[i][j];
}
}
void Transposition(Type L[][N], int n)
{
for(int i = 0; i < n; i++)
{
for(int j = 0; j < i; j++)
swap(L[i][j], L[j][i]);
}
}
void Multi(Type A[][N], Type B[][N], int n)
{
Type **C = new Type*[n];
for(int i = 0; i < n; i++)
C[i] = new Type[n];
for(int i = 0; i < n; i++)
{
for(int j = 0; j < n; j++)
{
C[i][j] = 0;
for(int k = 0; k < n; k++)
C[i][j] += A[i][k] * B[k][j];
}
}
for(int i = 0; i < n; i++)
{
for(int j = 0; j < n; j++)
B[i][j] = C[i][j];
}
for(int i = 0; i < n; i++)
{
delete[] C[i];
C[i] = NULL;
}
delete C;
C = NULL;
}
/** 回带过程 */
vector Solve(Type L[][N], Type D[][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];
}
/** DL^TX = Y => X */
Transposition(L, n);
Multi(D, L, n);
for(int k = n - 1; k >= 0; k--)
{
for(int i = k + 1; i < n; i++)
X[k] -= X[i] * L[k][i];
X[k] /= L[k][k];
}
return X;
}
void Print(Type L[][N], Type D[][N], const vector B, int n)
{
for(int i = 0; i < n; i++)
{
for(int j = 0; j < n; j++)
cout<
cout<
}
cout<
vector X = Solve(L, D, B, n);
vector::iterator it;
for(it = X.begin(); it != X.end(); it++)
cout<<*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, D, n);
Print(L, D, 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
*/

参考资料：http://class.htu.cn/nla/cha1/sect3.htm

转自：http://blog.csdn.net/acdreamers/article/details/44656847
【矩阵分解 LDL^T分解】 http://blog.csdn.net/zhouliyang1990/article/details/21952485