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