计算数学

content

正交变换与最小二乘

Cholesky分解

如果实矩阵A 对称正定，则存在一个实的下三角矩阵L ,使

{A = L L}^{T}

。若限定L的对角元素为正时，则这种分解是唯一的，这就是Cholesky分解，即

[\begin{matrix} a_{11} & a_{12} & \dots & a_{1 n} \\ a_{21} & a_{22} & \dots & a_{2 n} \\ ⋮ & ⋱ \\ a_{n 1} & a_{n 2} & \dots & a_{n n} \end{matrix}] = [\begin{matrix} l_{11} \\ l_{21} & l_{22} \\ \dots & ⋱ \\ l_{n 1} & l_{n 2} & \dots & l_{n n} \end{matrix}] [\begin{matrix} l_{11} & l_{21} & \dots & l_{n 1} \\ l_{22} & \dots & l_{n 2} \\ ⋱ \\ l_{n n} \end{matrix}]

计算时先令

L = 0

，然后计算
[STEP1]

l_{11} = \sqrt{a_{11}}

[STEP2]

l_{i 1} = \frac{a_{i 1}}{l_{11}}, i = 2, N

[STEP3] 对j = 2,N，做STEP4-STEP5。
[STEP4]

l_{j j} = \sqrt{a_{j j} - \sum_{k = 1}^{j - 1} l_{j k}^{2}}

[STEP5]

l_{i j} = \frac{(a_{i j} - \sum_{k = 1}^{j - 1} l_{i k} l_{j k})}{l_{j j}}, 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');
		}
	}
}