计算数学
正交变换与最小二乘
Cholesky分解
如果实矩阵A 对称正定， 则存在一个实的下三角矩阵L ,使$${\bf{A = LL}}^T$$ 。若限定L的对角元素为正时， 则这种分解是唯一的，这就是Cholesky分解，即 \[ \left[ {\begin{array}{*{20}c} {a_{11} } & {a_{12} } & \cdots & {a_{1n} } \\ {a_{21} } & {a_{22} } & \cdots & {a_{2n} } \\ \vdots & {} & \ddots & {} \\ {a_{n1} } & {a_{n2} } & \cdots & {a_{nn} } \\ \end{array}} \right] = \left[ {\begin{array}{*{20}c} {l_{11} } & {} & {} & {} \\ {l_{21} } & {l_{22} } & {} & {} \\ \cdots & {} & \ddots & {} \\ {l_{n1} } & {l_{n2} } & \cdots & {l_{nn} } \\ \end{array}} \right]\left[ {\begin{array}{*{20}c} {l_{11} } & {l_{21} } & \cdots & {l_{n1} } \\ {} & {l_{22} } & \cdots & {l_{n2} } \\ {} & {} & \ddots & {} \\ {} & {} & {} & {l_{nn} } \\ \end{array}} \right] \] 计算时先令 $${\bf{L=0}}$$，然后计算[STEP1] $$l_{11} = \sqrt {a_{11} }$$ [STEP2] $$l_{i1} = \frac{{a_{i1} }}{{l_{11} }},i = 2,N$$ [STEP3] 对j = 2,N，做STEP4-STEP5。
[STEP4] $$l_{jj} = \sqrt {a_{jj} - \sum\limits_{k = 1}^{j - 1} {l_{jk}^2 } }$$ [STEP5] $$l_{ij} = \frac{{\left( {a_{ij} - \sum\limits_{k = 1}^{j - 1} {l_{ik} l_{jk} } } \right)}}{{l_{jj} }},\quad i = j + 1,N$$ 由以上几步就已经算出了L 矩阵，Cholesky分解经常用于对称正定方程的求解， 许多软件包如MATLAB等矩阵分析软件都提供了Cholesky分解的函数。
下面给出cholesky分解的Java版本代码。
package cholesky;
/**
* @author song Yezhi
* @version 1.0
* @since 2019.05.01
*/
class cholesky
{
/**
* cholesky分解
* @param A --对称正定矩阵
* @return L--下三角矩阵
*/
static double[][] chol(double[][] A)
{
int N = A.length;
double L[][] = new double[N][N];
int i, j;
for (i = 0; i < N; i++)
for (j = 0; j < N; j++)
L[i][j] = 0.0;
L[0][0] = Math.sqrt(A[0][0]);
for (i = 1; i < N; i++)
L[i][0] = A[i][0] / L[0][0];
double sum;
int k;
for (j = 1; j < N; j++)
{
sum = 0.0;
for (k = 0; k <= j - 1; k++)
sum = sum + L[j][k] * L[j][k];
L[j][j] = Math.sqrt(A[j][j] - sum);
for (i = j + 1; i < N; i++)
{
sum = 0;
for (k = 0; k <= j - 1; k++)
sum = sum + L[i][k] * L[j][k];
L[i][j] = (A[i][j] - sum) / L[j][j];
}
}
return L;
}
}
package cholesky;
import java.util.Scanner;
import java.io.FileInputStream;
import java.io.FileNotFoundException;
import java.io.PrintWriter;
/**
* 计算cholesky分解的算例
* @author song
*/
public class pro
{
public static void main(String[] args)
{
java.io.File file = new java.io.File("report.txt");
PrintWriter fout = null;
try
{
fout = new PrintWriter(file);
} catch (FileNotFoundException e)
{
System.out.println("Erro opening the file");
}
double[][] A; // = new double [][];
A = readMat();
System.out.println("A=");
fout.println("A=");
printMat(A, fout);
// print mat to the screen
double L[][];
L = cholesky.chol(A);
System.out.println("L=");
fout.println("L=");
printMat(L, fout);
fout.close();
}
/**
* 从文件读取矩阵
* @return 读取的矩阵
*/
static double[][] readMat()
{
Scanner fin = null;
try
{
fin = new Scanner(new FileInputStream("fin.txt"));
} catch (FileNotFoundException e)
{
System.exit(0);
}
int M, N;
M = fin.nextInt();
N = fin.nextInt();
double[][] A = new double[M][N];
int i, j;
for (i = 0; i < M; i++)
for (j = 0; j < N; j++)
A[i][j] = fin.nextDouble();
fin.close();
return A;
}
/**
* 输出矩阵到屏幕和文件中
* @param A
* @param fout
*/
static void printMat(double[][] A, PrintWriter fout)
{
int M = A.length;
int N = A[0].length;
int i, j;
for (i = 0; i < M; i++)
{
for (j = 0; j < N; j++)
{
System.out.printf("%18.8f", A[i][j]);
fout.printf("%18.8f", A[i][j]);
}
System.out.print('\n');
fout.print('\n');
}
}
}