线性代数教程2: Cholesky分解


#1

没看第一章的人可以点击这里,里面介绍了些数学符号和矩阵乘法。

Cholesky分解

知道了矩阵乘法以后我们就可以开始学习矩阵分解。一般来说在数值线性代数里,解线性方程组都是通过分解矩阵来做到的。即把矩阵分解成好解的矩阵,例如酉矩阵或者三角矩阵。

解一个非常一般的线性系统是比较复杂的,所以我们先来介绍一个最简单特殊情况。如果矩阵是对称而且正定的话,我们可以用一个非常简单的递归关系来算出Cholesky分解。

假设我们有一个任意的对称正定矩阵 A ,那么我们可以把它写成

\begin{aligned} A&=\begin{bmatrix} a_{11} & w^\dagger \\ w & K \end{bmatrix}\\ &= \begin{bmatrix} \alpha & 0 \\ \frac{w}{\alpha} & I \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0 & K - \frac{ww^\dagger}{a_{11}} \end{bmatrix} \begin{bmatrix} \alpha & \frac{w^\dagger}{\alpha} \\ 0 & I \end{bmatrix} \\ &= L_1 A_1 L_1^\dagger \end{aligned}

其中 KA[2:end, 2:end]wA[2:end, 1],当然 w^\dagger 就是 A[2:end, 1]' 了。最后, \alpha\sqrt{a_{11}} 。我们可以看出这里的 A_1 还是一个正定并且对称的矩阵,并且 L_1 是一个下三角矩阵。这样,我们就得到了一个递归的关系。我们可以像剥洋葱那样把 A 分解成 A=LL^\dagger ,即递归的计算 L 的对角线下的第 i 列。

A = \underbrace{L_1, L_2 ,\cdots L_m}_{L}\underbrace{L_m^\dagger ,\cdots L_2^\dagger L_1^\dagger}_{L^\dagger}

这种写法利用了下三角矩阵的如下性质:

julia> L₁ = [1 0 0; 2 1 0; 3 0 1]
3×3 Array{Int64,2}:
 1  0  0
 2  1  0
 3  0  1

julia> L₂ = [1 0 0; 0 1 0; 0 4 1]
3×3 Array{Int64,2}:
 1  0  0
 0  1  0
 0  4  1

julia> L₁ * L₂
3×3 Array{Int64,2}:
 1  0  0
 2  1  0
 3  4  1

最后,Cholesky 分解的 Julia 代码就是:

using LinearAlgebra
function chol_!(A::AbstractMatrix{T}) where T
    @views @inbounds begin
        m, n = size(A)
        @assert m==n
        α = A[1, 1] = sqrt(A[1, 1])
        if n > 1
            w = rmul!(A[2:n, 1], inv(α))
            K = A[2:n, 2:n]
            for j in eachindex(w), i in j:n-1
                K[i, j] -= w[i]*w[j]'
            end
            chol_!(K)
        end
        return LowerTriangular(A)
    end
end
chol_(A) = chol_!(copy(A))

using Test
@testset "Cholesky Factorization Tests" begin
    for siz in 2:15:100, T in (Float64, ComplexF64)
        A = randn(T, siz, siz)
        A = A'A # 对称正定矩阵
        L = chol_(A)
        @test L*L' ≈ A
    end
end

如果代码看不懂的话,可以在下面提问。谢谢!


#2

有一些数学上的细节(比如说下三角矩阵乘法的特性)和背景知识我还没有写。之后我和 @Gnimuc 会找时间补上的。


#3

学习学习,多谢分享


#4

支持一个。我觉得楼主讲的非常清楚。未来能不能做个finite difference 和 finite element的专题。不好意思扯远了


#5

那还是得看时间了。如果时间允许的话,我可以介绍一下解PDE的基础。