用Python的Numpy求解线性方程组

在本文中，您将看到如何使用Python的Numpy库解决线性方程组。

由Kaizong Ye，Liao Bao撰写

维基百科将线性方程组定义为：

可下载资源

在数学中，线性方程组（或线性系统）是两个或多个涉及同一组变量的线性方程的集合。

预备工作：线性方程组类型的判断

对于 $Ax=b$ ，(a) $\operatorname{rank}[A]=\operatorname{rank}[A,b]=n$ ，则解唯一，称为相容方程组；(b) $\operatorname{rank}[A]=\operatorname{rank}[A,b]<n$ ，则解无穷，称为超定方程组；(c) $\operatorname{rank}[A]<\operatorname{rank}[A,b]$ ，则无解，称为不相容方程组。

此外，判断方程组是否有解还可以用以下2个充要条件：(a) 如果 $A\in \mathbb{C}^{m \times n}$ 是左可逆矩阵， $A_L^{-1}$ 是它的一个左逆矩阵，若 $(E_m-AA_L^{-1})b=0$ 成立，则 $Ax=b$ 有唯一解，且解为 $x=(A^HA)^{-1}A^Hb$ ；(b) 如果 $A\in \mathbb{C}^{m \times n}$ 是右可逆矩阵，则 $Ax=b$ 对任何 $b \in \mathbb{C}^m$ 都有解，若 $b \neq 0$ 时，则方程组的解可表示成 $x=A_R^{-1}b$ ，其中 $A_R^{-1}$ 是 $A$ 的一个右逆矩阵。（单边逆的定义、判别和求法，参见【1】P178-182）

首先讨论相容线性方程组的解法。

一. 三角分解法

思路： $A$ 被三角分解成 $UR$ 后，将 $Ax=b$ 转化为 $Rx=U^Hb$ ，算出 $U^Hb$ 。因为 $R$ 是上三角矩阵，容易求解出最后一行的 $x_n$ ；将 $x_n$ 代入倒数第二行，容易求解出 $x_{n-1}$ ；将 $x_n$ , $x_{n-1}$ 代入倒数第三行，容易求解出 $x_{n-2}$ ，以此类推。

1.1 QR分解

1.1.1 形式一: 满秩方阵

(a) 对 $A$ 进行列分块，得列分块向量 $(\alpha_1,\alpha_2,...,\alpha_n)$

(b) 对列分块向量施密特正交化，得到正交向量 $(\beta_1,\beta_2,...,\beta_n)$ ，其中

$\begin{equation} \beta_1=\frac{\alpha_1}{||\alpha_1||},\beta_2=\frac{\alpha_2-(\alpha_2,\beta_1)\beta_1}{||\alpha_2-(\alpha_2,\beta_1)\beta_1||},\beta_3=\frac{\alpha_3-(\alpha_3,\beta_1)\beta_1-(\alpha_3,\beta_2)\beta_2}{||\alpha_3-(\alpha_3,\beta_1)\beta_1-(\alpha_3,\beta_2)\beta_2||},...,\beta_i=\frac{\alpha_i-\sum_{j=1}^{i-1}(\alpha_i,\beta_j)\beta_j}{||\alpha_i-\sum_{j=1}^{i-1}(\alpha_i,\beta_j)\beta_j||} \end{equation}$

(c) 令系数 $\begin{aligned} \left\{ \begin{array}{1} k_{ij}=(\alpha_i,\beta_j) \quad i=2,3,...,n;j=1,2,...,i-1\\ k_{ii}=\left\{\begin{array}{1}||a_{ij}|| \quad i=1 \\ ||a_i-\sum_{j=1}^{i-1}(\alpha_i,\beta_j)\beta_j|| \quad i=2,3,...,n \end{array}\right. \end{array} \right. \end{aligned}$

(d) 最后得到 $A=(k_{11}\beta_1,k_{21}\beta_1+k_{22}\beta_2,...,k_{n1}\beta_1+k_{n2}\beta_2+...+k_{nn}\beta_n)=(\beta_1,\beta_2,...,\beta_n) \begin{pmatrix} k_{11} &k_{21}&...&k_{n1}\\ 0 &k_{22} &...&k_{n2}\\ &...\\ 0 &0 &... &k_{nn} \end{pmatrix}$

1. 拓展：上述过程中 $A \in \mathbb{R}_n^{n \times n}$ ，同时还可进行 $A=$ $LQ$ 分解，在步骤(a)中行分块即可。类似地， $A \in \mathbb{C}_n^{n \times n}$ 时可进行 $A=UR$ 或 $A=LU$ 分解。区别在于 $Q$ 是正交矩阵，而 $U$ 是酉矩阵。下面开始仅讨论复矩阵的情形，相应的实矩阵情形很好类推。

1.1.2 形式二: 行满秩矩阵

$A=[L\,O]U,A \in \mathbb{C}_m^{m \times n}$ ，其中 $L$ 为 $m$ 阶正线下三角复矩阵， $U$ 为 $n$ 阶酉矩阵， $O$ 为零矩阵。（【1】P92证明）

1.1.3 形式三: 列满秩矩阵

$A=U \begin{bmatrix}R\\O \end{bmatrix},A \in \mathbb{C}_n^{m \times n}$ ，其中 $U$ 为 $m$ 阶酉矩阵， $R$ 为 $n$ 阶正线上三角复矩阵， $O$ 为零矩阵。（【1】P92证明）

1.1.4 形式四: 任意矩阵

$A=U\begin{bmatrix}L&O\\O&O\end{bmatrix} V,A\in \mathbb{C}_r^{m \times n}$ ，其中 $U$ 为 $m$ 阶酉矩阵， $V$ 为 $n$ 阶酉矩阵， $L$ 为 $r$ 阶正线下三角复矩阵， $O$ 为零矩阵。还可以分解为 $A=U\begin{bmatrix}R&O\\O&O\end{bmatrix} V$ ， $R$ 为 $r$ 阶正线下三角复矩阵。（【1】P94证明）

1.2 LU分解

1.2.1 形式一: 满秩方阵

$A=LR^*=L^*R=L^*DR^*,A\in \mathbb{C}_n^{n \times n}$ ， $A$ 的各阶顺序主子式行列式不为零， $L,R,D$ 的主对角线上元素不为0。（【1】P89-92证明，4个命题等价）

1. 参数含义： $R$ 为上三角复矩阵， $L$ 为下三角复矩阵， $R^*$ 为单位上三角复矩阵， $L^*$ 为单位下三角复矩阵， $D$ 为对角矩阵。
2. 拓展：同样地， $A\in \mathbb{R}_n^{n \times n}$ 的分解形式类推，相应的所得则为实矩阵。
3. 分解方法：
(1) 待定系数法
(a) 写出含待定系数的分解形式

$\begin{pmatrix} 2&4&4&2 \\3&3&12&6\\2&4&-1&2\\4&2&1&1 \end{pmatrix}= \begin{pmatrix} 1\\l_{21}&1\\l_{31}&l_{32}&1\\l_{41} &l_{42}&l_{43}&1 \end{pmatrix} \begin{pmatrix} u_{11}&u_{12}&u_{13}&u_{14}\\&u_{22}&u_{23}&u_{24}\\&&u_{33}&u_{34}\\&&&u_{44} \end{pmatrix}$
(b) 按照规则“ 先行后列”、“列除行不除”、“旧元素减去其所在行和列前(k-1)个元素的对应乘积然后求和”，对照上式一步步地写出待定系数

$\begin{align} \begin{pmatrix} 2&4&4&2 \\3&3&12&6\\2&4&-1&2\\4&2&1&1 \end{pmatrix}&= \begin{pmatrix} 1\\3/2&1\\1&l_{32}&1\\2 &l_{42}&l_{43}&1 \end{pmatrix} \begin{pmatrix} 2&4&4&2\\&u_{22}&u_{23}&u_{24}\\&&u_{33}&u_{34}\\&&&u_{44} \end{pmatrix}=\begin{pmatrix} 1\\3/2&1\\1&0&1\\2&-2&l_{43}&1 \end{pmatrix} \begin{pmatrix} 2&4&4&2\\&-3&6&3\\&&-5&0\\&&&u_{44} \end{pmatrix}\\&=\begin{pmatrix} 1\\3/2&1\\1&0&1\\2&2&19/5&1 \end{pmatrix} \begin{pmatrix} 2&4&4&2\\&-3&6&3\\&&-5&0\\&&&-9 \end{pmatrix} \end{align}$
(2) 高斯消元法
(a) 搭建框架
$\begin{pmatrix}1\\&1\\&&1\\&&&1\end{pmatrix}\begin{pmatrix}2&4&4&2\\3&3&12&6\\2&4&-1&2\\4&2&1&1\end{pmatrix}=\begin{pmatrix}&\\&\\&\\&&&&&\end{pmatrix}$
(b) 提取高斯消元（即初等行变换）过程中的乘子和结果。(注意行列的对应，第 $i$ 行对第 $j$ 行做初等行变换的乘子，对应到结果矩阵的 $A_{ji}$ 元素）

$\begin{pmatrix}1\\-3/2&1\\-1&&1\\-2&&&1\end{pmatrix}\begin{pmatrix}2&4&4&2\\3&3&12&6\\2&4&-1&2\\4&2&1&1\end{pmatrix}=\begin{pmatrix}2&4&4&2\\0&-3&6&3\\0&0&-5&0\\0&-6&-7&-3\end{pmatrix}\Rightarrow\begin{pmatrix}1\\-3/2&1\\-1&0&1\\-2&-2&&1\end{pmatrix}\begin{pmatrix}2&4&4&2\\3&3&12&6\\2&4&-1&2\\4&2&1&1\end{pmatrix}=\begin{pmatrix}2&4&4&2\\0&-3&6&3\\0&0&-5&0\\0&0&-19&-9\end{pmatrix}\\ \Rightarrow\begin{pmatrix}1\\-3/2&1\\-1&0&1\\-2&-2&-19/5&1\end{pmatrix}\begin{pmatrix}2&4&4&2\\3&3&12&6\\2&4&-1&2\\4&2&1&1\end{pmatrix}=\begin{pmatrix}2&4&4&2\\0&-3&6&3\\0&0&-5&0\\0&0&0&-9\end{pmatrix}\Rightarrow\begin{pmatrix}2&4&4&2\\3&3&12&6\\2&4&-1&2\\4&2&1&1\end{pmatrix}=\begin{pmatrix}1\\3/2&1\\1&0&1\\2&2&19/5&1\end{pmatrix}\begin{pmatrix}2&4&4&2\\0&-3&6&3\\0&0&-5&0\\0&0&0&-9\end{pmatrix}\\$
4. 注意：
(a) 利用高斯消元法对矩阵 $A$ 进行 $LU$ 矩阵分解时，绝对值小的主元可能产生麻烦，因此一般在进行前先判断是否对行进行调整。假设调整后的矩阵是 $B$ ,分解为 $B=L^*U$ ，那么， $A$ 的分解又该表示为什么形式呢？我们用一个单位矩阵 $I$ 来表征这种变换，如果要交换矩阵 $A$ 中的一些行，那我们就对单位矩阵 $I$ 中相应的行进行交换，得到矩阵 $P$ 。因此有 $PA=L^*U\Rightarrow A=P^{-1}L^*U$ 。
(b) 得到 $A=L^*U$ 后，我们可以很容易地基于此，将 $U$ 拆成 $DU^*$ 的形式从而得到 $A=L^*DU^*$ ；再将 $L^*D$ 合并为 $L$ 从而得到 $A=LU^*$ ，过程如下：

$\begin{align} \begin{pmatrix} 2&4&4&2 \\3&3&12&6\\2&4&-1&2\\4&2&1&1 \end{pmatrix} &=\begin{pmatrix} 1\\3/2&1\\1&0&1\\2&2&19/5&1 \end{pmatrix} \begin{pmatrix} 2&4&4&2\\&-3&6&3\\&&-5&0\\&&&-9 \end{pmatrix} =\begin{pmatrix} 1\\3/2&1\\1&0&1\\2&2&19/5&1 \end{pmatrix} \begin{pmatrix} 2&&&\\&-3&&\\&&-5&\\&&&-9 \end{pmatrix}\begin{pmatrix}1&2&2&1\\&1&-2&-1\\&&1&0\\&&&1\end{pmatrix} \\ &=\begin{pmatrix} 2\\3&-3\\2&0&-5\\4&-6&-19&-9\end{pmatrix}\begin{pmatrix}1&2&2&1\\&1&-2&-1\\&&1&0\\&&&1\end{pmatrix} \end{align}$
(c) 对于三对角矩阵，如果它的非零系数在主对角线上和两条次对角线上，那么可以对其进行 $LU$ 分解，上述的两种分解方法都可以使用。

1.2.2 形式二: 行满秩矩阵

$A=LU,A\in \mathbb{C}_m^{m \times n}$ ，其中 $L$ 是 $m$ 阶正线下三角矩阵， $U\in U_m^{m\times n }$ ，表示以 $m$ 个两两正交的单位向量为行组成的矩阵的集合。（【1】P93证明，证明过程给出了分解步骤）

1.2.3 形式三: 列满秩矩阵

$A=UR,A\in \mathbb{C}_n^{m \times n}$ ，其中， $R$ 是 $n$ 阶正线上三角矩阵， $U\in U_n^{m\times n }$ ，表示以 $n$ 个两两正交的单位向量为列组成的矩阵的集合。（【1】P93证明，证明过程给出了分解步骤）

1.3 Cholesky分解

$A=R^HR,A\in \mathbb{C}_n^{n \times n}$ ， $A$ 是正定Hermite矩阵， $R$ 是正线上三角复矩阵。类似的， $A \in\mathbb{R}_n^{n \times n}$ 则要求 $A$ 为实对称正定矩阵。

1. 拓展：对于 $A$ 同时还是正定Hermite矩阵的情况，它还可以分解为 $A=R^{*H}DR^*$ ，这里的 $R^*$ 为单位上三角复矩阵， $D$ 为对角阵。
2. 分解方法（即待定系数法）：
$\begin{pmatrix}4&-2&2\\-2&2&-4\\2&-4&11\end{pmatrix}=\begin{pmatrix}l_{11}\\l_{12}&l_{22}\\l_{13}&l_{23}&l_{33}\end{pmatrix}\begin{pmatrix}l_{11}&l_{12}&l_{13}\\&l_{22}&l_{23}\\&&l_{33}\end{pmatrix}=\begin{pmatrix}l_{11}^2&l_{11}l_{12}&l_{11}l_{13}\\l_{11}l_{12}&l_{12}^2+l_{22}^2&l_{12}l_{13}+l_{22}l_{23}\\l_{11}l_{13}&l_{12}l_{13}+l_{22}l_{23}&l_{13}^2+l_{23}^2+l_{33}^2\end{pmatrix}$

二. 迭代解法

2.1 迭代思想

从最简单的一元函数零点问题开始考虑，假设 $f(x)$ 具有一个零点 $x^*$ 使得 $f(x^*)=0$ 成立，通过变换我们可以得到等式 $x^*=\varphi(x^*)$ 。通常情况下，我们可以直接把 $f(x)=0$ 变换为 $x=\varphi(x)$ 的形式，容易发现，仅当 $x=x^*$ 时该形式才成立，这时我们把通过 $f(x)=0$ 求解 $x^*$ 的问题转换成了通过 $x=\varphi(x)$ 求解 $x^*$ 的问题，显然后者更容易求解。为什么呢？我们选择一个适当的初始值 $x_0$ 代入到等式右边，可以在等式左边得到 $x_1$ ，如果 $x_1 \neq x_0$ ，我们继续将 $x_1$ 代入到等式右边，继续可以在等式左边得到 $x_2$ ，重复进行该操作，我们得到一系列的 $x_1,x_2,...,x_n$ 组成一个数列 $\{x_n\}$ 。这就是格式为 $x_{n+1}=\varphi(x_n)$ 的迭代计算，如果该数列的极限 $\lim_{n \to \infty}x_n$ 存在且等于 $x^*$ ，则称该迭代格式收敛， $x^*$ 就是我们求解的零点，也叫不动点。这个例子是最简单的“不动点迭代法”，其实我们还可以自己构造出各式各样的迭代格式，然而不一定每种迭代格式都是可行的。那我们又该怎么判断一个迭代格式是否可行呢？

首先判断迭代格式的不动点的存在唯一性，然后再判断其收敛性，下面以“不动点迭代法”为例，可以证明以下定理：

(1-a) （存在唯一性）设 $\varphi(x) \in C[a,b],$ 且 $a \leq \varphi(x) \leq b$ 对一切 $x \in [a,b]$ 成立，则 $\varphi(x)$ 在 $[a,b]$ 上一定有不动点。进一步设 $\varphi(x) \in C^1(a,b),$ 且存在常数 $0<L<1$ 使 $|\varphi'(x)| \leq L$ 对一切 $x\in(a,b)$ 成立，则 $\varphi(x)$ 在 $[a,b]$ 上的不动点是唯一的（【2】P25-26证明）。

(1-b) （全局收敛性）设 $\varphi(x) \in C^1[a,b]$ 且满足（1） $a \leq \varphi(x) \leq b$ 对一切 $x \in [a,b]$ 成立；（2）存在常数 $0<L<1，$ 使 $|\varphi'(x)|\leq L,$ 对一切 $x \in (a,b)$ 成立，则对任意的 $x_0 \in [a,b],x_n=\varphi(x_{n-1})$ 产生的序列 $\{x_n\}$ 必收敛到 $\varphi(x)$ 的不动点（【2】P26证明）。

(1-c) （局部收敛性）设 $x^*$ 为 $\varphi(x)$ 的不动点， $\varphi '(x)$ 在 $x^*$ 连续且 $|\varphi '(x^*)|<1$ ，则存在 $x^*$ 的某邻域对任意 $x_0$ 属于该邻域，迭代格式 $x_{n+1}=\varphi(x_n)$ 产生的序列 $\{x_n\}$ 收敛到不动点 $x^*$ 。

当然，一个迭代格式就像极限那样，无穷逼近与不动点，因此我们要对其设置终止准则。同样以“不动点迭代法”为例，可以证明以下定理：

(1-d)（全局收敛性的误差估计）设 $\varphi(x) \in C^1[a,b]$ 且满足（1） $a \leq \varphi(x) \leq b$ 对一切 $x \in [a,b]$ 成立；（2）存在常数 $0<L<1，$ 使 $|\varphi'(x)|\leq L,$ 对一切 $x \in (a,b)$ 成立，则对任意 $x_0 \in [a,b],x_{n+1}=\varphi(x_n)$ 产生的序列满足下面两式： $|x^*-x_n|\leq \frac{1}{1-L}|x_{n+1}-x_n| \leq \frac{L}{1-L}|x_n-x_{n-1}|,\,|x^*-x_n|\leq \frac{L^n}{1-L}|x_1-x_0|$ （【2】P26-27证明）。

“不动点迭代法”是在原式的基础上进行构造的，除此之外，我们还可以采用“牛顿迭代法”，它先在原式的基础上近似后再进行构造的，具体做法为：用 $f(x)$ 的泰勒展开式 $f(x) \approx f(x_0)+f'(x_0)(x-x_0)+\frac{f''(\xi)(x-x_0)^2}{2!}$ 中的线性函数 $f(x_0)+f'(x_0)(x-x_0)$ 近似代替函数 $f(x)$ ，则将求解 $f(x)=0$ 的问题转换为求解 $f(x_0)+f'(x_0)(x-x_0)=0$ 的问题。由此很容易得到牛顿迭代格式为式 $x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}\quad (n=0,1,2,...)$ 。（【3】P32-34证明，牛顿迭代法的存在唯一性、收敛性等）

最后，不同迭代格式之间如何进行优劣比较呢？考察各迭代格式的收敛速度即可。收敛速度的定义为：

设 $\lim_{n \to \infty} x_n=x^*$ （即序列 $\{x_n\}$ 收敛于 $x^*$ ），如果存在 $a>0,p \geq1$ 使得 $\lim_{n \to \infty} \frac{|x_{n+1}-x^*|}{|x_n-x^*|^p}=a$ ，则称数列 $\{x_n\}$ 为 $p$ 阶收敛。特别地：当 $p=1(a<1)$ 时，称为线性收敛；当 $p>1$ 时，称为超线性收敛；当 $p=2$ 时，称为平方收敛。序列的收敛阶数越高，则收敛速度越快。

判断一个迭代格式的收敛阶数可用以下定理：

设 $x^*$ 为方程 $x=\varphi(x)$ 的根。如果迭代函数 $\varphi(x)$ 满足条件：（1） $\varphi(x)$ 在 $x^*$ 邻近是 $p$ 次连续可微的 $(p \geq 2)$ ；（2） $\varphi'(x^*)=\varphi''(x^*)=...=\varphi^{(p-1)}(x^*)=0$ 而 $\varphi^{(p)}(x^*) \neq0$ ，则当初值 $x_0$ 取得充分靠近 $x^*$ 时，迭代格式 $x_{n+1}=\varphi(x_n)$ 是 $p$ 阶收敛的。

至此，迭代思想的核心大致总结为4个方面：“根的存在唯一性”、“收敛性”、“终止准则”和“收敛速度”。

1. 注意：(1) 需要构造迭代格式时，如果没指明需要构造的形式，往往我们可以首先考虑牛顿迭代法；(2) 牛顿迭代法至少二阶收敛的前提是没有重根，但如果遇到有重根的情况则需要更加仔细地代入公式进去算，因为在计算过程中，能够导致为0的元素可能分子分母消去，比如 $f(x)=(x^2-a)^2$ ;(3) 牛顿迭代法擅长处理非线性方程，所以在对于实际问题构造时往往要往非线性方程构造，比如对于同一个问题构造为 $f(x)=x-1/a$ 是无法使用牛顿迭代法的，应该构造为 $f(x)=1/x-a$ 。预备工作：线性方程组类型的判断

首先讨论相容线性方程组的解法。

一. 三角分解法

1.1 QR分解

1.1.1 形式一: 满秩方阵

(a) 对 $A$ 进行列分块，得列分块向量 $(\alpha_1,\alpha_2,...,\alpha_n)$

(b) 对列分块向量施密特正交化，得到正交向量 $(\beta_1,\beta_2,...,\beta_n)$ ，其中

1. 拓展：上述过程中 $A \in \mathbb{R}_n^{n \times n}$ ，同时还可进行 $A=$ $LQ$ 分解，在步骤(a)中行分块即可。类似地， $A \in \mathbb{C}_n^{n \times n}$ 时可进行 $A=UR$ 或 $A=LU$ 分解。区别在于 $Q$ 是正交矩阵，而 $U$ 是酉矩阵。下面开始仅讨论复矩阵的情形，相应的实矩阵情形很好类推。

1.1.2 形式二: 行满秩矩阵

$A=[L\,O]U,A \in \mathbb{C}_m^{m \times n}$ ，其中 $L$ 为 $m$ 阶正线下三角复矩阵， $U$ 为 $n$ 阶酉矩阵， $O$ 为零矩阵。（【1】P92证明）

1.1.3 形式三: 列满秩矩阵

$A=U \begin{bmatrix}R\\O \end{bmatrix},A \in \mathbb{C}_n^{m \times n}$ ，其中 $U$ 为 $m$ 阶酉矩阵， $R$ 为 $n$ 阶正线上三角复矩阵， $O$ 为零矩阵。（【1】P92证明）

1.1.4 形式四: 任意矩阵

1.2 LU分解

1.2.1 形式一: 满秩方阵

$A=LR^*=L^*R=L^*DR^*,A\in \mathbb{C}_n^{n \times n}$ ， $A$ 的各阶顺序主子式行列式不为零， $L,R,D$ 的主对角线上元素不为0。（【1】P89-92证明，4个命题等价）

1. 参数含义： $R$ 为上三角复矩阵， $L$ 为下三角复矩阵， $R^*$ 为单位上三角复矩阵， $L^*$ 为单位下三角复矩阵， $D$ 为对角矩阵。
2. 拓展：同样地， $A\in \mathbb{R}_n^{n \times n}$ 的分解形式类推，相应的所得则为实矩阵。
3. 分解方法：
(1) 待定系数法
(a) 写出含待定系数的分解形式

$\begin{pmatrix} 2&4&4&2 \\3&3&12&6\\2&4&-1&2\\4&2&1&1 \end{pmatrix}= \begin{pmatrix} 1\\l_{21}&1\\l_{31}&l_{32}&1\\l_{41} &l_{42}&l_{43}&1 \end{pmatrix} \begin{pmatrix} u_{11}&u_{12}&u_{13}&u_{14}\\&u_{22}&u_{23}&u_{24}\\&&u_{33}&u_{34}\\&&&u_{44} \end{pmatrix}$
(b) 按照规则“ 先行后列”、“列除行不除”、“旧元素减去其所在行和列前(k-1)个元素的对应乘积然后求和”，对照上式一步步地写出待定系数

$\begin{align} \begin{pmatrix} 2&4&4&2 \\3&3&12&6\\2&4&-1&2\\4&2&1&1 \end{pmatrix}&= \begin{pmatrix} 1\\3/2&1\\1&l_{32}&1\\2 &l_{42}&l_{43}&1 \end{pmatrix} \begin{pmatrix} 2&4&4&2\\&u_{22}&u_{23}&u_{24}\\&&u_{33}&u_{34}\\&&&u_{44} \end{pmatrix}=\begin{pmatrix} 1\\3/2&1\\1&0&1\\2&-2&l_{43}&1 \end{pmatrix} \begin{pmatrix} 2&4&4&2\\&-3&6&3\\&&-5&0\\&&&u_{44} \end{pmatrix}\\&=\begin{pmatrix} 1\\3/2&1\\1&0&1\\2&2&19/5&1 \end{pmatrix} \begin{pmatrix} 2&4&4&2\\&-3&6&3\\&&-5&0\\&&&-9 \end{pmatrix} \end{align}$
(2) 高斯消元法
(a) 搭建框架
$\begin{pmatrix}1\\&1\\&&1\\&&&1\end{pmatrix}\begin{pmatrix}2&4&4&2\\3&3&12&6\\2&4&-1&2\\4&2&1&1\end{pmatrix}=\begin{pmatrix}&\\&\\&\\&&&&&\end{pmatrix}$
(b) 提取高斯消元（即初等行变换）过程中的乘子和结果。(注意行列的对应，第 $i$ 行对第 $j$ 行做初等行变换的乘子，对应到结果矩阵的 $A_{ji}$ 元素）

$\begin{pmatrix}1\\-3/2&1\\-1&&1\\-2&&&1\end{pmatrix}\begin{pmatrix}2&4&4&2\\3&3&12&6\\2&4&-1&2\\4&2&1&1\end{pmatrix}=\begin{pmatrix}2&4&4&2\\0&-3&6&3\\0&0&-5&0\\0&-6&-7&-3\end{pmatrix}\Rightarrow\begin{pmatrix}1\\-3/2&1\\-1&0&1\\-2&-2&&1\end{pmatrix}\begin{pmatrix}2&4&4&2\\3&3&12&6\\2&4&-1&2\\4&2&1&1\end{pmatrix}=\begin{pmatrix}2&4&4&2\\0&-3&6&3\\0&0&-5&0\\0&0&-19&-9\end{pmatrix}\\ \Rightarrow\begin{pmatrix}1\\-3/2&1\\-1&0&1\\-2&-2&-19/5&1\end{pmatrix}\begin{pmatrix}2&4&4&2\\3&3&12&6\\2&4&-1&2\\4&2&1&1\end{pmatrix}=\begin{pmatrix}2&4&4&2\\0&-3&6&3\\0&0&-5&0\\0&0&0&-9\end{pmatrix}\Rightarrow\begin{pmatrix}2&4&4&2\\3&3&12&6\\2&4&-1&2\\4&2&1&1\end{pmatrix}=\begin{pmatrix}1\\3/2&1\\1&0&1\\2&2&19/5&1\end{pmatrix}\begin{pmatrix}2&4&4&2\\0&-3&6&3\\0&0&-5&0\\0&0&0&-9\end{pmatrix}\\$
4. 注意：
(a) 利用高斯消元法对矩阵 $A$ 进行 $LU$ 矩阵分解时，绝对值小的主元可能产生麻烦，因此一般在进行前先判断是否对行进行调整。假设调整后的矩阵是 $B$ ,分解为 $B=L^*U$ ，那么， $A$ 的分解又该表示为什么形式呢？我们用一个单位矩阵 $I$ 来表征这种变换，如果要交换矩阵 $A$ 中的一些行，那我们就对单位矩阵 $I$ 中相应的行进行交换，得到矩阵 $P$ 。因此有 $PA=L^*U\Rightarrow A=P^{-1}L^*U$ 。
(b) 得到 $A=L^*U$ 后，我们可以很容易地基于此，将 $U$ 拆成 $DU^*$ 的形式从而得到 $A=L^*DU^*$ ；再将 $L^*D$ 合并为 $L$ 从而得到 $A=LU^*$ ，过程如下：

$\begin{align} \begin{pmatrix} 2&4&4&2 \\3&3&12&6\\2&4&-1&2\\4&2&1&1 \end{pmatrix} &=\begin{pmatrix} 1\\3/2&1\\1&0&1\\2&2&19/5&1 \end{pmatrix} \begin{pmatrix} 2&4&4&2\\&-3&6&3\\&&-5&0\\&&&-9 \end{pmatrix} =\begin{pmatrix} 1\\3/2&1\\1&0&1\\2&2&19/5&1 \end{pmatrix} \begin{pmatrix} 2&&&\\&-3&&\\&&-5&\\&&&-9 \end{pmatrix}\begin{pmatrix}1&2&2&1\\&1&-2&-1\\&&1&0\\&&&1\end{pmatrix} \\ &=\begin{pmatrix} 2\\3&-3\\2&0&-5\\4&-6&-19&-9\end{pmatrix}\begin{pmatrix}1&2&2&1\\&1&-2&-1\\&&1&0\\&&&1\end{pmatrix} \end{align}$
(c) 对于三对角矩阵，如果它的非零系数在主对角线上和两条次对角线上，那么可以对其进行 $LU$ 分解，上述的两种分解方法都可以使用。

1.2.2 形式二: 行满秩矩阵

1.2.3 形式三: 列满秩矩阵

1.3 Cholesky分解

$A=R^HR,A\in \mathbb{C}_n^{n \times n}$ ， $A$ 是正定Hermite矩阵， $R$ 是正线上三角复矩阵。类似的， $A \in\mathbb{R}_n^{n \times n}$ 则要求 $A$ 为实对称正定矩阵。

1. 拓展：对于 $A$ 同时还是正定Hermite矩阵的情况，它还可以分解为 $A=R^{*H}DR^*$ ，这里的 $R^*$ 为单位上三角复矩阵， $D$ 为对角阵。
2. 分解方法（即待定系数法）：
$\begin{pmatrix}4&-2&2\\-2&2&-4\\2&-4&11\end{pmatrix}=\begin{pmatrix}l_{11}\\l_{12}&l_{22}\\l_{13}&l_{23}&l_{33}\end{pmatrix}\begin{pmatrix}l_{11}&l_{12}&l_{13}\\&l_{22}&l_{23}\\&&l_{33}\end{pmatrix}=\begin{pmatrix}l_{11}^2&l_{11}l_{12}&l_{11}l_{13}\\l_{11}l_{12}&l_{12}^2+l_{22}^2&l_{12}l_{13}+l_{22}l_{23}\\l_{11}l_{13}&l_{12}l_{13}+l_{22}l_{23}&l_{13}^2+l_{23}^2+l_{33}^2\end{pmatrix}$

二. 迭代解法

2.1 迭代思想

首先判断迭代格式的不动点的存在唯一性，然后再判断其收敛性，下面以“不动点迭代法”为例，可以证明以下定理：

当然，一个迭代格式就像极限那样，无穷逼近与不动点，因此我们要对其设置终止准则。同样以“不动点迭代法”为例，可以证明以下定理：

最后，不同迭代格式之间如何进行优劣比较呢？考察各迭代格式的收敛速度即可。收敛速度的定义为：

判断一个迭代格式的收敛阶数可用以下定理：

至此，迭代思想的核心大致总结为4个方面：“根的存在唯一性”、“收敛性”、“终止准则”和“收敛速度”。

1. 注意：(1) 需要构造迭代格式时，如果没指明需要构造的形式，往往我们可以首先考虑牛顿迭代法；(2) 牛顿迭代法至少二阶收敛的前提是没有重根，但如果遇到有重根的情况则需要更加仔细地代入公式进去算，因为在计算过程中，能够导致为0的元素可能分子分母消去，比如 $f(x)=(x^2-a)^2$ ;(3) 牛顿迭代法擅长处理非线性方程，所以在对于实际问题构造时往往要往非线性方程构造，比如对于同一个问题构造为 $f(x)=x-1/a$ 是无法使用牛顿迭代法的，应该构造为 $f(x)=1/x-a$ 。

解决线性方程组的最终目标是找到未知变量的值。这是带有两个未知变量的线性方程组的示例：

等式1：

4x  + 3y = 20
-5x + 9y = 26

为了解决上述线性方程组，我们需要找到x和y变量的值。解决此类系统的方法有多种，例如消除变量，克莱默规则，矩阵解决方案。在本文中，我们将介绍矩阵解决方案。

课程

R语言数据分析挖掘必知必会

从数据获取和清理开始，有目的的进行探索性分析与可视化。让数据从生涩的资料，摇身成为有温度的故事。

立即参加 ➜

在矩阵解中，要求解的线性方程组以矩阵形式表示AX = B。例如，我们可以用矩阵形式表示等式1，如下所示：

A = [[ 4   3]
     [-5   9]]
 
X = [[x]
     [y]]
 
B = [[20]
     [26]]

要查找的值x和y变量方程1，我们需要找到在矩阵中的值X。为此，我们可以采用矩阵逆的点积A和矩阵B，如下所示：

X = inverse(A).B

用numpy求解线性方程组

要求解线性方程组，我们需要执行两个操作：矩阵求逆和矩阵点积。Python的Numpy库支持这两种操作。如果尚未安装Numpy库，则可以使用以下pip命令：

$ pip install numpy

现在让我们看看如何使用Numpy库解决线性方程组。

使用inv（）和dot（）方法

首先，我们将找到A在上一节中定义的矩阵逆。

首先让我们A在Python中创建矩阵。要创建矩阵，array可以使用Numpy模块的方法。矩阵可以视为列表列表，其中每个列表代表一行。

在以下脚本中，我们创建一个名为的列表m_list，其中进一步包含两个列表：[4,3]和[-5,9]。这些列表是矩阵中的两行A。要A使用Numpy 创建矩阵，将m_list传递给array方法，如下所示：

import numpy as np
 
m_list = [[4, 3], [-5, 9]]
A = np.array(m_list)

为了找到矩阵的逆，将矩阵传递给linalg.inv()Numpy模块：

inv_A = np.linalg.inv(A)
 
print(inv_A)

下一步是找出矩阵的逆矩阵之间的点积A和矩阵B。重要的是要提一下，只有在矩阵的维度相等的情况下，才可能在矩阵之间获得矩阵点积，即，左矩阵的列数必须与右矩阵的行数匹配。

最受欢迎的见解

1.用R语言模拟混合制排队随机服务排队系统

2.R语言中使用排队论预测等待时间

3.R语言中实现马尔可夫链蒙特卡罗MCMC模型

4.R语言中的马尔科夫机制转换(Markov regime switching)模型

5.matlab贝叶斯隐马尔可夫hmm模型

6.用R语言模拟混合制排队随机服务排队系统

7.Python基于粒子群优化的投资组合优化

8.R语言马尔可夫转换模型研究交通伤亡人数事故预测

9.用机器学习识别不断变化的股市状况——隐马尔可夫模型的应用

要使用Numpy库查找点积，使用linalg.dot()函数。

B = np.array([20, 26])
X = np.linalg.inv(A).dot(B)
 
print(X)

输出：

[2. 4.]

验证一下，如果在方程式中插入x并4替换未知数，您将看到结果为20。

现在，让我们解决由三个线性方程组成的系统，如下所示：

4x + 3y + 2z = 25
-2x + 2y + 3z = -10
3x -5y + 2z = -4

可以使用Numpy库按以下方式求解以上方程式：

公式2：

  print(X)

在上面的脚本中，linalg.inv()和linalg.dot()方法链接在一起。该变量X包含方程式2的解，并输出如下：

[ 5.  3. -2.]

未知数x，y和的值分别是5、3 z和-2。您可以将这些值代入公式2并验证其正确性。

使用solve（）方法

在前两个示例中，我们使用linalg.inv()和linalg.dot()方法来找到方程组的解。但是，Numpy库包含该linalg.solve()方法，该方法可用于直接找到线性方程组的解：

  print(X2)

输出：

[ 5.  3. -2.]

您可以看到输出与以前相同。

一个真实的例子

让我们看看如何使用线性方程组来解决实际问题。

假设有一个卖水果的人一天就卖出了20个芒果和10个橘子，总价为350元。第二天，他以500元的价格出售了17个芒果和22个橙子。如果这两天的水果价格都保持不变，那么一个芒果和一个橙子的价格是多少？

使用两个线性方程组可以轻松解决此问题。

假设一个芒果x的价格为，一个橙子的价格为y。上面的问题可以这样转换：

20x + 10y = 350
17x + 22y = 500

上面的方程组的解决方案如下所示：

 
X = np.linalg.solve(A,B)
 
print(X)

这是输出：

[10. 15.]

输出显示，一个芒果的价格为10元，一个橙子的价格为15元。

结论

本文介绍了如何使用Python的Numpy库解决线性方程组。您可以使用linalg.inv()和linalg.dot()方法来求解线性方程组，也可以简单地使用solve()方法。solve()方法是首选方法。

可下载资源

关于作者

Kaizong Ye是拓端研究室（TRL）的研究员。在此对他对本文所作的贡献表示诚挚感谢，他在上海财经大学完成了统计学专业的硕士学位，专注人工智能领域。擅长Python.Matlab仿真、视觉处理、神经网络、数据分析。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

非常感谢您阅读本文，如需帮助请联系我们！

用Python的Numpy求解线性方程组

一. 三角分解法

1.1 QR分解

1.2 LU分解

1.3 Cholesky分解

二. 迭代解法

2.1 迭代思想

一. 三角分解法

1.1 QR分解

1.2 LU分解

1.3 Cholesky分解

二. 迭代解法

2.1 迭代思想

用numpy求解线性方程组

结论

相关文章

关注有关新文章的微信公众号