song.yz@foxmail.com wechat: math-box

计算数学



LDL'分解

对于对系数是对称正定矩阵的方程,虽然也可以用普通的高斯消去法法解算, 但是没有充分利用其矩阵的特点。处理这一问题较为常用的方法是Cholesky分解法或称为LDL'分解, 在Cholesky分解过程中需要进行多次开平方运算,故而该方法又称平方根法。

可以对Cholesky分解进行修正可以得到不开平方的Cholesky分解法即LDL'。 其中L为单位下三角阵,D为对角矩阵,形如

$$ {\bf{A}} = \left[ {\begin{array}{*{20}{c}} 1&{}&{}&{}\\ {{l_{21}}}&1&{}&{}\\ \vdots & \vdots & \ddots &{}\\ {{l_{n1}}}&{{l_{n2}}}& \cdots &1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{d_1}}&{}&{}&{}\\ {}&{{d_2}}&{}&{}\\ {}&{}& \ddots &{}\\ {}&{}&{}&{{d_n}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} 1&{{l_{21}}}& \cdots &{{l_{n1}}}\\ {}&1& \cdots &{{l_{n2}}}\\ {}&{}& \ddots & \vdots \\ {}&{}&{}&1 \end{array}} \right] $$

这里给出修正的不开平方Cholesky分解算法。

$$ \begin{array}{l} {d_1} = {a_{11}}\\ \left. \begin{array}{l} {g_{ij}} = {a_{ij}} - \sum\limits_{k = 1}^{j - 1} {{g_{ik}}{l_{jk}}} \left( {j = 1, \cdots ,i - 1} \right)\\ {l_{ij}} = \frac{{{g_{ij}}}}{{{d_j}}}\left( {j = 1, \cdots ,i - 1} \right)\\ {d_i} = {a_{ii}} - \sum\limits_{k = 1}^{i - 1} {{g_{ik}}{l_{ik}}} \end{array} \right\}i = 2, \cdots ,n \end{array} $$

以下为go语言版本算法,整个过程没有开平方运算。


// matlib.go
// 导出函数首字母需要大写
package main

import (
	"fmt"
)


func Vec(N int) (V []float64) {
	/*------------------------------------------------------
	!  Author  : Song Yezhi      
	!  verison : 2020-5-24 17:50
	!
	!
	!  -----------------------------------------------------
	!  Input  Parameters :
	!
	!  Output Parameters :
	!
	------------------------------------------------------*/
	V = make([]float64, N)
	for i := 0; i < N; i++ {
		V[i] = 0.0
	}
	return
}

func Mat(M int, N int) (A [][]float64) {
	/*------------------------------------------------------
	  !  Author  : Song Yezhi      
	  !  verison : 2020-5-23 23:43
	  !
	  !  creat a two dimention slices (matrix)
	  !  -----------------------------------------------------
	  !  Input  Parameters :
	  !
	  !  Output Parameters :
	  !
	  ------------------------------------------------------*/
	A = make([][]float64, M)
	for i := 0; i < M; i++ {
		A[i] = make([]float64, N)
	}

	for i := 0; i < M; i++ {
		for j := 0; j < N; j++ {
			A[i][j] = 0.0
		}
	}

	return
}

func MatDim(A [][]float64) (M int, N int) {
	/*------------------------------------------------------
	!  Author  : Song Yezhi      
	!  verison : 2020-5-24 12:45
	------------------------------------------------------*/
	M = len(A)
	N = len(A[M-1])
	return
}

func VecOutput(V []float64) {
	N := len(V)
	for i := 0; i < N; i++ {
		fmt.Printf("%18.6f\n", V[i])
	}
}

func VecDot(x1 []float64, x2 []float64) (dot float64) {
	/*------------------------------------------------------
	!  Author  : Song Yezhi      
	!  verison : 2020-5-19 22:53
	!
	!     dot product
	------------------------------------------------------*/
	N := len(x1)
	var i int
	dot = 0.0
	for i = 0; i < N; i++ {
		dot = dot + x1[i]*x2[i]
	}
	return
}

func MatVec(A [][]float64, x []float64) (b []float64) {
	/*------------------------------------------------------
	!  Author  : Song Yezhi      
	!  verison : 2020-5-19 22:25
	!      b = Ax
	!
	!  -----------------------------------------------------
	!  Input  Parameters :
	!     A[M][N] ----MATRIX
	!  Output Parameters :
	!     b[M]-----vector
	------------------------------------------------------*/
	M := len(A)
	N := len(A[0])
	//the size of mat

	//b = make([]float64, M)
  
  b= Vec(M)

	var sum float64

	for i := 0; i < M; i++ {
		sum = 0.0
		for j := 0; j < N; j++ {
			sum = sum + A[i][j]*x[j]
		}
		b[i] = sum
	}

	return
}

func MatMul(A [][]float64, B [][]float64) (C [][]float64) {
	/*------------------------------------------------------
	!  Author  : Song Yezhi      
	!  verison : 2020-5-23 23:27
	!    C = A*B
	!
	!  -----------------------------------------------------
	!  Input  Parameters :
	!     A[M][N]
	!     B[N][P]
	!  Output Parameters :
	!     C[M][P]
	------------------------------------------------------*/
	M := len(A)
	N := len(A[M-1])
	P := len(B[N-1])

	var i, j, k int

	C = Mat(M, P)
	var sum float64

	for i = 0; i < M; i++ {
		for j = 0; j < P; j++ {
			sum = 0.0
			for k = 0; k < N; k++ {
				sum = sum + A[i][k]*B[k][j]
				C[i][j] = sum
			}
		}
	}

	return
}



func MatOutput(A [][]float64) {
	/*------------------------------------------------------
	!  Author  : Song Yezhi      
	!  verison : 2020-5-24 16:54
	!
	------------------------------------------------------*/
	M, N := MatDim(A)
	for i := 0; i < M; i++ {
		for j := 0; j < N; j++ {
			fmt.Printf("%18.6f   ", A[i][j])
		}
		fmt.Printf("\n")
	}
}

func LDL(A [][]float64) (L [][]float64, D []float64) {
	/*------------------------------------------------------
	  !  Author  : Song Yezhi      
	  !  verison : 2020-5-24 17:24
	  !
	  !    A=LDL'
	  !  -----------------------------------------------------
	  !  Input  Parameters :
	  !    A
	  !  Output Parameters :
	  !    L
	  !    D
	  ------------------------------------------------------*/
	M, N := MatDim(A)
	if M != N {
		fmt.Printf("the rows and cols is not equal")
	}

	G := Mat(N, N)

	L = Mat(N, N)
	D = Vec(N)

	var i, j, k int

	D[0] = A[0][0]

	var tmp1 float64

	for i = 1; i < N; i++ {
		for j = 0; j <= i-1; j++ {
			tmp1 = 0.0
			for k = 0; k <= j-1; k++ {
				tmp1 = tmp1 + G[i][k]*L[j][k]
			}
			G[i][j] = A[i][j] - tmp1
		}

		for j = 0; j <= i-1; j++ {
			L[i][j] = G[i][j] / D[j]
		}

		//comput D[i]
		tmp1 = 0.0

		for k = 0; k <= i-1; k++ {
			tmp1 = tmp1 + G[i][k]*L[i][k]
		}

		D[i] = A[i][i] - tmp1

	}

	//matOutput(L)

	//set the diag elements
	for i = 0; i < N; i++ {
		L[i][i] = 1.0
	}
	return
}

测试代码如下

// mattest.go
package main

import (
	"bufio"
	"fmt"
	"io"
	"os"
	"strconv"
	"strings"
)

func main() {
	//test_vDot_matVec()
	//
	//testmatmul()

	testLDL()

}

func test_vDot_matVec() {
	// vdot
	// matvec
	x := []float64{1, 2, 3}
	y := []float64{1, 2, 3}

	z := VecDot(x, y)

	KK := [][]int{{1, 2, 3}, {4, 5, 6}}

	fmt.Printf("%d\n", len(KK))
	fmt.Printf("%d", len(KK[0]))

	fmt.Printf("%12.5f \n", z)

	//----------------------------

	A := [][]float64{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}, {10, 11, 12}}
	b := MatVec(A, x)
	K := len(b)
	for i := 0; i < K; i++ {
		fmt.Printf("%12.5f \n", b[i])
	}
	//-------------------------------

}

func testmatmul() {
	//test the matrix multiply

	fi, _ := os.Open("matmul.txt")
	defer fi.Close()

	var M, N, P int
	var A [][]float64
	var B [][]float64
	var C [][]float64
	var word []string
	var i, j int

	br := bufio.NewReader(fi)
	for {
		var line string
		line, ioerr := br.ReadString('\n')
		if ioerr == io.EOF {
			break
		}

		if strings.TrimSpace(line) == "" {
			continue
		}
		// 这很重要,否则空格行分割要报错

		word = strings.Fields(line)

		if word[0] == "DIM" {
			line, _ = br.ReadString('\n')
			word = strings.Fields(line)
			M, _ = strconv.Atoi(word[0])
			N, _ = strconv.Atoi(word[1])
			P, _ = strconv.Atoi(word[2])
			fmt.Printf("M=%3d N=%3d P=%3d \n", M, N, P)
		}

		break

	}

	A = Mat(M, N)
	B = Mat(N, P)

	for {

		//var line string
		line, ioerr := br.ReadString('\n')
		if ioerr == io.EOF {
			break
		}

		if strings.TrimSpace(line) == "" {
			continue
		}
		// 这很重要,否则空格行分割要报错

		word = strings.Fields(line)

		if word[0] == "A" {

			//fmt.Printf("A = \n")

			for i = 0; i < M; i++ {

				line, ioerr = br.ReadString('\n')
				if ioerr == io.EOF {
					break
				}

				//fmt.Printf("i=%d  %s", i, line)

				word = strings.Fields(line)

				//fmt.Printf("%s", word)

				for j = 0; j < N; j++ {
					A[i][j], _ = strconv.ParseFloat(word[j], 64)

					//fmt.Printf("%18.6f ", A[i][j])

				}
				//fmt.Printf("\n")
			}

		}

		if word[0] == "B" {
			//fmt.Printf("\nB= \n")
			for i = 0; i < N; i++ {
				line, ioerr = br.ReadString('\n')
				if ioerr == io.EOF {
					break
				}
				word = strings.Fields(line)
				for j = 0; j < P; j++ {
					B[i][j], _ = strconv.ParseFloat(word[j], 64)
					//fmt.Printf(" %18.6f ", B[i][j])
				}
				//fmt.Printf("\n")
			}

		}

	}

	fmt.Printf("A= \n")
	MatOutput(A)
	fmt.Printf("B= \n")
	MatOutput(B)
	fmt.Printf("C= \n")
	C = MatMul(A, B)
	for i = 0; i < M; i++ {
		for j = 0; j < P; j++ {
			fmt.Printf("%18.6f", C[i][j])
		}
		fmt.Printf("\n")
	}

}

func testLDL() {

	fi, _ := os.Open("ldl.txt")
	defer fi.Close()

	//fout, _ := os.Open("ldl.out")

	var N int
	//the dimention of matrix A
	var A [][]float64

	var word []string

	br := bufio.NewReader(fi)

	for {

		//var line string
		line, ioerr := br.ReadString('\n')
		if ioerr == io.EOF {
			break
		}

		if strings.TrimSpace(line) == "" {
			continue
		}
		// 这很重要,否则空格行分割要报错

		word = strings.Fields(line)

		if word[0] == "dim" {
			N, _ = strconv.Atoi(word[1])
			fmt.Printf("N=%3d  \n", N)
		}
		break
		//it is important
	}

	A = Mat(N, N)

	for {
		var line string
		line, ioerr := br.ReadString('\n')
		if ioerr == io.EOF {
			break
		}

		if strings.TrimSpace(line) == "" {
			continue
		}

		word = strings.Fields(line)

		if word[0] == "A" {

			for i := 0; i < N; i++ {
				line, ioerr = br.ReadString('\n')
				if ioerr == io.EOF {
					break
				}
				word = strings.Fields(line)
				for j := 0; j < N; j++ {
					A[i][j], _ = strconv.ParseFloat(word[j], 64)
				}

			}
		}

	}

	fmt.Printf("A= \n")
	MatOutput(A)

	var L [][]float64
	var D []float64

	L, D = LDL(A)

	fmt.Printf("L=\n")
	MatOutput(L)

	fmt.Printf("D=\n")
	VecOutput(D)

	//write to file
	fout, _ := os.OpenFile("ldl_out.txt", os.O_RDWR|os.O_CREATE, os.ModePerm)
	defer fout.Close()

	var Nbyte int
	Nbyte, _ = fout.WriteString("L=\n")
	//fmt.Printf("%d", Nb)

	for i := 0; i < N; i++ {
		for j := 0; j < N; j++ {
			Nbyte, _ = fout.WriteString(fmt.Sprintf("%12.6f", L[i][j]))
		}
		Nbyte, _ = fout.WriteString("\n")
	}

	Nbyte, _ = fout.WriteString("D=\n")

	for i := 0; i < N; i++ {
		Nbyte, _ = fout.WriteString(fmt.Sprintf("%12.6f\n", D[i]))
	}

	fmt.Printf("%d", Nbyte)

}

测试数据如下,其文件名为"ldl.txt"

dim   5
    
A
   2.009812444224590   2.060104739664038   2.332961237400925   2.026604634542785   1.534406401310821
   2.060104739664038   2.801229204101148   2.819142292276107   2.334663017207164   1.741469892899938
   2.332961237400925   2.819142292276107   3.320797559925783   2.763883678670930   1.829446907264691
   2.026604634542785   2.334663017207164   2.763883678670930   2.940684174766706   1.647703419445861
   1.534406401310821   1.741469892899938   1.829446907264691   1.647703419445861   1.347530322398982 
   
ans   
L 
   1.000000000000000                   0                   0                   0                   0
   1.025023377471848   1.000000000000000                   0                   0                   0
   1.160785546982226   0.620386905632990   1.000000000000000                   0                   0
   1.008355103167188   0.373195635460327   0.724897356264742   1.000000000000000                   0
   0.763457508545188   0.244596718902562  -0.162120943577746   0.126656663999231   1.000000000000000
D
   2.009812444224590   0.689573685904954   0.347326813862890   0.618594720146081   0.115768425216063