线性方程组与本征值问题(02)：Cholesky 分解

Cholesky 分解基于这样一个事实: 对于一个正定自伴矩阵\(\mathbf{A}\in \mathbb{C}^{n\times n}\), 可以找到一个矩阵\(\mathbf{H}\)使得\(\mathbf{A}=\mathbf{H}^\dagger\mathbf{H}\). 事实上, 可以要求\(\mathbf{H}\)是一个上三角矩阵, 则这个分解称为 Cholesky 分解. Cholesky 分解同样把复杂的线性方程组转化为两个较简单的方程组: \[ \mathbf{A}\pmb{x}=\pmb{b}\Rightarrow\begin{cases}\displaystyle \mathbf{H}^\dagger\pmb{y} = \pmb{b} \\\\\displaystyle \mathbf{H}\pmb{x}=\pmb{y} \end{cases} \] 一个显而易见的引理是, 正定自伴矩阵的主子矩阵仍是正定自伴的. 这可以直接在如下求和式中 \[ \sum_{j=1}^n\sum_{k=1}^na_{jk}x_jx_k \] 将\(j,k=m+1,m+2,\cdots,n\)的\(x_j,x_k\)取作\(0\), 不应影响正定二次型的结论, 从而得证.

有了这个引理, 对 Cholesky 分解的研究仍可以按照归纳法的思路进行. 假设主子矩阵\(\mathbf{A}_{m}\) 存在分解矩阵\(\mathbf{H}_m\), 更高阶的主子矩阵可以进行分块 \[ \mathbf{A}_{m+1}=\begin{pmatrix} \mathbf{A}_m & \pmb{v} \\ \pmb{v}^\dagger & \lambda \end{pmatrix} \] 我们希望它仍存在分解矩阵\(\mathbf{H}_{m+1}\), 同样可以分块为 \[ \mathbf{H}_{m+1}=\begin{pmatrix} \mathbf{H}_m & \pmb{u} \\ \pmb{0}^\dagger & \mu \end{pmatrix} \\ \mathbf{H}_{m+1}^\dagger=\begin{pmatrix} \mathbf{H}_m^\dagger & \pmb{0} \\ \pmb{u}^\dagger & \mu \end{pmatrix} \] 则有 \[ \mathbf{H}_m^\dagger\pmb{u}=\pmb{v} \\ \pmb{u}^\dagger\pmb{u}+\mu^2=\lambda \] 设\(\mathbf{H}^\dagger\)的\((j,k)\)元为\(h_{jk}\), 则 \[ \sum_{l=1}^{k}h_{jl}h_{kl}^\dagger=a_{jk} \\(k=1,2,\cdots,j) \] 写成\(h_{jk}\)的显式 \[ \begin{align} h_{jk}&=\frac{\displaystyle a_{jk}-\sum_{l=1}^{k-1}h_{jl}h^\dagger_{kl}}{h^\dagger_{kk}} \newline(k&=1,2,\cdots,j-1) \newline h_{jj}&=\sqrt{a_{jj}-\sum_{l=1}^{j-1}|h_{jl}|^2} \end{align} \]

算法总结

如下为算法流程:

import numpy as np # 仅作为示意, 指标起点取0
h[0,0] = np.sqrt(a[0,0])
for j in range(1, n):
    for k in range(j):
        h[j,k]=(a[j, k]-h[j, :k]@np.conj(h[k, :k])) / h[j, j]
       
    h[j, j] = np.sqrt(a[j, j] - np.sum(np.abs(h[j, :j]) ** 2))

Cholesky 分解的计算量约为\(\displaystyle\frac{n^3}{3}\), 比 LU 分解通常少一半. 它的特殊形式为正定实对称矩阵的分解, 这种矩阵相当常见, 比如协方差矩阵.