计算数学
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