没看第一章的人可以点击这里，里面介绍了些数学符号和矩阵乘法。
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}
其中 K 是 A[2:end, 2:end]， w 是A[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
如果代码看不懂的话，可以在下面提问。谢谢！