线性代数 中,科列斯基分解 (英语:Cholesky decomposition 或 Cholesky factorization )是指将一个正定 的埃尔米特矩阵 分解成一个下三角矩阵 与其共轭转置 之乘积 。这种分解方式在提高代数运算效率、蒙特卡罗方法 等场合中十分有用。实数 矩阵 的科列斯基分解由安德烈-路易·科列斯基 最先发明。实际应用中,科列斯基分解在求解线性方程组 中的效率约两倍于LU分解 。[ 1]
描述
对正定埃尔米特矩阵
A
{\displaystyle \mathbf {A} }
进行科列斯基分解,即求矩阵
L
{\displaystyle \mathbf {L} }
使下式成立
A
=
L
L
∗
{\displaystyle \mathbf {A} =\mathbf {LL} ^{*}}
其中,
L
{\displaystyle \mathbf {L} }
是一个下三角矩阵且所有对角元素均为正实数 ,
L
∗
{\displaystyle \mathbf {L} ^{*}}
表示
L
{\displaystyle \mathbf {L} }
的共轭转置。每一个正定埃尔米特矩阵都有一个唯一的科列斯基分解。[ 2]
当矩阵
A
{\displaystyle \mathbf {A} }
是一个半正定的埃尔米特矩阵,若允许
L
{\displaystyle \mathbf {L} }
的对角线元素为零,则
A
{\displaystyle \mathbf {A} }
也存在上述形式的分解。[ 3]
当
A
{\displaystyle \mathbf {A} }
为实数矩阵,则
L
{\displaystyle \mathbf {L} }
也为实数矩阵且科列斯基分解可改写成
A
=
L
L
T
{\displaystyle \mathbf {A} =\mathbf {LL} ^{\mathbf {T} }}
。[ 4]
当
A
{\displaystyle \mathbf {A} }
是正定矩阵 时,科列斯基分解是唯一的,即只存在一个对角元素均严格大于零的下三角矩阵,使
A
=
L
L
∗
{\displaystyle \mathbf {A} =\mathbf {LL} ^{*}}
成立。然而,当
A
{\displaystyle \mathbf {A} }
是半正定时,分解则不一定是唯一的。
定理的逆命题 自然成立:对于某些可逆矩阵
L
{\displaystyle \mathbf {L} }
(下三角矩阵或其他矩阵),如果
A
{\displaystyle \mathbf {A} }
可被写成
L
L
∗
{\displaystyle \mathbf {LL} ^{*}}
,则
A
{\displaystyle \mathbf {A} }
是一个正定的埃尔米特矩阵。
LDL分解
经典科列斯基分解的一个变形是LDL分解,即
A
=
L
D
L
∗
{\displaystyle \mathbf {A} =\mathbf {LDL} ^{*}}
其中,
L
{\displaystyle \mathbf {L} }
是一个单位下三角矩阵 ,
D
{\displaystyle \mathbf {D} }
是一个对角矩阵 。
该分解与经典科列斯基分解犹有关系,如下:
A
=
L
D
L
∗
=
L
D
1
2
(
D
1
2
)
∗
L
∗
=
L
D
1
2
(
L
D
1
2
)
∗
{\displaystyle \mathbf {A} =\mathbf {LDL} ^{*}=\mathbf {LD} ^{\frac {1}{2}}(\mathbf {D} ^{\frac {1}{2}})^{*}\mathbf {L} ^{*}=\mathbf {LD} ^{\frac {1}{2}}(\mathbf {LD} ^{\frac {1}{2}})^{*}}
或者,由于
L
{\displaystyle \mathbf {L} }
的对角元素必须为
1
{\displaystyle 1}
,且
L
C
h
o
l
e
s
k
y
{\displaystyle \mathbf {L} ^{Cholesky}}
与
L
{\displaystyle \mathbf {L} }
都是下三角矩阵,故如果已知经典科列斯基分解
L
C
h
o
l
e
s
k
y
{\displaystyle \mathbf {L} ^{Cholesky}}
,则
L
D
L
T
{\displaystyle \mathbf {LDL} ^{T}}
形式可依下式求出,设S 是包含
L
C
h
o
l
e
s
k
y
{\displaystyle \mathbf {L} ^{Cholesky}}
的对角元素的对角矩阵,
D
=
S
2
{\displaystyle \mathbf {D} =\mathbf {S} ^{2}}
L
=
L
C
h
o
l
e
s
k
y
S
−
1
{\displaystyle \mathbf {L} =\mathbf {L} ^{Cholesky}\mathbf {S} ^{-1}}
LDL变形 如果得以有效执行,构造及使用时所需求的空间及计算的复杂性与经典科列斯基分解是相同的,但是可避免提取平方根 。[ 5] 某些不存在科列斯基分解的不定矩阵 ,也可以执行LDL分解,此时矩阵
D
{\displaystyle \mathbf {D} }
中会出现负数 元素。因此人们更倾向于使用LDL分解。对于实数矩阵,该种分解的形式可被改写成
A
=
L
D
L
T
{\displaystyle \mathbf {A} =\mathbf {LDL} ^{\mathbf {T} }}
此形式通常称为LDLT分解 (或LDLT 分解 )。它与实对称矩阵的特征分解 密切相关,因为对于实对称矩阵 ,存在特征分解
A
=
Q
Λ
Q
T
{\displaystyle \mathbf {A} =\mathbf {Q\Lambda Q} ^{\mathbf {T} }}
。
实例
以下乃一个实对称矩阵 的科列斯基分解︰
(
4
12
−
16
12
37
−
43
−
16
−
43
98
)
=
(
2
0
0
6
1
0
−
8
5
3
)
(
2
6
−
8
0
1
5
0
0
3
)
{\displaystyle {\begin{aligned}\left({\begin{array}{*{3}{r}}4&12&-16\\12&37&-43\\-16&-43&98\\\end{array}}\right)=\left({\begin{array}{*{3}{r}}2&0&0\\6&1&0\\-8&5&3\\\end{array}}\right)\left({\begin{array}{*{3}{r}}2&6&-8\\0&1&5\\0&0&3\\\end{array}}\right)\end{aligned}}}
以下乃其LDLT 分解︰
(
4
12
−
16
12
37
−
43
−
16
−
43
98
)
=
(
1
0
0
3
1
0
−
4
5
1
)
(
4
0
0
0
1
0
0
0
9
)
(
1
3
−
4
0
1
5
0
0
1
)
{\displaystyle {\begin{aligned}\left({\begin{array}{*{3}{r}}4&12&-16\\12&37&-43\\-16&-43&98\\\end{array}}\right)&=\left({\begin{array}{*{3}{r}}1&0&0\\3&1&0\\-4&5&1\\\end{array}}\right)\left({\begin{array}{*{3}{r}}4&0&0\\0&1&0\\0&0&9\\\end{array}}\right)\left({\begin{array}{*{3}{r}}1&3&-4\\0&1&5\\0&0&1\\\end{array}}\right)\end{aligned}}}
应用
科列斯基分解主要被用于线性方程组
A
x
=
b
{\displaystyle \mathbf {Ax} =\mathbf {b} }
的求解。如果A 是对称正定 的,我们可以先求出
A
=
L
L
T
{\displaystyle \mathbf {A} =\mathbf {LL} ^{\mathbf {T} }}
,随后借向后替换法 对y 求解
L
y
=
b
{\displaystyle \mathbf {Ly} =\mathbf {b} }
,再以向前替换法 对x 求解
L
T
x
=
y
{\displaystyle \mathbf {L} ^{\mathbf {T} }\mathbf {x} =\mathbf {y} }
即得最终解。
另一种可避免在计算
L
L
T
{\displaystyle \mathbf {LL} ^{\mathbf {T} }}
时需要解平方根的方法就是计算
A
=
L
D
L
T
{\displaystyle \mathbf {A} =\mathbf {LDL} ^{\mathrm {T} }}
,然后对y 求解
L
y
=
b
{\displaystyle \mathbf {Ly} =\mathbf {b} }
,最后求解
D
L
T
x
=
y
{\displaystyle \mathbf {DL} ^{\mathrm {T} }\mathbf {x} =\mathbf {y} }
对于可以被改写成对称矩阵的线性方程组,科列斯基分解及其LDL变形是一个较高效率及较高数值稳定性的求解方法。相比之下,其效率几近为LU分解 的两倍。[来源请求]
矩阵求逆
若欲对埃尔米特矩阵 直接求逆,可以通过科列斯基分解,类似求解线性方程组一般求出逆矩阵,这需要
n
3
{\displaystyle n^{3}}
次运算(
1
2
n
3
{\displaystyle {\frac {1}{2}}n^{3}}
次乘法运算)。 该方法即便要求逐步计算也非常有效率。
非埃尔米特矩阵
B
{\displaystyle \mathbf {B} }
也可以通过下例性质求逆,因为其中
B
B
∗
{\displaystyle \mathbf {BB} ^{*}}
总是埃尔米特矩阵︰
B
−
1
=
B
∗
(
B
B
∗
)
−
1
{\displaystyle \mathbf {B} ^{-1}=\mathbf {B} ^{*}(\mathbf {BB} ^{*})^{-1}}
计算
有各种各样的方法用于计算科列斯基分解。 常用的算法的计算复杂度在一般情况下为
O
(
n
3
)
{\displaystyle O(n^{3})}
。[来源请求]
下面的算法何者较缓存决于执行时的具体条件。总体而言,第一个算法会稍慢,因为其以一种不太寻常的方式读取数据。
科列斯基算法
用于计算矩阵
L
{\displaystyle \mathbf {L} }
的科列斯基算法 ,是以高斯消元法 为基础而调整得来的。
递归算法由
i
:=
1
{\displaystyle i:=1}
开始,令
A
(
1
)
:=
A
{\displaystyle \mathbf {A} ^{(1)}:=\mathbf {A} }
在步骤
i
{\displaystyle i}
,矩阵
A
(
i
)
{\displaystyle A^{(i)}}
的形式如下:
A
(
i
)
=
(
I
i
−
1
0
0
0
a
i
,
i
b
i
∗
0
b
i
B
(
i
)
)
,
{\displaystyle \mathbf {A} ^{(i)}={\begin{pmatrix}\mathbf {I} _{i-1}&0&0\\0&a_{i,i}&\mathbf {b} _{i}^{*}\\0&\mathbf {b} _{i}&\mathbf {B} ^{(i)}\end{pmatrix}},}
其中,
I
i
−
1
{\displaystyle \mathbf {I} _{i-1}}
代表
i
−
1
{\displaystyle i-1}
维的单位矩阵 。
此时定义矩阵
L
i
{\displaystyle \mathbf {L} _{i}}
为
L
i
:=
(
I
i
−
1
0
0
0
a
i
,
i
0
0
1
a
i
,
i
b
i
I
n
−
i
)
,
{\displaystyle \mathbf {L} _{i}:={\begin{pmatrix}\mathbf {I} _{i-1}&0&0\\0&{\sqrt {a_{i,i}}}&0\\0&{\frac {1}{\sqrt {a_{i,i}}}}\mathbf {b} _{i}&\mathbf {I} _{n-i}\end{pmatrix}},}
随即
A
(
i
)
{\displaystyle \mathbf {A} ^{(i)}}
可以被改写成
A
(
i
)
=
L
i
A
(
i
+
1
)
L
i
∗
{\displaystyle \mathbf {A} ^{(i)}=\mathbf {L} _{i}\mathbf {A} ^{(i+1)}\mathbf {L} _{i}^{*}}
其中
A
(
i
+
1
)
=
(
I
i
−
1
0
0
0
1
0
0
0
B
(
i
)
−
1
a
i
,
i
b
i
b
i
∗
)
.
{\displaystyle \mathbf {A} ^{(i+1)}={\begin{pmatrix}\mathbf {I} _{i-1}&0&0\\0&1&0\\0&0&\mathbf {B} ^{(i)}-{\frac {1}{a_{i,i}}}\mathbf {b} _{i}\mathbf {b} _{i}^{*}\end{pmatrix}}.}
注意︰ 在此
b
i
b
i
∗
{\displaystyle \mathbf {b} _{i}\mathbf {b} _{i}^{*}}
是一个外积 。
重复此步骤直到
i
{\displaystyle i}
从
1
{\displaystyle 1}
到
n
{\displaystyle n}
。
n
{\displaystyle n}
步之后,我们得到
A
(
n
+
1
)
=
I
{\displaystyle A^{(n+1)}=\mathbf {I} }
。因此,所求的下三角矩阵
L
{\displaystyle L}
为
L
:=
L
1
L
2
…
L
n
.
{\displaystyle \mathbf {L} :=\mathbf {L} _{1}\mathbf {L} _{2}\dots \mathbf {L} _{n}.}
科列斯基-巴那齐耶维茨及科列斯基-克劳特算法
科列斯基-巴那齐耶维茨算法以一个 5×5 矩阵执行的读取顺序(白色)及写入顺序(黄色)
写出等式
A
=
L
L
∗
{\displaystyle \mathbf {A} =\mathbf {LL} ^{*}}
,
A
=
L
L
T
=
(
L
11
0
0
L
21
L
22
0
L
31
L
32
L
33
)
(
L
11
L
21
L
31
0
L
22
L
32
0
0
L
33
)
=
(
L
11
2
(
symmetric
)
L
21
L
11
L
21
2
+
L
22
2
L
31
L
11
L
31
L
21
+
L
32
L
22
L
31
2
+
L
32
2
+
L
33
2
)
{\displaystyle {\begin{aligned}{\mathbf {A=LL^{T}} }&={\begin{pmatrix}\mathbf {L} _{11}&0&0\\\mathbf {L} _{21}&\mathbf {L} _{22}&0\\\mathbf {L} _{31}&\mathbf {L} _{32}&\mathbf {L} _{33}\\\end{pmatrix}}{\begin{pmatrix}\mathbf {L} _{11}&\mathbf {L} _{21}&\mathbf {L} _{31}\\0&\mathbf {L} _{22}&\mathbf {L} _{32}\\0&0&\mathbf {L} _{33}\end{pmatrix}}\\&={\begin{pmatrix}\mathbf {L} _{11}^{2}&&({\text{symmetric}})\\\mathbf {L} _{21}\mathbf {L} _{11}&\mathbf {L} _{21}^{2}+\mathbf {L} _{22}^{2}&\\\mathbf {L} _{31}\mathbf {L} _{11}&\mathbf {L} _{31}\mathbf {L} _{21}+\mathbf {L} _{32}\mathbf {L} _{22}&\mathbf {L} _{31}^{2}+\mathbf {L} _{32}^{2}+\mathbf {L} _{33}^{2}\end{pmatrix}}\end{aligned}}}
则得到矩阵
L
{\displaystyle \mathbf {L} }
的元素的计算公式如下︰
L
j
,
j
=
A
j
,
j
−
∑
k
=
1
j
−
1
L
j
,
k
2
{\displaystyle \mathbf {L} _{j,j}={\sqrt {A_{j,j}-\sum _{k=1}^{j-1}\mathbf {L} _{j,k}^{2}}}}
L
i
,
j
=
1
L
j
,
j
(
A
i
,
j
−
∑
k
=
1
j
−
1
L
i
,
k
L
j
,
k
)
,
for
i
>
j
{\displaystyle \mathbf {L} _{i,j}={\frac {1}{\mathbf {L} _{j,j}}}\left(A_{i,j}-\sum _{k=1}^{j-1}\mathbf {L} _{i,k}\mathbf {L} _{j,k}\right),\qquad {\text{for }}i>j}
只要
A
{\displaystyle \mathbf {A} }
是实数正定矩阵,则平方根 下的表达式恒为正。
对于复埃尔米特矩阵,则适用如下公式:
L
j
,
j
=
A
j
,
j
−
∑
k
=
1
j
−
1
L
j
,
k
L
j
,
k
∗
.
{\displaystyle \mathbf {L} _{j,j}={\sqrt {A_{j,j}-\sum _{k=1}^{j-1}\mathbf {L} _{j,k}\mathbf {L} _{j,k}^{*}}}.}
L
i
,
j
=
1
L
j
,
j
(
A
i
,
j
−
∑
k
=
1
j
−
1
L
i
,
k
L
j
,
k
∗
)
,
for
i
>
j
.
{\displaystyle \mathbf {L} _{i,j}={\frac {1}{\mathbf {L} _{j,j}}}\left(\mathbf {A} _{i,j}-\sum _{k=1}^{j-1}\mathbf {L} _{i,k}\mathbf {L} _{j,k}^{*}\right),\qquad {\text{for }}i>j.}
因此,要计算
L
i
,
j
(
i
≠
j
)
{\displaystyle \mathbf {L} _{i,j}(i\neq j)}
只需利用其的左、上方元素的值。计算通常是以以下其中一种顺序进行的。
科列斯基-巴那齐耶维茨算法 从矩阵L的左上角开始,依行进行计算。
科列斯基-克劳特算法 从矩阵L的左上角开始,依列进行计算。
若有需要,整个矩阵可以逐个元素计算得出 ,无论使用何种顺序读取。
计算的稳定性
假设我们要求解一个良置 的线性方程组。采用了LU分解 的算法,除非进行某种绕轴旋转,否则是不稳定的;如果算法进行了绕轴旋转,则其误差取决于矩阵的增长因子 ,这个数通常(但非总是)很小。
现在,假设算法适用科列斯基分解。如上所述,算法的效率将会是原来的两倍。此外,无必要进行绕轴旋转且误差总是很小。具体而言,若要求解
A
x
=
b
{\displaystyle \mathbf {Ax} =\mathbf {b} }
,
y
{\displaystyle \mathbf {y} }
表示已计算出的解,然后能解出干扰方程组
(
A
+
E
)
y
=
b
{\displaystyle (\mathbf {A} +\mathbf {E} )\mathbf {y} =\mathbf {b} }
,其中
‖
E
‖
2
≤
c
n
ε
‖
A
‖
2
.
{\displaystyle \|\mathbf {E} \|_{2}\leq c_{n}\varepsilon \|\mathbf {A} \|_{2}.}
在这里,
‖
‖
2
{\displaystyle \|\quad \|_{2}}
是指矩阵2-范数 ,而
c
n
{\displaystyle c_{n}}
是一个取决于
n
{\displaystyle n}
的小常量,
ε
{\displaystyle \varepsilon }
表示单位舍入 。
值得一提的是,科列斯基分解与平方根的使用有关。如果被分解的矩阵是正定的,那么只要运算精确 ,矩阵中带有平方根的元素的平方根下的数字永远是正数。不幸的是,由于存在舍入误差 ,这些数字可能为负数 ,并使算法搁浅。然而,此种情况仅当矩阵为病态 时才会出现。一种可解决此问题的方法,是增添一个对角修正矩阵至待分解矩阵,以增加矩阵的正定性。[ 6] 此法虽或将减少分解的准确性,但有在某些情况却颇有作为;譬如,执行应用于最优化的牛顿法 时,若初期值距最优值较远,则此时引入对角矩阵可以提高算法的稳定性。
LDL分解
科列斯基分解的另一种形式——LDL分解的计算方式如下所示,
A
=
L
D
L
T
=
(
1
0
0
L
21
1
0
L
31
L
32
1
)
(
D
1
0
0
0
D
2
0
0
0
D
3
)
(
1
L
21
L
31
0
1
L
32
0
0
1
)
=
(
D
1
(
s
y
m
m
e
t
r
i
c
)
L
21
D
1
L
21
2
D
1
+
D
2
L
31
D
1
L
31
L
21
D
1
+
L
32
D
2
L
31
2
D
1
+
L
32
2
D
2
+
D
3
.
)
{\displaystyle {\begin{aligned}{\mathbf {A=LDL} ^{\mathrm {T} }}&={\begin{pmatrix}1&0&0\\L_{21}&1&0\\L_{31}&L_{32}&1\\\end{pmatrix}}{\begin{pmatrix}D_{1}&0&0\\0&D_{2}&0\\0&0&D_{3}\\\end{pmatrix}}{\begin{pmatrix}1&L_{21}&L_{31}\\0&1&L_{32}\\0&0&1\\\end{pmatrix}}\\&={\begin{pmatrix}D_{1}&&(\mathrm {symmetric} )\\L_{21}D_{1}&L_{21}^{2}D_{1}+D_{2}&\\L_{31}D_{1}&L_{31}L_{21}D_{1}+L_{32}D_{2}&L_{31}^{2}D_{1}+L_{32}^{2}D_{2}+D_{3}.\end{pmatrix}}\end{aligned}}}
如果
A
{\displaystyle \mathbf {A} }
是实数矩阵,下述之递归计算式适用于矩阵
D
{\displaystyle \mathbf {D} }
及
L
{\displaystyle \mathbf {L} }
中的所有元素︰
D
j
=
A
j
j
−
∑
k
=
1
j
−
1
L
j
k
2
D
k
{\displaystyle D_{j}=A_{jj}-\sum _{k=1}^{j-1}L_{jk}^{2}D_{k}}
L
i
j
=
1
D
j
(
A
i
j
−
∑
k
=
1
j
−
1
L
i
k
L
j
k
D
k
)
,
for
i
>
j
.
{\displaystyle L_{ij}={\frac {1}{D_{j}}}\left(A_{ij}-\sum _{k=1}^{j-1}L_{ik}L_{jk}D_{k}\right),\qquad {\text{for }}i>j.}
如果
A
{\displaystyle \mathbf {A} }
是复埃尔米特矩阵 ,则矩阵
D
{\displaystyle \mathbf {D} }
及
L
{\displaystyle \mathbf {L} }
的计算适用以下公式:
D
j
=
A
j
j
−
∑
k
=
1
j
−
1
L
j
k
L
j
k
∗
D
k
{\displaystyle D_{j}=A_{jj}-\sum _{k=1}^{j-1}L_{jk}L_{jk}^{*}D_{k}}
L
i
j
=
1
D
j
(
A
i
j
−
∑
k
=
1
j
−
1
L
i
k
L
j
k
∗
D
k
)
,
for
i
>
j
.
{\displaystyle L_{ij}={\frac {1}{D_{j}}}\left(A_{ij}-\sum _{k=1}^{j-1}L_{ik}L_{jk}^{*}D_{k}\right),\qquad {\text{for }}i>j.}
同样地,若有需求,该读取顺序可以逐一计算矩阵中的每一个元素。
分块矩阵变形
当
A
{\displaystyle \mathbf {A} }
是一个不定矩阵时,LDL分解 在未经过谨慎的绕轴旋转的情况下是不稳定的;[ 7] 特别是,分解出的矩阵的元素可能是任意的。一个可取的改进手段是执行矩阵分块后再执行分解,通常将原矩阵分为包含
2
×
2
{\displaystyle 2\times 2}
的小矩阵的分块矩阵:[ 8]
A
=
L
D
L
T
=
(
I
0
0
L
21
I
0
L
31
L
32
I
)
(
D
1
0
0
0
D
2
0
0
0
D
3
)
(
I
L
21
T
L
31
T
0
I
L
32
T
0
0
I
)
=
(
D
1
(
s
y
m
m
e
t
r
i
c
)
L
21
D
1
L
21
D
1
L
21
T
+
D
2
L
31
D
1
L
31
D
1
L
21
T
+
L
32
D
2
L
31
D
1
L
31
T
+
L
32
D
2
L
32
T
+
D
3
)
{\displaystyle {\begin{aligned}{\mathbf {A=LDL} ^{\mathrm {T} }}&={\begin{pmatrix}\mathbf {I} &0&0\\\mathbf {L} _{21}&\mathbf {I} &0\\\mathbf {L} _{31}&\mathbf {L} _{32}&\mathbf {I} \\\end{pmatrix}}{\begin{pmatrix}\mathbf {D} _{1}&0&0\\0&\mathbf {D} _{2}&0\\0&0&\mathbf {D} _{3}\\\end{pmatrix}}{\begin{pmatrix}\mathbf {I} &\mathbf {L} _{21}^{\mathrm {T} }&\mathbf {L} _{31}^{\mathrm {T} }\\0&\mathbf {I} &\mathbf {L} _{32}^{\mathrm {T} }\\0&0&\mathbf {I} \\\end{pmatrix}}\\&={\begin{pmatrix}\mathbf {D} _{1}&&(\mathrm {symmetric} )\\\mathbf {L} _{21}\mathbf {D} _{1}&\mathbf {L} _{21}\mathbf {D} _{1}\mathbf {L} _{21}^{\mathrm {T} }+\mathbf {D} _{2}&\\\mathbf {L} _{31}\mathbf {D} _{1}&\mathbf {L} _{31}\mathbf {D} _{1}\mathbf {L} _{21}^{\mathrm {T} }+\mathbf {L} _{32}\mathbf {D} _{2}&\mathbf {L} _{31}\mathbf {D} _{1}\mathbf {L} _{31}^{\mathrm {T} }+\mathbf {L} _{32}\mathbf {D} _{2}\mathbf {L} _{32}^{\mathrm {T} }+\mathbf {D} _{3}\end{pmatrix}}\end{aligned}}}
在此,矩阵中每一个元素都是一个子矩阵。故此,可依照上述递归公式模拟计算:
D
j
=
A
j
j
−
∑
k
=
1
j
−
1
L
j
k
D
k
L
j
k
T
{\displaystyle \mathbf {D} _{j}=\mathbf {A} _{jj}-\sum _{k=1}^{j-1}\mathbf {L} _{jk}\mathbf {D} _{k}\mathbf {L} _{jk}^{\mathrm {T} }}
L
i
j
=
(
A
i
j
−
∑
k
=
1
j
−
1
L
i
k
D
k
L
j
k
T
)
D
j
−
1
{\displaystyle \mathbf {L} _{ij}=\left(\mathbf {A} _{ij}-\sum _{k=1}^{j-1}\mathbf {L} _{ik}\mathbf {D} _{k}\mathbf {L} _{jk}^{\mathrm {T} }\right)\mathbf {D} _{j}^{-1}}
上述计算涉及矩阵相乘同精确的求逆,故而实践中不能使用过大的分块矩阵。
修正 分解
在实践中经常有修正 科列斯基分解的需求。即,经已计算一些矩阵
A
{\displaystyle \mathbf {A} }
的科列斯基分解
A
=
L
L
∗
{\displaystyle \mathbf {A} =\mathbf {LL} ^{*}}
,然后在
A
{\displaystyle \mathbf {A} }
上稍作修改而产生的矩阵
A
~
{\displaystyle {\tilde {\mathbf {A} }}}
,欲对其进行一个修正的科列斯基分解
A
~
=
L
~
L
~
∗
{\displaystyle {\tilde {\mathbf {A} }}={\tilde {\mathbf {L} }}{\tilde {\mathbf {L} }}^{*}}
。问题是,能否用已知的
A
{\displaystyle \mathbf {A} }
的分解去修正
A
~
{\displaystyle {\tilde {\mathbf {A} }}}
的分解。
秩为1的修正(相加)
如果修正 矩阵
A
~
=
A
+
x
x
∗
{\displaystyle {\tilde {\mathbf {A} }}=\mathbf {A} +\mathbf {xx} ^{*}}
,则其修正的分解被称为秩为1的修正(相加)。。
以下是一个 MATLAB 语言写的实现秩为1的修正的小程序:
function [L] = cholupdate ( L,x)
n = length ( x );
for k = 1 : n
r = sqrt ( L ( k , k ) ^2 + x ( k ) ^2 );
c = r / L ( k , k );
s = x ( k ) / L ( k , k );
L ( k , k ) = r ;
L ( k + 1 : n , k ) = ( L ( k + 1 : n , k ) + s * x ( k + 1 : n )) / c ;
x ( k + 1 : n ) = c * x ( k + 1 : n ) - s * L ( k + 1 : n , k );
end
end
秩为1的修正(相减)
如果
A
~
=
A
−
x
x
∗
{\displaystyle {\tilde {\mathbf {A} }}=\mathbf {A} -\mathbf {x} \mathbf {x} ^{*}}
,那么只有当
A
~
{\displaystyle {\tilde {\mathbf {A} }}}
仍然是正定的时候该方法才适用。
此条件下,上面的代码也可用于相减的情况,只需要将r
和L(k+1:n,k)
的赋值式的加号替换成减号。
证明半正定矩阵的科列斯基分解唯一
上面的算法表明,每一个正定矩阵
A
{\displaystyle A}
都有一个科列斯基分解。此结论可以加入一些限制条件 后,推广到半正定的情况。但该条件 并未被完全建立,例如,它未给出明确的数值算法以计算科列斯基因子。
如果
A
{\displaystyle A}
是一个
n
×
n
{\displaystyle n\times n}
的半正定矩阵,则序列
{
A
}
=
{
A
+
(
1
/
k
)
I
n
}
{\displaystyle \{\mathbf {A} \}=\{\mathbf {A} +(1/k)\mathbf {I} _{n}\}}
是一个由正定矩阵构成的序列。而且,在算子范数 中
A
k
→
A
{\displaystyle \mathbf {A} _{k}\rightarrow \mathbf {A} }
在正定的情形下,每一个
A
k
{\displaystyle \mathbf {A} _{k}}
都有其科列斯基分解
A
k
=
L
k
L
k
∗
{\displaystyle \mathbf {A} _{k}=\mathbf {L} _{k}\mathbf {L} _{k}^{*}}
。根据算子范数的性质,
‖
L
k
‖
2
≤
‖
L
k
L
k
∗
‖
=
‖
A
k
‖
{\displaystyle \|\mathbf {L} _{k}\|^{2}\leq \|\mathbf {L} _{k}\mathbf {L} _{k}^{*}\|=\|\mathbf {A} _{k}\|}
因而
{
L
k
}
{\displaystyle \{\mathbf {L} _{k}\}}
是算子巴拿赫空间 上的一个有界集 ,因此是相对紧空间 (因为基础向量空间是有限维的)。因此,它有一个带有限制
L
{\displaystyle \mathbf {L} }
收敛子序列,也记作
{
L
k
}
{\displaystyle \{\mathbf {L} _{k}\}}
。容易验证,矩阵
L
{\displaystyle \mathbf {L} }
具有所需的特性,例如,满足
A
=
L
L
∗
{\displaystyle \mathbf {A} =\mathbf {LL} ^{*}}
以及
L
{\displaystyle \mathbf {L} }
是下三角矩阵且其对角元素非负。对于所有的
x
{\displaystyle x}
和
y
{\displaystyle y}
都有,
⟨
A
x
,
y
⟩
=
⟨
lim
A
k
x
,
y
⟩
=
⟨
lim
L
k
L
k
∗
x
,
y
⟩
=
⟨
L
L
∗
x
,
y
⟩
{\displaystyle \langle \mathbf {A} x,y\rangle =\langle \lim \mathbf {A} _{k}x,y\rangle =\langle \lim \mathbf {L} _{k}\mathbf {L} _{k}^{*}x,y\rangle =\langle \mathbf {L} \mathbf {L} ^{*}x,y\rangle }
因此,
A
=
L
L
∗
{\displaystyle \mathbf {A} =\mathbf {LL} ^{*}}
。因为基础向量空间是有限维的,所以算子空间的所有拓扑结构都是等价的。故此,范数上
L
k
{\displaystyle \mathbf {L} _{k}}
收敛于
L
{\displaystyle \mathbf {L} }
,意味着,
L
k
{\displaystyle \mathbf {L} _{k}}
的每个元素都独立地收敛于
L
{\displaystyle \mathbf {L} }
。此同时暗示,由于每个
L
k
{\displaystyle \mathbf {L} _{k}}
都是对角元素非负的下三角矩阵,
L
{\displaystyle \mathbf {L} }
亦如此。
推广
科列斯基分解可以推广到元素中包含算子 的矩阵 (称为算子矩阵)。[来源请求] 令
{
H
n
}
{\displaystyle \{{\mathcal {H}}_{n}\}}
是希尔伯特空间 上的一个序列。考虑如下算子矩阵
A
=
[
A
11
A
12
A
13
A
12
∗
A
22
A
23
A
13
∗
A
23
∗
A
33
⋱
]
{\displaystyle \mathbf {A} ={\begin{bmatrix}\mathbf {A} _{11}&\mathbf {A} _{12}&\mathbf {A} _{13}&\;\\\mathbf {A} _{12}^{*}&\mathbf {A} _{22}&\mathbf {A} _{23}&\;\\\mathbf {A} _{13}^{*}&\mathbf {A} _{23}^{*}&\mathbf {A} _{33}&\;\\\;&\;&\;&\ddots \end{bmatrix}}}
满足直和
H
=
⊕
n
H
n
,
{\displaystyle {\mathcal {H}}=\oplus _{n}{\mathcal {H}}_{n},}
其中每一个
A
i
j
:
H
j
→
H
i
{\displaystyle \mathbf {A} _{ij}:{\mathcal {H}}_{j}\rightarrow {\mathcal {H}}_{i}}
都是一个有界算子。如果
A
{\displaystyle \mathbf {A} }
是正定或半正定的,即对于有限的
k
{\displaystyle k}
及任意的
h
∈
⊕
n
=
1
k
H
k
,
{\displaystyle h\in \oplus _{n=1}^{k}{\mathcal {H}}_{k},}
都有
⟨
h
,
A
h
⟩
≥
0
{\displaystyle \langle h,\mathbf {A} h\rangle \geq 0}
,则存在一个下三角的算子矩阵
L
{\displaystyle \mathbf {L} }
使得
A
=
L
L
∗
{\displaystyle \mathbf {A} =\mathbf {LL} ^{*}}
。此外也可以把
L
{\displaystyle \mathbf {L} }
的对角元素化为正数。
用编程语言实现
参见
脚注
^ Press, William H.; Saul A. Teukolsky; William T. Vetterling; Brian P. Flannery. Numerical Recipes in C: The Art of Scientific Computing (second edition) . Cambridge University England EPress. 1992: 994 [2017-08-30 ] . ISBN 0-521-43108-5 . (原始内容 存档于2018-08-25).
^ Golub & Van Loan (1996 , p. 143), Horn & Johnson (1985 , p. 407), Trefethen & Bau (1997 , p. 174)
^ Golub & Van Loan (1996 , p. 147)
^ Horn & Johnson (1985 , p. 407)
^ Krishnamoorthy, Aravindh; Menon, Deepak. Matrix Inversion Using Cholesky Decomposition 1111 : 4144. 2011. Bibcode:2011arXiv1111.4144K . arXiv:1111.4144 .
^ Fang, Haw-ren; O’Leary, Dianne P. Modified Cholesky Algorithms: A Catalog with New Approaches (PDF) . 8 August 2006 [2017-08-30 ] . (原始内容存档 (PDF) 于2017-05-16).
^ Nocedal, Jorge. Numerical Optimization. Springer. 2000.
^ Fang, Haw-ren. Analysis of Block LDLT Factorizations for Symmetric Indefinite Matrices. 24 August 2007.
参考文献
Dereniowski, Dariusz; Kubale, Marek. 5th International Conference on Parallel Processing and Applied Mathematics (PDF) . Lecture Notes on Computer Science 3019 . Springer-Verlag: 985–992. 2004 [2017-08-30 ] . ISBN 978-3-540-21946-0 . doi:10.1007/978-3-540-24669-5_127 . (原始内容 (PDF) 存档于2011-07-16). .
Golub, Gene H. ; Van Loan, Charles F. Matrix Computations 3rd. Baltimore: Johns Hopkins. 1996. ISBN 978-0-8018-5414-9 .
Horn, Roger A.; Johnson, Charles R. Matrix Analysis. Cambridge University Press. 1985. ISBN 0-521-38632-2 . .
S. J. Julier and J. K. Uhlmann. "A General Method for Approximating Nonlinear Transformations of ProbabilityDistributions".
S. J. Julier and J.K. Uhlmann, "A new extension of the Kalman filter to nonlinear systems," in Proc. AeroSense: 11th Int. Symp. Aerospace/Defence Sensing, Simulation and Controls, 1997, pp. 182–193.
Trefethen, Lloyd N. ; Bau, David. Numerical linear algebra . Philadelphia: Society for Industrial and Applied Mathematics. 1997. ISBN 978-0-89871-361-9 . .
外部链接
科学史
Sur la résolution numérique des systèmes d'équations linéaires , Cholesky's 1910 manuscript, online and analyzed on BibNum (页面存档备份 ,存于互联网档案馆 ) (法文) (英文) [for English, click 'A télécharger']
资讯
电脑代码
模拟中矩阵的利用
在线电脑