1. 全栈程序员必看首页

矩阵分解 Cholesky分解

全栈程序员-站长 • • 未分类 • 阅读 2

矩阵分解 Cholesky分解Cholesky 分解法又叫平方根法 是求解对称正定线性方程组最常用的方法之一 对于一般矩阵 为了消除 LU 分解的局限性和误差的过分积累 采用了选主元的方法 但对于对称正定矩阵而言 选主元是不必要的 nbsp 定理 若对称正定 则存在一个对角元为正数的下三角矩阵 使得成立 nbsp 假设现在要求解线性方程组 其中为对称正定矩阵 那么可通过下面步骤求解 nbsp 1 求的

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

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

 

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

 



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

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

 

（1）求矩阵分解 Cholesky分解的Cholesky分解，得到矩阵分解 Cholesky分解

（2）求解矩阵分解 Cholesky分解，得到矩阵分解 Cholesky分解

（3）求解矩阵分解 Cholesky分解，得到矩阵分解 Cholesky分解

 

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

 

       矩阵分解 Cholesky分解

 

通过矩阵分解 Cholesky分解比较两边的关系，首先由矩阵分解 Cholesky分解，再由

 

       矩阵分解 Cholesky分解

 

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

 

       矩阵分解 Cholesky分解

 

可以得到

 

       矩阵分解 Cholesky分解

 

进一步再由

 

                 矩阵分解 Cholesky分解

 

最终得到

 

       矩阵分解 Cholesky分解

 

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

 

代码：

[cpp]  view plain   copy

  1. #include 
      
  2. #include 
      
  3. #include 
      
  4. #include 
      
  5. #include 
      
  6.    
  7. using namespace std;  
  8. const int N = 1005;  
  9. typedef double Type;  
  10.    
  11. Type A[N][N], L[N][N];  
  12.    
  13. / 分解A得到A = L * L^T */  
  14. void Cholesky(Type A[][N], Type L[][N], int n)  
  15. {  
  16.     for(int k = 0; k < n; k++)  
  17.     {  
  18.         Type sum = 0;  
  19.         for(int i = 0; i < k; i++)  
  20.             sum += L[k][i] * L[k][i];  
  21.         sum = A[k][k] – sum;  
  22.         L[k][k] = sqrt(sum > 0 ? sum : 0);  
  23.         for(int i = k + 1; i < n; i++)  
  24.         {  
  25.             sum = 0;  
  26.             for(int j = 0; j < k; j++)  
  27.                 sum += L[i][j] * L[k][j];  
  28.             L[i][k] = (A[i][k] – sum) / L[k][k];  
  29.         }  
  30.         for(int j = 0; j < k; j++)  
  31.             L[j][k] = 0;  
  32.     }  
  33. }  
  34.    
  35. / 回带过程 */  
  36. vector

     Solve(Type L[][N], vector

     X, 
    int
     n)  

  37. {  
  38.     / LY = B  => Y */  
  39.     for(int k = 0; k < n; k++)  
  40.     {  
  41.         for(int i = 0; i < k; i++)  
  42.             X[k] -= X[i] * L[k][i];  
  43.         X[k] /= L[k][k];  
  44.     }  
  45.     / L^TX = Y => X */  
  46.     for(int k = n – 1; k >= 0; k–)  
  47.     {  
  48.         for(int i = k + 1; i < n; i++)  
  49.             X[k] -= X[i] * L[i][k];  
  50.         X[k] /= L[k][k];  
  51.     }  
  52.     return X;  
  53. }  
  54.    
  55. void Print(Type L[][N], const vector

     B, 
    int n)  
  56. {  
  57.     for(int i = 0; i < n; i++)  
  58.     {  
  59.         for(int j = 0; j < n; j++)  
  60.             cout<
    ” ”
    ;  
  61.         cout<
  62.     }  
  63.     cout<
  64.     vector

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

    ::iterator it;  
  66.     for(it = X.begin(); it != X.end(); it++)  
  67.         cout<<*it<<” “;  
  68.     cout<
  69. }  
  70.    
  71. int main()  
  72. {  
  73.     int n;  
  74.     cin>>n;  
  75.     memset(L, 0, sizeof(L));  
  76.     for(int i = 0; i < n; i++)  
  77.     {  
  78.         for(int j = 0; j < n; j++)  
  79.             cin>>A[i][j];  
  80.     }  
  81.     vector

     B;  
  82.     for(int i = 0; i < n; i++)  
  83.     {  
  84.         Type y;  
  85.         cin>>y;  
  86.         B.push_back(y);  
  87.     }  
  88.     Cholesky(A, L, n);  
  89.     Print(L, B, n);  
  90.     return 0;  
  91. }  
  92.    
  93. /data 
  94. 4 
  95. 4 -2 4 2 
  96. -2 10 -2 -7 
  97. 4 -2 8 4 
  98. 2 -7 4 7 
  99. 8 2 16 6 
  100. */  

 

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



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

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


版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请联系我们举报，一经查实，本站将立刻删除。

发布者：全栈程序员-站长，转载请注明出处：https://javaforall.net/214974.html原文链接：https://javaforall.net

赞 (0)
0 0

关于作者

全栈程序员-站长的头像

全栈程序员-站长

133.5K 文章
3 粉丝
本网站汇聚当前互联网主流语音，持续更新，欢迎关注公众号“全栈程序员社区”
上一篇 2026年3月18日 下午2:59
下一篇 2026年3月18日 下午3:00


相关推荐

发表回复

您的邮箱地址不会被公开。 必填项已用 * 标注

关注全栈程序员社区公众号