矩阵分解 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)求解,得到
代码:
- #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], 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<Type> Solve(Type L[][N], Type D[][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];
- }
- /** 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<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, D, 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, 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
http://blog.csdn.net/zhouliyang1990/article/details/21952485
矩阵分解 LDL^T分解相关推荐
- 解方程AX=b与矩阵分解:奇异值分解(SVD分解) 特征值分解 QR分解 三角分解 LLT分解
文章目录 1. 前言 2. LU三角分解 3. Cholesky分解 - LDLT分解 4. Cholesky分解 - LLT分解 5. QR分解 6. 奇异值分解 7. 特征值分解 1. 前言 本博 ...
- 视觉SLAM中的数学——解方程AX=b与矩阵分解:奇异值分解(SVD分解) 特征值分解 QR分解 三角分解 LLT分解
前言 本博客主要介绍在SLAM问题中常常出现的一些线性代数相关的知识,重点是如何采用矩阵分解的方法,求解线性方程组AX=B.主要参考了<计算机视觉--算法与应用>附录A以及Eigen库的方 ...
- 三阶矩阵的lu分解详细步骤_数学 - 线性代数导论 - #4 矩阵分解之LU分解的意义、步骤和成立条件...
线性代数导论 - #4 矩阵分解之LU分解的意义.步骤和成立条件 目前我们用于解线性方程组的方法依然是Gauss消元法.在Gauss消元法中,我们将右侧向量b与A写在一起作为一个增广矩阵进行同步的操作 ...
- 7.4.3 矩阵极分解和平方根分解
7.4.3 矩阵极分解和平方根分解 当矩阵 AAA 是方阵时 A=UΣVT=UVTVΣVT=(UVT)(VΣVT)=QSQ=UVT是正交矩阵,S=VΣVT是对称半正定矩阵,即对任意向量x,有xTSx≥ ...
- 机器学习(十一)——机器学习中的矩阵方法(1)LU分解、QR分解
http://antkillerfarm.github.io/ 因子分析的EM估计(续) 去掉和各参数无关的部分后,可得: ∑i=1mE[logp(x(i)|z(i);μ,Λ,Ψ)]=∑i=1mE[1 ...
- 将投影矩阵P利用QR分解分解出摄像机内外参数(Opencv)
将投影矩阵P利用QR分解分解出摄像机内外参数(Opencv) /*************************************************************** ...
- 几种矩阵分解算法: LU分解,Cholesky分解,QR分解,SVD分解,Jordan分解
目录 1.LU分解 2. LDLT分解法 3. Cholesky分解的形式 4. QR分解 5.SVD分解 5.1 SVD与广义逆矩阵 6. Jordan 分解 参考文章: ---------我只是搬 ...
- 矩阵的LU分解,LU分解的推广,LU分解有什么意义,为什么要用LU分解。
一点点数学!开干! 参考书籍:<矩阵分析与计算>李继根 张新发编著 矩阵的LU分解: LU分解定理:如果n阶方阵A的各阶顺序主子式≠0(K=1.2.3,-,n),即A的各阶顺序主子式矩阵都 ...
- 矩阵的三角分解(LU分解)
矩阵的三角分解将矩阵分解成一个下三角矩阵L和一个上三角矩阵U的乘积. 定义:如果n阶矩阵A能够分解成一个下三角矩阵L和一个上三角矩阵U的乘积,则称这种分解为三角分解或LU分解,如果n阶矩阵A能够分解为 ...
最新文章
- 链表问题11——两个单链表相交的系列问题(一):找到有环链表的环入口节点
- 【AIX 命令学习】加载与卸载文件系统!
- 【数字信号处理】序列傅里叶变换 ( 狄义赫利条件 | 序列傅里叶变换定义 )
- HTML中通过CSS方式隐藏元素
- Drools学习 入门实例
- sql中索引不会被用到的几种情况
- Kafka核心源码解析 - KafkaController源码解析
- h5页面调用相机功能
- Sudo bug 可导致非权限 Linux 和 MacOS 用户以根身份运行命令
- VIM之taglist
- Java设计模式之观察者模式(发布-订阅模式)
- 《unity2021》如何改成中文
- easyui图标对照
- 基于Android的健康打卡系统,基于Android平台的个人健康管理系统
- 服务器数据恢复的两种方法
- ietester测试本地html,网站浏览器兼容测试软件–IETester
- 计算机本地网络给手机使用吗,手机网络也能共享给电脑(台式or笔记本)使用吗?...
- VB完全控制IE浏览器,操作ie对象,响应ie事件
- 基金从业考试一般要准备多长时间?
- kinhdown稳定版无法连接服务器,KinhDown稳定版