数学之美：cholesky矩阵分解

Cholesky分解法又叫平方根法，是求解对称正定线性方程组最常用的方法之一。对于一般矩阵，为了消除LU分

解的局限性和误差的过分积累，采用了选主元的方法，但对于对称正定矩阵而言，选主元是不必要的。

定理：若对称正定，则存在一个对角元为正数的下三角矩阵，使得成立。

Cholesky分解的条件(这里针对复数矩阵)

一、Hermitian matrix（埃尔米特矩阵）：矩阵中的元素共轭对称（复数域的定义，类比于实数对称矩阵）。

Hermitian意味着对于任意向量x和y，(x*)Ay共轭相等

二、Positive-definite：正定（矩阵域，类比于正实数的一种定义）。正定矩阵A意味着，对于任何向量x，(x^T)Ax总是大于零(复数域是(x*)Ax>0)

Cholesky分解的形式

可记作A = L L*。其中L是下三角矩阵。L*是L的共轭转置矩阵。

可以证明，只要A满足以上两个条件，L是唯一确定的，而且L的对角元素肯定是正数。反过来也对，即存在L把A分解的话，A满足以上两个条件。

如果A是半正定的（semi-definite），也可以分解，不过这时候L就不唯一了。

特别的，如果A是实数对称矩阵，那么L的元素肯定也是实数。

另外，满足以上两个条件意味着A矩阵的特征值都为正实数，因为Ax = lamda * x,

(x*)Ax = lamda * (x*)x > 0, lamda > 0

假设现在要求解线性方程组，其中为对称正定矩阵，那么可通过下面步骤求解

（1）求的Cholesky分解，得到

（2）求解，得到

（3）求解，得到

现在的关键问题是对进行Cholesky分解。假设

通过比较两边的关系，首先由，再由

这样便得到了矩阵的第一列元素，假定已经算出了的前列元素，通过

可以得到

进一步再由

最终得到

这样便通过的前列求出了第列，一直递推下去即可求出，这种方法称为平方根法。

代码：

#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
*/

用上述的方法需要进行开方，这有可能损失精度和增加运算量，为了避免开方，Cholesky分解有个改进的版本分解。

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

转自：http://blog.csdn.net/acdreamers/article/details/44656847

http://blog.csdn.net/zhouliyang1990/article/details/21952485

数学之美：cholesky矩阵分解相关推荐

TensorFlow通过Cholesky矩阵分解实现线性回归
这里将用TensorFlow实现矩阵分解,对于一般的线性回归算法,求解Ax=b,则x=(A'A)^(-1)A'b,然而该方法在大部分情况下是低效率的,特别是当矩阵非常大时效率更低.另一种方式则是通过矩 ...
三阶矩阵的lu分解详细步骤_数学 - 线性代数导论 - #4 矩阵分解之LU分解的意义、步骤和成立条件...
线性代数导论 - #4 矩阵分解之LU分解的意义.步骤和成立条件目前我们用于解线性方程组的方法依然是Gauss消元法.在Gauss消元法中,我们将右侧向量b与A写在一起作为一个增广矩阵进行同步的操作 ...
[Math]常见矩阵分解及复杂度 Cholesky QR LU
本篇尝试介绍几个常见的矩阵分解及其数值方法,探究其算法复杂度及不稳定性来源,能够让我们更好的理解工具,使用工具.常见的矩阵分解可以在维基中找到文章目录 Cholesky 定义数值方法数值方法的数 ...
Cholesky和LU矩阵分解
1.Cholesky分解在线性代数中,矩阵分解是将矩阵分解为矩阵的乘积.有许多不同的矩阵分解.其中之一就是Cholesky分解. Cholesky分解是将Hermitian正定矩阵分解为下三角矩阵及 ...
视觉SLAM中的数学——解方程AX=b与矩阵分解：奇异值分解（SVD分解）特征值分解 QR分解三角分解 LLT分解
前言本博客主要介绍在SLAM问题中常常出现的一些线性代数相关的知识,重点是如何采用矩阵分解的方法,求解线性方程组AX=B.主要参考了<计算机视觉--算法与应用>附录A以及Eigen库的方 ...
【高等工程数学】南理工研究生课程突击笔记5 矩阵分解与广义逆矩阵
矩阵分解与广义逆矩阵文章目录矩阵分解与广义逆矩阵写在前面一.矩阵分解三角分解(Doolittle分解) 二.满秩分解 Hermite标准型广义逆矩阵A+ 解总结写在前面第三章主要内容 ...
矩阵分解——8.2 乔里斯基（Cholesky）分解
8.2 乔里斯基(Cholesky)分解
推荐算法矩阵分解实战——keras算法练习
当今这个信息爆炸的社会,每个人都会面对无数的商品,无数的选择.而推荐算法的目的帮助大家解决选择困难症的问题,在大千世界中推荐专属于你的商品. 推荐系统算法简介这里简单介绍下推荐系统中最为主要的协同过 ...
Eigen密集矩阵求解 1 - 线性代数及矩阵分解
简介这里介绍线性系统的解析,如何进行各种分解计算,如LU,QR,SVD,特征值分解等. 简单线性求解在一个线性系统,常如下表示,其中A,b分别是一个矩阵,需要求x: Ax=bAx \:= \: b ...

数学之美：cholesky矩阵分解

数学之美：cholesky矩阵分解相关推荐

最新文章

热门文章