Cholesky分解法

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

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

 

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

 

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

 

（1）求的Cholesky分解，得到

（2）求解，得到

（3）求解，得到

 

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

 

       

 

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

 

       

 

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

 

       

 

可以得到

 

       

 

进一步再由

 

                 

 

最终得到

 

       

 

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

 

代码：

[cpp] view plain copy

1. #include <iostream>  
2. #include <string.h>  
3. #include <stdio.h>  
4. #include <vector>  
5. #include <math.h>  
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<Type> Solve(Type L[][N], vector<Type> 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<Type> B, int n)  
56. {  
57.     for(int i = 0; i < n; i++)  
58.     {  
59.         for(int j = 0; j < n; j++)  
60.             cout<<L[i][j]<<" ";  
61.         cout<<endl;  
62.     }  
63.     cout<<endl;  
64.     vector<Type> X = Solve(L, B, n);  
65.     vector<Type>::iterator it;  
66.     for(it = X.begin(); it != X.end(); it++)  
67.         cout<<*it<<" ";  
68.     cout<<endl;  
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<Type> 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. */  

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

 

将对称正定矩阵通过分解成，其中是单位下三角矩阵，是对角均为正数的对角矩阵。把这

一分解叫做分解，是Cholesky分解的变形。对应两边的元素，很容易得到

 

          

 

由此可以确定计算和的公式如下

 

          

 

在实际计算时，是将的严格下三角元素存储在的对应位置上，而将的对角元存储在的对应的对角位置上。

类似地求解线性方程组的解步骤如下

 

（1）对矩阵进行分解得到

（2）求解，得到

（3）求解，得到

 

代码：

[cpp] view plain copy

1. #include <iostream>  
2. #include <string.h>  
3. #include <stdio.h>  
4. #include <vector>  
5. #include <math.h>  
6.    
7. using namespace std;  
8. const int N = 1005;  
9. typedef double Type;  
10.    
11. Type A[N][N], L[N][N], D[N][N];  
12.    
13. /** 分解A得到A = LDL^T */  
14. void Cholesky(Type A[][N], Type L[][N], Type D[][N], int n)  
15. {  
16.     for(int k = 0; k < n; k++)  
17.     {  
18.         for(int i = 0; i < k; i++)  
19.             A[k][k] -= A[i][i] * A[k][i] * A[k][i];  
20.         for(int j = k + 1; j < n; j++)  
21.         {  
22.             for(int i = 0; i < k; i++)  
23.                 A[j][k] -= A[j][i] * A[i][i] * A[k][i];  
24.             A[j][k] /= A[k][k];  
25.         }  
26.     }  
27.     memset(L, 0, sizeof(L));  
28.     memset(D, 0, sizeof(D));  
29.     for(int i = 0; i < n; i++)  
30.     {  
31.         D[i][i] = A[i][i];  
32.         L[i][i] = 1;  
33.     }  
34.     for(int i = 0; i < n; i++)  
35.     {  
36.         for(int j = 0; j < i; j++)  
37.             L[i][j] = A[i][j];  
38.     }  
39. }  
40.    
41. void Transposition(Type L[][N], int n)  
42. {  
43.     for(int i = 0; i < n; i++)  
44.     {  
45.         for(int j = 0; j < i; j++)  
46.             swap(L[i][j], L[j][i]);  
47.     }  
48. }  
49.    
50. void Multi(Type A[][N], Type B[][N], int n)  
51. {  
52.     Type **C = new Type*[n];  
53.     for(int i = 0; i < n; i++)  
54.         C[i] = new Type[n];  
55.     for(int i = 0; i < n; i++)  
56.     {  
57.         for(int j = 0; j < n; j++)  
58.         {  
59.             C[i][j] = 0;  
60.             for(int k = 0; k < n; k++)  
61.                 C[i][j] += A[i][k] * B[k][j];  
62.         }  
63.     }  
64.     for(int i = 0; i < n; i++)  
65.     {  
66.         for(int j = 0; j < n; j++)  
67.             B[i][j] = C[i][j];  
68.     }  
69.     for(int i = 0; i < n; i++)  
70.     {  
71.         delete[] C[i];  
72.         C[i] = NULL;  
73.     }  
74.     delete C;  
75.     C = NULL;  
76. }  
77.    
78. /** 回带过程 */  
79. vector<Type> Solve(Type L[][N], Type D[][N], vector<Type> X, int n)  
80. {  
81.     /** LY = B  => Y */  
82.     for(int k = 0; k < n; k++)  
83.     {  
84.         for(int i = 0; i < k; i++)  
85.             X[k] -= X[i] * L[k][i];  
86.         X[k] /= L[k][k];  
87.     }  
88.     /** DL^TX = Y => X */  
89.     Transposition(L, n);  
90.     Multi(D, L, n);  
91.     for(int k = n - 1; k >= 0; k--)  
92.     {  
93.         for(int i = k + 1; i < n; i++)  
94.             X[k] -= X[i] * L[k][i];  
95.         X[k] /= L[k][k];  
96.     }  
97.     return X;  
98. }  
99.    
100. void Print(Type L[][N], Type D[][N], const vector<Type> B, int n)  
101. {  
102.     for(int i = 0; i < n; i++)  
103.     {  
104.         for(int j = 0; j < n; j++)  
105.             cout<<L[i][j]<<" ";  
106.         cout<<endl;  
107.     }  
108.     cout<<endl;  
109.     vector<Type> X = Solve(L, D, B, n);  
110.     vector<Type>::iterator it;  
111.     for(it = X.begin(); it != X.end(); it++)  
112.         cout<<*it<<" ";  
113.     cout<<endl;  
114. }  
115.    
116. int main()  
117. {  
118.     int n;  
119.     cin>>n;  
120.     memset(L, 0, sizeof(L));  
121.     for(int i = 0; i < n; i++)  
122.     {  
123.         for(int j = 0; j < n; j++)  
124.             cin>>A[i][j];  
125.     }  
126.     vector<Type> B;  
127.     for(int i = 0; i < n; i++)  
128.     {  
129.         Type y;  
130.         cin>>y;  
131.         B.push_back(y);  
132.     }  
133.     Cholesky(A, L, D, n);  
134.     Print(L, D, B, n);  
135.     return 0;  
136. }  
137.    
138. /**data** 
139. 4 
140. 4 -2 4 2 
141. -2 10 -2 -7 
142. 4 -2 8 4 
143. 2 -7 4 7 
144. 8 2 16 6 
145. */  

 

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