Cholesky矩阵分解

TOC

• 介绍(introductions)
• 算法(algorithm)
• 代码(code)

1. 介绍

Cholesky 分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解。它要求矩阵的所有特征值必须大于零，故分解的下三角的对角元也是大于零的。Cholesky分解法又称平方根法，是当A为实对称正定矩阵时，LU三角分解法的变形。其中L称之为矩阵A的Cholesky 因子。

2. 算法

以3x3的矩阵为例，算法如下

3. 代码

fortran 代码如下：

!  Cholesky_decomposition.f90
!
!****************************************************************************
!
!  PROGRAM: Cholesky_decomposition
!
!  PURPOSE:  Use the Cholesky method to decopose a positive definite matrix ,S=L*L'
!  here L' is transpose matrix of L. We then obtain the inverse matrix of S and L
!  matrix as an output.
!
!****************************************************************************

    program Cholesky_decomposition

    implicit none

    ! Variables
    integer::i,j,k
    real(kind=8)::summ
    real(kind=8)::S(3,3)=(/ 1.0,1.0,1.0,1.0,5.0,3.0,1.0,3.0,11.0 /)
    real(kind=8)::L(3,3),L_inv(3,3),L_inv_T(3,3),S_inv(3,3)
    real(kind=8),intrinsic::sqrt,transpose,matmul
    ! Body of Cholesky_decomposition for the positive definite matrix S
    write(*,*),"S matrix"
    write(*,"(3(f12.6,2x))"),S
 !   L(1,1)=sqrt(S(1,1))
    do i=1,3
         do j=1,3
                     
                 if(i>j)then
                     summ=0.0
                     do k=1,j-1
                         summ=summ+L(i,k)*L(j,k)
                     end do
                     L(i,j)=(S(i,j)-summ)/L(j,j)
                 elseif(i==j)then
                     summ=0.0
                     do k=1,i-1
                         summ=summ+L(i,k)**2
                     end do
                     L(i,j)=sqrt(S(i,j)-summ)
                 else
                     L(i,j)=0.0
                 end if
                 
           
         end do
    end do
    write(*,*)"L matrix"
    write(*,"(3(f12.6,2x))"),((L(j,i),i=1,3),j=1,3)
    !!calculate the inverse matrix of L matrix ,at the same time
    !!  S matrix by be sovled with sqare root method
    do i=1,3
         do j=1,3
            if(i==j)then
                L_inv(i,j)=1/L(i,j)
            elseif(i>j)then
                summ=0.0
                do k=1,i-1
                    summ=summ+L(i,k)*L_inv(k,j)
                end do
                L_inv(i,j)=-1*summ/L(i,i)
            else
                L_inv(i,j)=0.0
            end if
         end do
    end do
    write(*,*)"L_inv matrix"
    write(*,"(3(f12.6,2x))"),((L_inv(j,i),i=1,3),j=1,3)
    L_inv_T=transpose(L_inv)
    S_inv=matmul(L_inv_T,L_inv)
    write(*,*)"S_inv matrix"
    write(*,"(3(f12.6,2x))"),((S_inv(j,i),i=1,3),j=1,3)
    end program Cholesky_decomposition

通过ifort进行编译,执行

ifort Cholesky_decomposition.f90 -o cd.x
./cd.x

可以得到如下结果：

S matrix
   1.000000      1.000000      1.000000
   1.000000      5.000000      3.000000
   1.000000      3.000000     11.000000
L matrix
   1.000000      0.000000      0.000000
   1.000000      2.000000      0.000000
   1.000000      1.000000      3.000000
L_inv matrix
   1.000000      0.000000      0.000000
  -0.500000      0.500000      0.000000
  -0.166667     -0.166667      0.333333
S_inv matrix
   1.277778     -0.222222     -0.055556
  -0.222222      0.277778     -0.055556
  -0.055556     -0.055556      0.111111