线性方程求解

这篇文章介绍了解线性方程的思路，包括 LU 分解、Cholesky 分解和 Iterative Methods。

计算经济学的基础问题是解线性方程，一般表示为 $$A x = b.$$ 其中 $A$ 是一个 $n \times n$ 的矩阵，$b$ 是 $n$ 维向量，两者给定求 $n$ 维向量 $x$ 。

Forward subsititution 链接到标题

如果 $A$ 是一个 lower trangular matrix, $$ A = \begin{bmatrix} a_{11} & 0 & 0 & \cdots & 0 \newline a_{21} & a_{22} & 0 & \cdots & 0 \newline a_{31} & a_{32} & a_{33} & \cdots & 0 \newline \vdots & \vdots & \vdots & \ddots & \vdots \newline a_{n1} & a_{n2} & a_{n3} & \cdots & a_{nn} \end{bmatrix}. $$

那么 $x$ 的各元素可以表示为 $$ x_i = \frac{b_i - \sum_{j=1}^{i-1}a_{ij}x_j}{a_{ii}}, \qquad\forall i. $$ 这一算法可以通过代码实现，如 Python。 类似地，如果 $A$ 是一个 upper triangular matrix，可以通过 backward subsititution 求解， 相关算法可以参考 Wiki。

LU 分解 链接到标题

但不是所有矩阵都是 triangular 矩阵，这种情况下可以使用 L-U 分解算法。

在分解阶段，利用高斯消除法 (Gaussian elimination) 把矩阵 $A$ 分解为两个三角矩阵的积 $$ A = LU, $$ 其中 $L$ 是 (row-permuted) lower triangular matrix，$U$ 是 upper triangular matrix。这样我们有 $$ A x = (LU) x = L (Ux) = b. $$

在求解阶段，先利用 forward subsititution 求 $y$, $$ L y = b. $$ 然后利用 backward subsititution 求 $x$, $$ U x = y. $$

Cholesky 分解 链接到标题

当 $A$ 是对称正定矩阵 (symetric positive definite) 时， 另一种分解方式——Cholesky 分解的速度更快， 需要的运算一般比 LU 分解少一半。

因为 $A$ 对称正定矩阵，因此 $A$ 可以表示为一个 upper triangular 矩阵和其转置的积 $$ A = U^T U. $$ 这里的 $U$ 叫做 $A$ 的 Cholesky factor 或 平方根。

同样利用 forward subsititution 和 backward subsititution 可以解出 $x$ $$ \begin{gather*} A x = U^T U x = U^T(Ux) = b \newline U^T y = b \newline U x = y \end{gather*} $$

在实际应用中，matlab 和 julia 的 \ 运算符能自动检测 $A$ 的 特性，并应用最合适的分解方法。

迭代法 Iterative Methods 链接到标题

LU 分解和 Cholesky 分解都能得出精确解，但在应对大型方程组的时候效率较低，另一种近似的方法能牺牲精确性，但获得较快的运算速度。

一个常用的近似方法是 iterative methods，其原理是找到一个可逆矩阵 $Q$ 满足 $$ Q x = b + (Q - A) x $$ 或者 $$ x = Q^{-1} b + (I - Q^{-1} A) x. $$

由此我们可以得到一个迭代的规则 $$ x^{(k+1)} \leftarrow Q^{-1} b + (I - Q^{-1} A) x^{(k)} $$ 如果该式能够收敛，则会收敛到我们想要得到的方程解。

两种最常用的迭代法分别是 Gauss-jacobi 和 Gauss-Seidel 法，其中 Gauss-Jacobi 法把 $Q$ 设定为 $A$ 的对角矩阵，Gauss-Seidel 法则把 $Q$ 定为 $A$ 的上三角矩阵。

下面我给出两种方法的 julia 实现。

# Gauss-Jacobi
function gaussjacobi(A, b, init, tol = 1e-10)
    d = diag(A)
    dist = 100
    x = init
    while dist >= tol
        dx = (b - A*x) ./ d
        x = x + dx
        dist = norm(dx)
    end
    return x
end

# Gauss-Seidel
function gaussseidel(A, b, init, tol = 1e-10, λ = 1.5)
    Q = tril(A)
    dist = 100
    x = init
    while dist >= tol
        dx = Q \ (b - A*x)
        x = x + λ * dx
        dist = norm(dx)
    end
    return x
end

更多的细节可以参考 Miranda and Fackler 的 Applied Computational Economics and Finance 第二章。