Cholesky矩阵分解¶

约 389 个字 82 行代码预计阅读时间 2 分钟

计算原理¶

对于对称正定矩阵 \(A\)，存在于一个实的下三角矩阵\(L\)，使得\(A=LL^T\)，若限定\(L\)的对角元素为正时，这种分解是唯一的。

\[ { \left( \begin{array} { l l l l } { a _ { 1 1 } } & { a _ { 1 2 } } & { \cdots } & { a _ { 1 n } } \\ { a _ { 2 1 } } & { a _ { 2 2 } } & { \cdots } & { a _ { 2 n } } \\ { \vdots } & & { \ddots } & \\ { a _ { n 1 } } & { a _ { n 2 } } & { \cdots } & { a _ { n n } } \end{array} \right) } = { \left( \begin{array} { l l l l } { l _ { 1 1 } } & & & \\ { l _ { 2 1 } } & { l _ { 2 2 } } & & \\ { \cdots } & & { \ddots } & \\ { l _ { n 1 } } & { l _ { n 2 } } & { \cdots } & { l _ { n n } } \end{array} \right) } { \left( \begin{array} { l l l l } { l _ { 1 1 } } & { l _ { 2 1 } } & { \cdots } & { l _ { n 1 } } \\ & { l _ { 2 2 } } & { \cdots } & { l _ { n 2 } } \\ & & { \ddots } & \\ & & & { l _ { n n } } \end{array} \right) } \]

\(l_{11} = \sqrt{a_{11}}\)
\(l_{i1} = \frac{a_{i1}}{l_{11}}, \quad i=2,N\)
对 \(j=2,N\)，做A~B处理：

A:

\[ l_{jj} = \sqrt{a_{jj} - \sum_{k=1}^{j-1} l_{jk}^2 } \]

B:

\[ l_{ij} = \frac{ \left( a_{ij} - \sum_{k=1}^{j-1} l_{ik} l_{jk} \right) }{ l_{jj} }, \quad i=j+1,N \]

Fortran数值实现¶

subroutine cholesky(a,l,n)
! cholesky矩阵分解(只适用于对称正定矩阵)
    ! A = L*L^T
    integer, intent(in) :: n
    real*8, intent(in) :: a(n,n)
    real*8, intent(out) :: l(n,n)

    integer :: i,j,k

    l = 0
    l(1,1) = sqrt(a(1,1))
    l(2:,1) = a(2:,1)/l(1,1)

    do j = 2,n
        s = 0
        do k = 1,j-1
            s = s + l(j,k)**2
        end do
        l(j,j) = sqrt(a(j,j)-s)
        do i = j+1,n
            s = 0
            do k = 1,j-1
                s = s + l(i,k) * l(j,k)
            end do
            l(i,j) = (a(i,j)-s)/l(j,j)
        end do
    end do
end subroutine cholesky

数值案例¶

\[ \mathbf{A} = \left( \begin{array}{cccc} 233.4615 & 113.8423 & 256.0623 & 145.0697 \\ 113.8423 & 78.6033 & 127.4298 & 95.3089 \\ 256.0623 & 127.4298 & 281.4721 & 164.8676 \\ 145.0697 & 95.3089 & 164.8676 & 181.2339 \end{array} \right) \]

主程序：

program main
    implicit none
    integer,  parameter :: n = 4
    real*8 :: a(n,n),l(n,n)
    integer :: i, j

    a = reshape([ &
        233.4615d0, 113.8423d0, 256.0623d0, 145.0697d0, &
        113.8423d0,  78.6033d0, 127.4298d0,  95.3089d0, &
        256.0623d0, 127.4298d0, 281.4721d0, 164.8676d0, &
        145.0697d0,  95.3089d0, 164.8676d0, 181.2339d0  &
    ], shape(a), order=[2,1])

    call cholesky(a, l, n)

    print *, "Lower Triangular Matrix L from Cholesky decomposition:"
    do i = 1, n
        write(*, 100) (l(i,j), j = 1, n)
    end do

100 format(4F12.6)
end program main

输出：

 Lower Triangular Matrix L from Cholesky decomposition:
   15.279447    0.000000    0.000000    0.000000
    7.450682    4.805272    0.000000    0.000000
   16.758610    0.534147    0.579450    0.000000
    9.494434    5.112904    5.217081    6.142468

求解矩阵方程¶

对于\(\mathbf { A x } = \mathbf { b }\)，进行cholesky分解：

\[ \mathbf { L } \mathbf { L } ^ { T } \mathbf { x } = \mathbf { b }，等价于\left\{ \begin{array} { l l } { \mathbf { L } \mathbf { y } = \mathbf { b } } \\ { \mathbf { L } ^ { \mathrm { T } } \mathbf { x } = \mathbf { y } } \end{array} \right. \]

对于：

\[ \mathbf { A } = { \left( \begin{array} { l l l } { 4 } & { - 1 } & { 1 } \\ { - 1 } & { 4 . 2 5 } & { 2 . 7 5 } \\ { 1 } & { 2 . 7 5 } & { 3 . 5 } \end{array} \right) } , \mathbf { b } = { \left( \begin{array} { l } { 4 } \\ { 6 } \\ { 7 . 2 5 } \end{array} \right) } \]

program main
    implicit none
    integer,  parameter :: n = 3
    real*8 :: a(n,n), l(n,n), b(n), x(n), y(n)
    integer :: i

    a = reshape([ &
        4.0d0,   -1.0d0,   1.0d0,  &
        -1.0d0,  4.25d0,  2.75d0, &
        1.0d0,   2.75d0,  3.5d0   &
    ], shape(a), order=[2,1])

    b = [4.0d0, 6.0d0, 7.25d0]

    call cholesky(a, l, n)
    ! 解 Ly = b
    call lowtri(l, b, y, n)
    ! 解 L^T x = y
    call uptri(transpose(l), y, x, n)

    print *, "Solution vector x:"
    write(*, 100) x

100 format(3F12.6)
end program main

输出：

 Solution vector x:
    1.000000    1.000000    1.000000

注意：子程序中默认：implicit real*8(a-z)

参考文献¶

宋叶志,茅永兴,赵秀杰.Fortran 95/2003科学计算与工程[M].清华大学出版社,2011.