计算数学

正交变换与最小二乘

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');
		}
	}
}

Astrodynamics

计算数学

目录

正交变换与最小二乘

Cholesky分解