Python实现修正cholesky分解

Cholesky分解针对的是对称非正定矩阵

主要目的包括将原矩阵强制转换为对称正定矩阵以及保证分解过程的稳定性
修正cholesky分解的做法
给定对称矩阵A∈Rn×nA \in R^{n\times n}A∈Rn×n：计算出
δ=maxi=1,…,n∣Aii∣\delta=max_{i=1,\ldots,n}|A_{ii}|δ=maxi=1,…,n∣Aii∣
ν=maxi≠j∣Ai,j∣\nu=max_{i\neq j}|A_{i,j}|ν=maxi=j∣Ai,j∣
β2=max(δ,ν(n2−1),μ)\beta^2=max(\delta,\frac{\nu}{\sqrt(n^2-1)},\mu)β2=max(δ,(n2−1)ν,μ)
1:dj′=max{δ,∣Aj,j−∑r=1j−1Cj,r2dr−1∣}d_{j}^{'}=max\{\delta,|A_{j,j}-\sum_{r=1}^{j-1}C_{j,r}^{2}d_{r}^{-1}|\}dj′=max{δ,∣Aj,j−∑r=1j−1Cj,r2dr−1∣}
2:Ci,j=Ai,j−∑r=1j−1Lj,rCi,r,i≥j+1,…,nC_{i,j}=A_{i,j}-\sum_{r=1}^{j-1}L_{j,r}C_{i,r},i\geq j+1,\ldots,nCi,j=Ai,j−∑r=1j−1Lj,rCi,r,i≥j+1,…,n
3:θj=max{∣Ci,j,i>j},(j=n,θj=0)\theta_{j}=max\{|C_{i,j},i>j\},(j=n,\theta_j=0)θj=max{∣Ci,j,i>j},(j=n,θj=0)
4:dj=max{dj′,θj2β2}d_j=max\{d_{j}^{'},\frac{\theta_{j}^{2}}{\beta^{2}}\}dj=max{dj′,β2θj2}
5:Li,j=Ci,jdj,i≥j+1,…,nL_{i,j}=\frac{C_{i,j}}{d_{j}},i\geq j+1,\ldots,nLi,j=djCi,j,i≥j+1,…,n
上面步骤中，目的是为了输出矩阵L和矩阵D，其中L是单位下三角矩阵，D是对角矩阵，D=diag{d1,d2,…,dn}D=diag\{d_1,d_2,\ldots,d_n\}D=diag{d1,d2,…,dn}
修正以后的矩阵Axiu=LDLTA_{xiu}=LDL^{T}Axiu=LDLT，可以证明Axiu−AA_{xiu}-AAxiu−A是一个对角矩阵。下面我们放代码：

#修正cholesky分解
def XC(A):n = A.shape[0]dia = np.diagonal(A.numpy())#以列表存储矩阵对角元素delta = torch.tensor(max(abs(dia)))#delta = (max(abs(dia))).clone().detach()U = np.triu(A.numpy()) - np.diag(dia)#取出非对角上三角函数#print(abs(U).max())be = max(delta,abs(U).max()/(n**2 - 1))beta = (1/be).clone().detach()#根据A给出beta^2分之一的值C = torch.zeros_like(A);L = torch.eye(n,n);D = torch.zeros_like(A)d = torch.zeros(n);theta = torch.zeros(n)#----------------------d[0] = max(delta,abs(A[0,0]))for i in range(1,n):C[i,0] = A[i,0]   theta[0] = max(abs(C[:,0]))#theta_0D[0,0] = max(beta*theta[0]**2,d[0])for i in range(1,n):L[i,0] = C[i,0]/D[0,0]#------------------------------for j in range(1,n):for i in range(j+1,n):for r in range(j):d[j] = A[j,j] - C[j,r]**2/d[r]C[i,j] = A[i,j] - L[j,r]*C[i,r]d[j] = max(delta,abs(d[j]))theta[j] = max(abs(C[:,j]))D[j,j] = max(d[j],beta*theta[j]**2)for m in range(j + 1,n):L[m,j] = C[m,j]/D[j,j]#print(beta)return L,D
A = torch.eye(3,3)
A[0,1] = 1.0;A[0,2] = 2.0
A[1,0] = 1;A[1,1] += 1e-7;A[1,2] = 3
A[2,0] = 2;A[2,1] = 3
print(A)
L,D = XC(A)
print(L@D@L.t() - A)

Python实现修正cholesky分解相关推荐

【数值分析】Doolittle分解和Cholesky分解的Python实现
Doolittle 分解 import numpy as np # A=[[1.0,2.0,-3.0], # [2.0,-1.0,3.0], # [3.0,-2.0,2.0]]# A=[ [2, 4, ...
Matlab与线性代数--矩阵的Cholesky分解
本图文介绍了Matlab对正交矩阵的Cholesky分解操作.
线性代数笔记： Cholesky分解
1 介绍当一个实矩阵A是对称正定矩阵的时候,它可以分解成一个下三角矩阵L以及它的转置的乘积,即: 1.1 矩阵半正定的情况如果矩阵是正定的话,那么L唯一确定:如果矩阵是半正定的话,那么也可以分解, ...
Python中矩阵SVD分解及还原
python中SVD分解及还原: import numpy as np from numpy import linalg as la S = np.zeros([5,5]) A=np.random.r ...
matlab和python中的svd分解的区别
matlab中的svd分解中得到的第三项是V, 代码如下: A=[1 2 3 4 5 6 7 8 9; 5 6 7 8 9 0 8 6 7; 9 0 8 7 1 4 3 2 1; 6 4 2 1 3 ...
cholesky分解java代码_Cholesky 分解（转）
Cholesky 分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解. 它要求矩阵的所有特征值必须大于零,故分解的下三角的对角元也是大于零的. Cholesky分解法又称平方根法,是 ...
矩阵的Cholesky分解
原文链接:矩阵的Cholesky分解首先来复习线性代数中几个重要的概念. 1)如果一个复矩阵A = A*(共轭转置),则A称为Hermitian矩阵.(注意,矩阵A转置后仍为其本身,显然A一定是方阵 ...
使用 uBLAS 进行实对称正定矩阵的 Cholesky 分解
Cholesky 分解理论矩阵分解--三角分解(Cholesky 分解) 矩阵分解--三角分解(二) 注:只有实对称矩阵才有 Cholesky 分解理论. 已知实对称正定矩阵 AA,其 Choles ...
矩阵分解——三角分解（Cholesky 分解）
(1)一个对角元素都是1的下三角矩阵,称为单位下三角矩阵. (2)上(下)三角矩阵的乘积仍是上(下)三角矩阵: (3)一般来说,矩阵的三角分解不唯一. (4)实对称正定矩阵 AA,Δk>0\De ...
cholesky分解_Time Series Analysis-1.2 LDL分解
最近考完两个小quiz,停了一段时间,今晚抽空继续来分享这门课的笔记. 1.前言上一期分享了Cholesky分解的基本步骤和伪代码,本期介绍另外一种矩阵分解的方法--LDL分解. 首先补充一下,近几 ...

Python实现修正cholesky分解

Cholesky分解针对的是对称非正定矩阵

Python实现修正cholesky分解相关推荐

最新文章

热门文章