一、简介
Cholesky 分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解。它要求矩阵的所有特征值必须大于零,故分解的下三角的对角元也是大于零的。Cholesky分解法又称平方根法,是当A为实对称正定矩阵时,LU三角分解法的变形。
如果矩阵A为n阶对称正定矩阵,则存在一个对角元素为正数的下三角实矩阵L,当限定L的对角元素为正时,这种分解是唯一的,称为Cholesky分解。
详解
二、实现
# -*- coding: utf-8 -*-
"""
Created on Wed Dec 14 19:43:58 2016
Cholesky分解
@author: Administrator
"""
from numpy import *
def cholesky_decomp(A):
m,n = A.shape
if m != n:
return
temp_A = A.copy()
L = zeros((n,n),dtype=float)
for i in range(0,n):
for j in range(0,n):
if i >= j + 1:
if j > 0:
L[i,j] = (temp_A[i,j] - sum(L[i,0:j] * L[j,0:j])) / L[j,j]
else:
L[i,j] = temp_A[i,j] / L[j,j]
else:
if j > 0:
L[j,j] = sqrt(temp_A[j,j] - sum(L[j,0:j] * L[j,0:j]))
else:
L[j,j] = sqrt(A[j,j])
return L
def cal_cholesky_decomp(A,b):
m,n = A.shape
L = cholesky_decomp(A)
y = array([0 for i in range(n)],dtype=float)
x = y.copy()
for i in range(0,n):
if i == 0:
y[i] = b[i] / L[i,i]
else:
y[i] = (b[i] - sum(L[i,0:i] * y[0:i].transpose())) / L[i,i]
for i in range(n-1,-1,-1):
if i == n-1:
x[i] = y[i] / L[i,i]
else:
x[i] = (y[i] - sum(L[i+1:n,i] * x[i+1:n].transpose())) / L[i,i]
return x
if __name__ == '__main__':
n = 10
A = ones((n,n),dtype=float)
b = arange(0,n)
for i in range(n):
for j in range(n):
A[i,j] = 1. / (i + j + 1) # i +j - 1 + 1 + 1
print 'n=',n,':',cal_cholesky_decomp(A,b)