本站原创文章,转载请说明来自《老饼讲解-机器学习》www.bbbdata.com
本文先介绍什么是Cholesky矩阵分解及分解的方法
然后讲解相关公式的推导及及代码实现
本节讲解什么是Cholesky矩阵分解和如何对矩阵进行Cholesky分解
Cholesky分解的定义
Cholesky分解就是将正定对称矩阵分解成两个三角
A是正定对称矩阵,
求一个下三角矩阵 L,使得
即将A分解成
,
Cholesky矩阵分解流程
L的求解流程,是逐列求解,
即先求第1列,再求第2列,再求第3列....
每列的求解分为两步:
👉 先求对角元素
👉 再求非对角元素
假设前i-1列已求出,求第i列
求解 L第i列对角元素
的求解公式如下
即第i个对角元素为:根号(Aii-L第i行前(i-1)个元素的平方和)
举例:
求解 L第i列所有元素
在求得后,第 i 列的求解公式如下
即L第i列为:(A的第i列- L前(i-1)列*各自第i行的元素)/
矩阵形式为:
✍️关于首列的计算公式
特别地,第1列的计算公式为
本节讲解上节Cholesky分解流程中的公式是如何推导出来的
Cholesky对角元素的公式推导
对于L的对角元素,
公式的推导流程如下
Aii等于L的 i 行乘的 i 列
的 i 列就是L的第i行的转置
i之后的元素是0,所以只要累加到i即可
根据上式,即可得到
Cholesky非角元素的公式推导
假设前(i-1)列已求出,且也已求得,
现求L第i列的所有元素
A的i列为 L乘以的 i 列
L与互为转置
A的i列为L每列乘以各自第i行元素后再相加
L的i行i列之后的元素为0,故不需i之后的列
独立抽出第i列
根据上式,即可得到
将上式写成矩阵形式,则为
本节展示如何用用python代码实现Cholesky分解
python代码实现Cholesky分解
# -*- coding: utf-8 -*-
"""
Cholesky分解的代码实现
"""
import numpy as np
def CholeskyDecompose(A):
r,c=A.shape
L=np.zeros((r,r))
L[0,0] = np.sqrt(A[0,0])
L[:,0] = A[:,0]/L[0,0]
i=4
for i in range(1,r):
L[i,i] = np.sqrt(A[i,i]-sum(L[i,:i+1]*L[i,:i+1]))
L[:,i] = (A[:,i]-L[:,:i]@L[i,:i].T)/L[i,i]
return L
# 测试样例
if __name__ == "__main__":
# 生成一个下三角矩阵L
np.random.seed(0)
rnd_mat = np.random.randint(0,10,(5,5))
L_true = np.tril(rnd_mat,0)
#将L*L^T组合成A
A = L_true@L_true.T
# 调用CholeskyDecompose对A进行Cholesky分解
L = CholeskyDecompose(A)
#打印结果
print('----生成的L-----')
print(L)
print('\n----LL^T组合成的A-----')
print(A)
print('\n----从A分解得到的L-------')
print(L)
代码运行结果如下
----生成的L-----
[[5. 0. 0. 0. 0.]
[9. 3. 0. 0. 0.]
[7. 6. 8. 0. 0.]
[6. 7. 7. 8. 0.]
[5. 9. 8. 9. 4.]]
----LL^T组合成的A-----
[[ 25 45 35 30 25]
[ 45 90 81 75 72]
[ 35 81 149 140 153]
[ 30 75 140 198 221]
[ 25 72 153 221 267]]
----从A分解得到的L-------
[[5. 0. 0. 0. 0.]
[9. 3. 0. 0. 0.]
[7. 6. 8. 0. 0.]
[6. 7. 7. 8. 0.]
[5. 9. 8. 9. 4.]]
参考
【1】Cholesky分解法 : https://blog.csdn.net/acdreamers/article/details/44656847
End