Introduction to Linear Algebra学习:第二章 2.4 矩阵运算的规则

0. 这一节到底在讲什么

前面几节我们一直在用矩阵做消元。这一节要回答一个更基础的问题:矩阵到底要怎么"乘"?
矩阵乘法不是把对应位置的数字相乘(那是另一种运算,叫"哈达玛积",这里不讨论)。矩阵乘法有自己的一套规则,而且这套规则是"设计"出来的——之所以这样定义,是因为这样定义之后,矩阵才能用来表示"先做一次变换,再做一次变换"这种连续的操作。
这一节会给出四种理解矩阵乘法 $AB$ 的方式,它们算出来的结果完全一样,只是看问题的角度不同。理解了这四种角度,后面学线性方程组、子空间、特征值的时候会轻松很多。

1. 能不能乘?——维度匹配规则

设矩阵 $A$ 是 $m$ 行 $n$ 列(写作 $m\times n$),矩阵 $B$ 是 $n$ 行 $p$ 列(写作 $n\times p$)。
规则:只有当 $A$ 的列数等于 $B$ 的行数时,才能计算 $AB$。
$$\underset{A}{\underbrace{(m\times n)}}\cdot \underset{B}{\underbrace{(n\times p)}}=\underset{AB}{\underbrace{(m\times p)}}$$
也就是说,中间的两个 $n$ 必须相等,乘完之后这两个 $n$ “消失”,剩下外面的 $m$ 和 $p$。
举个不能乘的例子:如果 $A$ 是 $3\times 2$,$B$ 也是 $3\times 2$,那么 $A$ 有 $2$ 列,$B$ 有 $3$ 行,$2\ne 3$,所以 $AB$ 没有意义。但反过来 $A$ 是 $3\times 2$、$B$ 是 $2\times 1$(列向量)、或者 $B$ 是 $2\times 2$、甚至 $B$ 是 $2\times 20$,都是可以乘的,因为 $A$ 的列数(2)和 $B$ 的行数都对上了。
一个简单的记忆口诀:“内侧数字要相等,外侧数字留下来”。

2. 矩阵加法——比乘法简单得多

在讲乘法之前先说一句加法。两个矩阵能相加,前提是形状完全相同(行数、列数都一样),加法就是对应位置的数字直接相加:
$$\left[\begin{array}{cc}1& 2\\ 3& 4\end{array}\right]+\left[\begin{array}{cc}5& 6\\ 7& 8\end{array}\right]=\left[\begin{array}{cc}6& 8\\ 10& 12\end{array}\right]$$
矩阵乘以一个常数 $c$,就是把每个数字都乘以 $c$:
$$2\left[\begin{array}{cc}1& 2\\ 3& 4\end{array}\right]=\left[\begin{array}{cc}2& 4\\ 6& 8\end{array}\right]$$
矩阵 $-A$ 就是 $c=-1$ 时的情况,所有符号取反。$A+(-A)$ 得到的是零矩阵(所有元素都是 0)。这部分纯粹是常识,记住"形状必须一样"就行。
矩阵里第 $i$ 行第 $j$ 列的那个数字,记作 $a_{ij}$ 或者 $A(i,j)$。比如第一行的元素依次是 $a_{11},a_{12},\dots ,a_{1n}$;左下角的元素是 $a_{m1}$,右下角是 $a_{mn}$。

3. 矩阵乘法的第一种方式:点积法(逐个元素地算)

这是教科书里最常见、最标准的算法。
规则:$AB$ 的第 $i$ 行第 $j$ 列那个数字,等于"$A$ 的第 $i$ 行"和"$B$ 的第 $j$ 列"做点积(对应位置相乘再求和)。
$$(AB{)}_{ij}=(\text{A 的第 }i\text{ 行})\cdot (\text{B 的第 }j\text{ 列})$$
写得更具体一点,如果 $A$ 是 $m\times n$、$B$ 是 $n\times p$:
$$(AB{)}_{ij}=a_{i1}{b}_{1j}+a_{i2}{b}_{2j}+\cdots +a_{in}{b}_{nj}=\sum _{k=1}^n{a}_{ik}{b}_{kj}$$

用一个具体例子感受一下

$$A=\left[\begin{array}{cc}1& 1\\ 0& 1\end{array}\right],B=\left[\begin{array}{cc}2& 1\\ 3& 0\end{array}\right]$$
算 $(AB{)}_{11}$(第1行第1列):取 $A$ 第1行 $(1,1)$,点乘 $B$ 第1列 $(2,3)$:
$$1\times 2+1\times 3=5$$
按照同样的方法把另外三个位置也算出来,最终结果是:
$$AB=\left[\begin{array}{cc}5& 1\\ 3& 0\end{array}\right]$$
这里一共算了 $4$ 个点积,每个点积要做 $2$ 次乘法,所以总共 $8$ 次乘法。

两个极端情况

• 行向量 $\times$ 列向量:$1\times n$ 乘以 $n\times 1$,结果是 $1\times 1$,也就是一个数字。这其实就是普通的"点积"(也叫"内积")。
• 列向量 $\times$ 行向量:$n\times 1$ 乘以 $1\times n$,结果是一整个 $n\times n$ 矩阵!这叫"外积"。很多初学者会惊讶——两个向量也能乘出一整个矩阵。

一般情况要算多少次乘法?

如果 $A$ 是 $m\times n$、$B$ 是 $n\times p$,那么 $AB$ 一共有 $m\times p$ 个位置,每个位置的点积需要 $n$ 次乘法,所以总共需要
$$m\times n\times p\text{ 次乘法}$$
特别地,当 $A,B$ 都是 $n\times n$ 的方阵时,需要 $n^3$ 次乘法。$n=100$ 时大约要算一百万次乘法;$n=2$ 时只要 $2^3=8$ 次,正好对上前面的例子。

拓展小知识:数学家曾经以为 $n=2$ 时 $8$ 次乘法已经是极限了,后来斯特拉森(Strassen)发现用 $7$ 次乘法(代价是多几次加法)也能算出同样的结果。把大矩阵切成 $2\times 2$ 的小块反复用这个技巧,理论上能把 $n^3$ 降到大约 $n^{2.376}$。但这类算法在实践中过于复杂,普通的科学计算还是老老实实用 $n^3$ 的方法。

4. 矩阵乘法的第二种方式:逐列地看(A 乘 B 的每一列)

把 $B$ 拆成一列一列来看。每一次,$A$ 去乘 $B$ 的一整列,得到的就是 $AB$ 的对应那一列。换句话说:

$AB$ 的每一列,都是 $A$ 的各列的某种"组合"(线性组合),组合的系数就来自 $B$ 的那一列。
举个例子,如果
$$A=\left[\begin{array}{cc}1& 0\\ 2& 1\\ 3& 5\end{array}\right],B=\left[\begin{array}{c}1\\ 2\end{array}\right]$$
那么 $AB$(也就是 $A$ 乘以这一列向量)等于:
$$1\times \left[\begin{array}{c}1\\ 2\\ 3\end{array}\right]+2\times \left[\begin{array}{c}0\\ 1\\ 5\end{array}\right]=\left[\begin{array}{c}1\\ 4\\ 13\end{array}\right]$$
这就是"$A$ 的第一列乘以 $1$,加上 $A$ 的第二列乘以 $2$"。这个角度在后面学习"列空间"“线性组合”"线性方程组的解是否存在"时会反复用到,非常重要。

5. 矩阵乘法的第三种方式:逐行地看(A 的每一行乘 B)

这是上面那个角度的镜像。这次把 $A$ 拆成一行一行,每一行(看成一个行向量)去乘整个矩阵 $B$,得到的就是 $AB$ 的对应那一行:

$AB$ 的每一行,都是 $B$ 的各行的某种线性组合,组合系数来自 $A$ 的那一行。
我们之前讲消元的时候,“用消元矩阵 $E$ 左乘 $A$"其实就是这个角度的应用——$E$ 的每一行决定了 $EA$ 的对应行是 $A$ 原来各行的什么组合(比如"第2行减去2倍第1行”)。

6. 矩阵乘法的第四种方式:列乘行,然后求和(很多人不知道的方式)

这是最容易被忽略、但其实非常优雅的一种方式。
做法:把 $A$ 按列拆开,把 $B$ 按行拆开。$A$ 的第 $k$ 列(一个 $m\times 1$ 的列向量)乘以 $B$ 的第 $k$ 行(一个 $1\times p$ 的行向量),这是一个"外积",结果是一整个 $m\times p$ 的矩阵。把 $k=1,2,\dots ,n$ 这些外积矩阵全部加起来,就得到 $AB$:
$$AB=(\text{A 的第1列})(\text{B 的第1行})+(\text{A 的第2列})(\text{B 的第2行})+\cdots +(\text{A 的第n列})(\text{B 的第n行})$$

用 $2\times 2$ 矩阵验证一下

$$A=\left[\begin{array}{cc}a& b\\ c& d\end{array}\right],B=\left[\begin{array}{cc}E& F\\ G& H\end{array}\right]$$
按"列乘行"的方式:
$$AB=\left[\begin{array}{c}a\\ c\end{array}\right]\left[\begin{array}{cc}E& F\end{array}\right]+\left[\begin{array}{c}b\\ d\end{array}\right]\left[\begin{array}{cc}G& H\end{array}\right]=\left[\begin{array}{cc}aE& aF\\ cE& cF\end{array}\right]+\left[\begin{array}{cc}bG& bH\\ dG& dH\end{array}\right]=\left[\begin{array}{cc}aE+bG& aF+bH\\ cE+dG& cF+dH\end{array}\right]$$
跟用点积法算出来的结果完全一样!这说明四种方式只是看问题的角度不同,本质上算的是同一件事,用的乘法次数也完全一样(都是 $mnp$ 次),只是把这些乘法重新归了类。

7. 四种方式的对比表

8. 矩阵乘法不满足交换律!

这是整节课最重要的提醒:一般来说 $AB\ne BA$。
先说一个最直观的原因:形状可能对不上。如果 $A$ 是 $3\times 2$、$B$ 是 $2\times 5$,那么 $AB$ 是 $3\times 5$,而 $BA$ 根本无法计算(B 的列数是5,A 的行数是3,对不上)。即使两边都能算,结果的形状也可能不同。
就算是两个同样大小的方阵,$AB$ 和 $BA$ 也几乎总是不相等。比如:
$$A=\left[\begin{array}{cc}0& 1\\ 0& 0\end{array}\right],B=\left[\begin{array}{cc}1& 0\\ 0& 0\end{array}\right]$$
$$AB=\left[\begin{array}{cc}0& 0\\ 0& 0\end{array}\right],BA=\left[\begin{array}{cc}0& 1\\ 0& 0\end{array}\right]$$
两者完全不同。所以做矩阵乘法的题目时,千万不能随便交换两个矩阵的相乘顺序。
但有一个特殊情况:单位矩阵 $I$ 跟任何方阵都满足交换律,即 $AI=IA=A$。也只有形如 $cI$(数乘单位矩阵)的矩阵能跟所有矩阵交换。这是单位矩阵之所以特殊的原因之一。

9. 矩阵运算的六条法则

矩阵运算遵守一些和普通数字类似的法则,但也有一条经典法则被打破了。

加法的三条法则(都成立)

$$A+B=B+A(\text{交换律})$$
$$c(A+B)=cA+cB(\text{分配律})$$
$$A+(B+C)=(A+B)+C(\text{结合律})$$

乘法的三条法则(交换律被打破!)

$$AB\ne BA\text{(乘法交换律一般不成立)}$$
$$A(B+C)=AB+AC(\text{左分配律})$$
$$(A+B)C=AC+BC(\text{右分配律})$$
$$A(BC)=(AB)C(\text{结合律,括号可以随便挪})$$
最后这条结合律 $A(BC)=(AB)C$ 非常重要,以至于书里专门给它起了名字:

矩阵乘法基本法则:$A$ 乘以 $BC$,等于 $AB$ 乘以 $C$。换句话说,先算 $BC$ 再左乘 $A$,和先算 $AB$ 再右乘 $C$,结果完全一样。
这条法则保证了我们可以放心地写 $ABC$ 而不用加括号——不管你先算哪两个的乘积,最终结果都一样。它的证明思路是:分配律 $A(B+C)=AB+AC$ 可以逐列证明(归结为 $A(b+c)=Ab+Ac$,这其实就是"线性"的核心定义),而结合律的证明虽然写起来比较繁琐,但用处极大,后面证明很多定理都要靠它。

10. 矩阵的幂

当 $A$ 是方阵时,可以定义它的幂:
$$A^p=\underset{p\text{ 个 }A\text{ 相乘}}{\underbrace{A\cdot A\cdot A\cdots A}}$$
矩阵的幂满足和普通数字幂一样的运算规则:
$$A^p\cdot A^q=A^{p+q},(A^p{)}^q=A^{pq}$$
比如 $A^3\cdot A^4=A^7$(7个 $A$ 连乘),而 $(A^3{)}^4=A^{12}$(12个 $A$ 连乘,不是 $3+4=7$,这一点要小心区分)。
当 $p,q$ 可以是零或负数时,上面的规则依然成立,前提是 $A$ 存在"$-1$ 次方",也就是逆矩阵 $A^{-1}$。这时约定 $A^0=I$(类比普通数字里 $2^0=1$)。

注意:矩阵的逆写作 $A^{-1}$,但绝不能写成 $1/A$ 或者 $\frac{1}{A}$(除了在 MATLAB 这种软件的特殊语法里)。不是所有矩阵都有逆,什么时候 $A$ 可逆是下一节(2.5节)要讨论的核心问题。

11. 分块矩阵与分块乘法

一个大矩阵可以被"切"成若干个小矩阵(称为"块",block),只要切法合理,运算时可以把每一块当成一个数字来处理。

直观例子

下面是一个 $4\times 6$ 的矩阵,按照 $2\times 2$ 切成小块,刚好每一块都是单位矩阵 $I$:
$$A=\left[\begin{array}{cccccc}1& 0& 1& 0& 1& 0\\ 0& 1& 0& 1& 0& 1\\ 1& 0& 1& 0& 1& 0\\ 0& 1& 0& 1& 0& 1\end{array}\right]=\left[\begin{array}{ccc}I& I& I\\ I& I& I\end{array}\right]$$
切完之后,这个 $4\times 6$ 的矩阵变成了一个 $2\times 3$ 的"块矩阵",每一块都是 $2\times 2$ 的单位矩阵。显然,这个块矩阵的表达比原始的一堆 0、1 数字清楚多了。

分块加法

如果 $A$ 和 $B$ 形状相同、并且切法(块的大小)也对应一致,那么可以逐块相加,就像普通矩阵加法一样。

分块乘法的规则

规则:只要 $A$ 的"列方向的切法"和 $B$ 的"行方向的切法"匹配,就可以把每一块当成一个数字来做矩阵乘法。
举个最重要的特例:把 $A$ 按"列"切成 $n$ 块(每块是一整列),把 $B$ 按"行"切成 $n$ 块(每块是一整行)。这正是前面讲的第四种乘法方式——分块乘法和"列乘行求和"本质上是一回事:
$$AB=\sum _{k=1}^n(\text{A 的第 }k\text{ 列})(\text{B 的第 }k\text{ 行})$$

应用:分块做消元(Schur 补)

假设一个矩阵被分成四块 $A,B,C,D$:
$$M=\left[\begin{array}{cc}A& B\\ C& D\end{array}\right]$$
我们想通过"消元"把左下角的块 $C$ 变成零块,方法是用 $C{A}^{-1}$(相当于普通消元里的"乘数 $c/a$")去乘第一行的块 $[AB]$,然后从第二行块中减去:
$$\left[\begin{array}{cc}A& B\\ C& D\end{array}\right]⟶\left[\begin{array}{cc}A& B\\ 0& D-C{A}^{-1}B\end{array}\right]$$
这里 $C$ 被消成了零块,而右下角变成了 $D-C{A}^{-1}B$,这个组合有个专门的名字,叫 Schur 补(Schur complement)。它和普通消元中"主元下面那个数字变成 $d-\frac{c}{a}b$"的公式是完全平行的,只不过这里的 $a,b,c,d$ 换成了矩阵块。

12. 应用实例:用矩阵数图中的路径条数

这是一个非常生动的应用,能帮你理解矩阵乘法"到底在干什么"。
考虑一个有 $4$ 个节点的无向图,节点之间用边连接:

用 mermaid 把这张图画出来:

这张图对应一个 $4\times 4$ 的邻接矩阵 $S$:如果节点 $i$ 和节点 $j$ 之间有边,那么 $S_{ij}=1$,否则 $S_{ij}=0$。因为是"无向图"(边没有方向),所以 $S$ 是对称矩阵($S_{ij}=S_{ji}$)。根据上面的连接关系:
$$S=\left[\begin{array}{cccc}0& 1& 1& 0\\ 1& 0& 1& 1\\ 1& 1& 0& 1\\ 0& 1& 1& 0\end{array}\right]$$

关键问题:$S^2$ 是什么意思?

按照点积法的定义:
$$(S^2{)}_{ij}=(\text{S 的第 }i\text{ 行})\cdot (\text{S 的第 }j\text{ 列})=\sum _{k=1}^4{S}_{ik}{S}_{kj}$$
每一项 $S_{ik}{S}_{kj}$ 只有当 $S_{ik}=1$ 且 $S_{kj}=1$ 时才等于 $1$,也就是"$i$ 到 $k$ 有边,并且 $k$ 到 $j$ 也有边"——这恰好就是一条"经过中间节点 $k$、长度为 $2$"的路径!把 $k$ 取遍所有节点再求和,$(S^2{)}_{ij}$ 算出来的就是:

从节点 $i$ 到节点 $j$、恰好走 $2$ 步的路径条数。
比如节点 $2$ 到节点 $3$:可以走 “2→1→3” 或者 “2→4→3”,所以应该有 $2$ 条长度为 2 的路径。
同样的道理可以推广:$S^N$ 的 $(i,j)$ 位置,就是从节点 $i$ 到节点 $j$ 恰好走 $N$ 步的路径条数,因为一条 $N$ 步的路径,可以看成"先走 $N-1$ 步到达某个中间节点 $k$,再从 $k$ 走最后一步到 $j$"——这正好对应了矩阵乘法 $S^N=S^{N-1}\cdot S$ 的点积定义。
这就是为什么"矩阵乘法"天然适合用来数图上的路径——无论是社交网络里"认识的人的朋友",还是公司里员工之间的沟通渠道,背后都是同一套数学。

13. 完整可运行的 C++ 代码

下面这份代码完整实现了本节讲到的核心内容,包括:

1. 用点积法(第一种方式)实现矩阵乘法;
2. 验证矩阵乘法结合律 $A(BC)=(AB)C$;
3. 验证矩阵乘法一般不满足交换律 $AB\ne BA$;
4. 用列乘行求和(第四种方式)实现矩阵乘法,并与点积法结果做对比;
5. 用邻接矩阵的幂 $S^2$ 来数图上长度为 2 的路径条数。
   代码不依赖 <bits/stdc++.h>,只用标准库里明确需要的头文件,可以直接用 g++ 编译运行。

// ===================================================================
// 文件名: matrix_rules_demo.cpp
// 功能  : 演示《线性代数》2.4节"矩阵运算的规则"中的核心知识点
//         1) 点积法做矩阵乘法
//         2) 验证结合律 A(BC) = (AB)C
//         3) 验证乘法一般不满足交换律 AB != BA
//         4) 列乘行求和的方式做矩阵乘法，并与点积法对比
//         5) 用邻接矩阵的平方 S^2 计算图中长度为2的路径条数
// 编译方式: g++ -std=c++17 -O2 matrix_rules_demo.cpp -o matrix_rules_demo
// 运行方式: ./matrix_rules_demo
// ===================================================================
#include <iostream>   // 用于标准输入输出 std::cout
#include <vector>     // 用于 std::vector 容器，存放矩阵的数据
#include <stdexcept>  // 用于抛出异常 std::invalid_argument
#include <iomanip>    // 用于控制输出格式（比如设置打印宽度）
#include <cmath>      // 用于 std::abs，比较浮点数是否相等时要用到
// -------------------------------------------------------------------
// 用一个类来表示矩阵，内部用 二维 vector<vector<double>> 存储数据
// 这样写的好处是：可以重载运算符，让矩阵乘法、加法看起来和数学公式一样自然
// -------------------------------------------------------------------
class Matrix {
public:
    // rows_: 矩阵的行数；cols_: 矩阵的列数
    // data_: 实际存放矩阵元素的二维数组，data_[i][j] 就是第 i 行第 j 列的元素
    int rows_;
    int cols_;
    std::vector<std::vector<double>> data_;
    // 构造函数：创建一个 r 行 c 列、所有元素初始化为 0 的矩阵
    Matrix(int r, int c) : rows_(r), cols_(c) {
        // 用 r 个长度为 c 的向量来初始化二维数组，初始值全部为 0.0
        data_.assign(r, std::vector<double>(c, 0.0));
    }
    // 用一个二维初始化列表直接构造矩阵，方便我们在 main 函数里快速写出矩阵
    // 例如: Matrix A({{1,2},{3,4}}); 就表示矩阵 [[1,2],[3,4]]
    Matrix(std::initializer_list<std::vector<double>> init) {
        data_ = init;                 // 把初始化列表转成 vector<vector<double>>
        rows_ = (int)data_.size();    // 行数就是外层 vector 的大小
        cols_ = rows_ > 0 ? (int)data_[0].size() : 0; // 列数取第一行的长度
    }
    // 重载 () 运算符，方便用 A(i, j) 的方式访问/修改第 i 行第 j 列的元素
    // 注意：这里用的是数学习惯，下标从 0 开始（程序员视角），
    // 而不是教材里从 1 开始的写法，使用时要留意这一点差异
    double& operator()(int i, int j) {
        return data_[i][j];
    }
    const double& operator()(int i, int j) const {
        return data_[i][j];
    }
    // 打印矩阵，方便调试和展示结果
    void print(const std::string& name) const {
        std::cout << name << " (" << rows_ << "x" << cols_ << ") = " << std::endl;
        for (int i = 0; i < rows_; ++i) {
            std::cout << "  [ ";
            for (int j = 0; j < cols_; ++j) {
                // setw(6) 让每个数字占6个字符宽度，方便对齐
                std::cout << std::setw(6) << data_[i][j] << " ";
            }
            std::cout << "]" << std::endl;
        }
    }
};
// -------------------------------------------------------------------
// 函数: multiplyDotProduct
// 功能: 用"点积法"（教材里的第一种方式）计算矩阵乘法 A * B
//       即 (AB)_{ij} = (A 的第 i 行) 点乘 (B 的第 j 列)
// -------------------------------------------------------------------
Matrix multiplyDotProduct(const Matrix& A, const Matrix& B) {
    // 先检查维度是否匹配：A 的列数必须等于 B 的行数
    if (A.cols_ != B.rows_) {
        throw std::invalid_argument("维度不匹配，无法相乘：A的列数必须等于B的行数");
    }
    // 结果矩阵 C 的大小是 (A的行数) x (B的列数)
    Matrix C(A.rows_, B.cols_);
    // 三层循环：i 遍历 C 的每一行，j 遍历 C 的每一列，k 用来做点积求和
    for (int i = 0; i < A.rows_; ++i) {
        for (int j = 0; j < B.cols_; ++j) {
            double sum = 0.0; // 用来累加点积的结果
            for (int k = 0; k < A.cols_; ++k) {
                // 点积公式: sum += a_{ik} * b_{kj}
                sum += A(i, k) * B(k, j);
            }
            C(i, j) = sum; // 把算好的点积结果放进 C 的第 i 行第 j 列
        }
    }
    return C;
}
// -------------------------------------------------------------------
// 函数: multiplyColumnTimesRow
// 功能: 用"列乘行求和"（教材里的第四种方式）计算矩阵乘法 A * B
//       即 AB = sum_{k=1}^{n} (A的第k列) * (B的第k行)
//       这里每一项 (A的第k列)*(B的第k行) 都是一个"外积"矩阵
// -------------------------------------------------------------------
Matrix multiplyColumnTimesRow(const Matrix& A, const Matrix& B) {
    if (A.cols_ != B.rows_) {
        throw std::invalid_argument("维度不匹配，无法相乘：A的列数必须等于B的行数");
    }
    int n = A.cols_;            // n 是 A 的列数，也是 B 的行数
    Matrix C(A.rows_, B.cols_); // 用来累加所有外积矩阵的结果，初始为全0
    // 对每一个 k，把 "A的第k列" 和 "B的第k行" 做外积，再累加进 C
    for (int k = 0; k < n; ++k) {
        for (int i = 0; i < A.rows_; ++i) {
            for (int j = 0; j < B.cols_; ++j) {
                // 外积矩阵的 (i,j) 元素就是 A(i,k) * B(k,j)
                // 把它累加到结果矩阵对应位置上
                C(i, j) += A(i, k) * B(k, j);
            }
        }
    }
    return C;
}
// -------------------------------------------------------------------
// 函数: isApproximatelyEqual
// 功能: 比较两个矩阵是否（在浮点误差范围内）相等，用于验证不同算法结果一致
// -------------------------------------------------------------------
bool isApproximatelyEqual(const Matrix& A, const Matrix& B, double eps = 1e-9) {
    if (A.rows_ != B.rows_ || A.cols_ != B.cols_) return false;
    for (int i = 0; i < A.rows_; ++i) {
        for (int j = 0; j < A.cols_; ++j) {
            // 如果对应位置的差值的绝对值超过 eps，认为不相等
            if (std::abs(A(i, j) - B(i, j)) > eps) return false;
        }
    }
    return true;
}
// -------------------------------------------------------------------
// 函数: matrixPower
// 功能: 计算方阵 A 的 N 次幂 A^N（N 是非负整数），用于后面数路径条数
//       做法很简单：从单位矩阵开始，循环乘 N 次 A
// -------------------------------------------------------------------
Matrix matrixPower(const Matrix& A, int N) {
    int n = A.rows_;
    // 构造一个 n x n 的单位矩阵，作为乘法的初始值（相当于数字里的"1"）
    Matrix result(n, n);
    for (int i = 0; i < n; ++i) {
        result(i, i) = 1.0;
    }
    // 循环做 N 次矩阵乘法： result = result * A
    for (int t = 0; t < N; ++t) {
        result = multiplyDotProduct(result, A);
    }
    return result;
}
// ===================================================================
// 主函数：依次演示前面讲到的知识点
// ===================================================================
int main() {
    std::cout << "========================================" << std::endl;
    std::cout << "演示1: 用点积法计算矩阵乘法 AB" << std::endl;
    std::cout << "========================================" << std::endl;
    // 对应讲解文档里第3节的例子: A = [[1,1],[0,1]], B = [[2,1],[3,0]]
    Matrix A({{1, 1},
              {0, 1}});
    Matrix B({{2, 1},
              {3, 0}});
    Matrix AB = multiplyDotProduct(A, B);
    A.print("A");
    B.print("B");
    AB.print("AB (点积法计算)");
    // 预期结果应该是 [[5,1],[3,0]]，与讲解文档里的手算结果一致
    std::cout << std::endl;
    std::cout << "========================================" << std::endl;
    std::cout << "演示2: 验证结合律 A(BC) = (AB)C" << std::endl;
    std::cout << "========================================" << std::endl;
    // 再随便造一个矩阵 C，用于验证结合律
    Matrix C({{1, 0, 2},
              {0, 1, 1}});
    // 先算 (AB)C：A是2x2，B是2x2，AB是2x2；再乘C(2x3)，得到(AB)C是2x3
    Matrix AB_then_C = multiplyDotProduct(multiplyDotProduct(A, B), C);
    // 再算 A(BC)：B是2x2，C是2x3，BC是2x3；再用A(2x2)左乘，得到A(BC)是2x3
    Matrix BC = multiplyDotProduct(B, C);
    Matrix A_then_BC = multiplyDotProduct(A, BC);
    AB_then_C.print("(AB)C");
    A_then_BC.print("A(BC)");
    if (isApproximatelyEqual(AB_then_C, A_then_BC)) {
        std::cout << ">>> 验证成功: (AB)C 和 A(BC) 完全相等，结合律成立！" << std::endl;
    } else {
        std::cout << ">>> 出现错误: 两者不相等，请检查代码。" << std::endl;
    }
    std::cout << std::endl;
    std::cout << "========================================" << std::endl;
    std::cout << "演示3: 验证矩阵乘法一般不满足交换律 AB != BA" << std::endl;
    std::cout << "========================================" << std::endl;
    Matrix BA = multiplyDotProduct(B, A);
    AB.print("AB");
    BA.print("BA");
    if (!isApproximatelyEqual(AB, BA)) {
        std::cout << ">>> 正如预期: AB 和 BA 不相等，矩阵乘法不满足交换律。" << std::endl;
    } else {
        std::cout << ">>> 这一组矩阵恰好满足 AB = BA（这是巧合，不是一般规律）。" << std::endl;
    }
    std::cout << std::endl;
    std::cout << "========================================" << std::endl;
    std::cout << "演示4: 用'列乘行求和'方式重新计算AB，并与点积法对比" << std::endl;
    std::cout << "========================================" << std::endl;
    Matrix AB_columnRow = multiplyColumnTimesRow(A, B);
    AB_columnRow.print("AB (列乘行求和法计算)");
    if (isApproximatelyEqual(AB, AB_columnRow)) {
        std::cout << ">>> 验证成功: 两种方法算出的AB完全一致！" << std::endl;
    } else {
        std::cout << ">>> 出现错误: 两种方法结果不一致，请检查代码。" << std::endl;
    }
    std::cout << std::endl;
    std::cout << "========================================" << std::endl;
    std::cout << "演示5: 用邻接矩阵的平方 S^2 数图中长度为2的路径条数" << std::endl;
    std::cout << "========================================" << std::endl;
    // 4个节点的无向图，对应讲解文档第12节的例子
    // 节点编号 0,1,2,3 对应教材里的节点 1,2,3,4
    // 边: 1-2, 1-3, 2-3, 2-4, 3-4
    Matrix S({{0, 1, 1, 0},
              {1, 0, 1, 1},
              {1, 1, 0, 1},
              {0, 1, 1, 0}});
    Matrix S2 = matrixPower(S, 2); // 计算 S 的平方
    S.print("邻接矩阵 S");
    S2.print("S^2 (每个位置是对应两点间长度为2的路径条数)");
    // 单独取出节点2到节点3的路径条数（注意下标从0开始，所以是 S2(1,2)）
    std::cout << "节点2到节点3之间，长度为2的路径条数 = "
              << S2(1, 2) << " (应该等于2: 2->1->3 和 2->4->3)" << std::endl;
    return 0;
}

代码运行结果说明

编译运行之后，你会看到：

• 演示1 会算出 $AB=\left[\begin{array}{cc}5& 1\\ 3& 0\end{array}\right]$，和文档第3节手算的结果完全一致；
• 演示2 会确认 $(AB)C$ 与 $A(BC)$ 是同一个矩阵，验证了结合律；
• 演示3 会展示 $AB\ne BA$（在大多数随便选的矩阵下都是这样）；
• 演示4 会用"列乘行求和"重新算一遍 $AB$，并确认和点积法结果一致，说明四种乘法方式只是计算顺序不同，结果相同；
• 演示5 会算出 $S^2$，并验证节点2到节点3之间恰好有2条长度为2的路径，对应"2→1→3"和"2→4→3"这两条路。

14. 用 ASCII 图直观感受"分块矩阵"的结构

下面用纯文本画出"分块矩阵"是怎么把一个大矩阵切成小矩阵的(对应第11节的例子,一个4×6矩阵切成2×3个小块,每块都是2×2的单位矩阵):

原始的 4 行 6 列矩阵 A：
      列1  列2  列3  列4  列5  列6
行1 [  1    0  |  1    0  |  1    0  ]
行2 [  0    1  |  0    1  |  0    1  ]
    ------------------------------------   <- 横向切割线（行方向）
行3 [  1    0  |  1    0  |  1    0  ]
行4 [  0    1  |  0    1  |  0    1  ]
              ^           ^
         竖向切割线   竖向切割线（列方向）
切完之后，看成一个 2 行 3 列的"块矩阵"：
           块列1   块列2   块列3
  块行1  [   I       I       I    ]
  块行2  [   I       I       I    ]
其中每个 I 都代表一个 2x2 的单位矩阵：
        I = [ 1  0 ]
            [ 0  1 ]

15. 小结(本节核心要点回顾)

2.5 逆矩阵（Inverse Matrices）—— 从零理解

一、什么是逆矩阵？先从"做和撤销"说起

想象你穿上了袜子，又穿上了鞋子。如果你要"撤销"这个过程，必须先脱鞋，再脱袜子——顺序是反过来的。
矩阵的逆，做的就是类似的事：如果矩阵 $A$ 对向量 $x$ 做了某种操作变成了 $b$（写成 $Ax=b$），那么逆矩阵 $A^{-1}$ 的作用就是把这个操作"撤销"，让你重新得到 $x$。
定义：对一个方阵 $A$，如果存在一个矩阵 $A^{-1}$，满足
$$A^{-1}A=I\text{且}A{A}^{-1}=I$$
那么就说 $A$ 是"可逆的"（invertible），$A^{-1}$ 叫做 $A$ 的逆矩阵。这里的 $I$ 是单位矩阵——它对任何向量什么都不做，乘了等于没乘：$Ix=x$。
直观地说：
$$A^{-1}Ax=x$$
无论 $A$ 把 $x$ 变成了什么，$A^{-1}$ 都能把它变回原来的 $x$。
类比数字：一个数 $a$ 的逆就是 $\frac{1}{a}$，只要 $a\ne 0$。矩阵的逆是同样的道理，只不过矩阵比数字复杂得多——不是所有方阵都有逆，这正是这一节要研究的核心问题。

二、不是所有矩阵都可逆！怎么判断？

这一节给出了好几种判断方法（书里称为"六个注记"），它们其实都在回答同一个问题：$A$ 有没有逆？

判断方法汇总表

下面逐一讲清楚每一种方法背后的道理。

方法 1：消元主元法（Note 1）

对一个 $n\times n$ 的方阵做消元（高斯消元），如果最终能得到 $n$ 个非零的主元（对角线上的关键数字），那么这个矩阵就一定可逆。换行（行交换）是允许的，不影响这个结论。
这是实际计算中最常用、最直接的判断方式——你甚至不需要算出逆矩阵，只要做消元看主元个数就行。

方法 2：逆矩阵是唯一的（Note 2）

一个矩阵不可能有两个不同的逆。证明很巧妙，靠的是"括号搬家"：
假设 $BA=I$（$B$ 是左逆）并且 $AC=I$（$C$ 是右逆）。考虑乘积 $BAC$，可以从两个方向去算括号：
$$B(AC)=(BA)C$$
左边 $B(AC)=BI=B$，右边 $(BA)C=IC=C$。所以：
$$B=C$$
这说明：左逆和右逆必须是同一个矩阵。这就是为什么我们可以放心地说"$A$ 的逆"，而不必区分"左边的逆"和"右边的逆"。

方法 3：用逆矩阵直接解方程（Note 3）

如果 $A$ 可逆，那么方程 $Ax=b$ 有且只有一个解：
$$x=A^{-1}b$$
这是逆矩阵最直接的用途：两边同时左乘 $A^{-1}$，左边 $A^{-1}Ax=Ix=x$，右边变成 $A^{-1}b$。
但实际计算中很少真的去求 $A^{-1}$——直接用消元法解 $Ax=b$ 通常更快（后面会算一下两者的计算量差异）。

方法 4：零向量检测法（Note 4，非常重要）

关键结论：如果存在一个非零向量 $x$，使得 $Ax=0$，那么 $A$ 绝对不可能有逆。
为什么？因为零向量 $0$ 经过 $A$ 变换后还是 $0$（$A\cdot 0=0$），但现在又有一个非零的 $x$ 也变成了 $0$。这就好比两条不同的路都通向了同一个终点——如果 $A^{-1}$ 存在，它要怎么知道"把 $0$ 变回去"应该变成 $x$ 还是变成 $0$ 呢？它做不到，所以 $A^{-1}$ 不存在。
反过来说：如果 $A$ 真的可逆，那么 $Ax=0$ 就只能有唯一解 $x=A^{-1}0=0$。

方法 5：2×2 矩阵的行列式公式（Note 5）

对于最简单的 2×2 矩阵：
$$A=\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]$$
它可逆的条件是 $ad-bc\ne 0$，并且逆矩阵有现成公式：
$$A^{-1}=\frac{1}{ad-bc}\left[\begin{array}{cc}d& -b\\ -c& a\end{array}\right]$$
这里 $ad-bc$ 这个数，就叫做 $A$ 的行列式（determinant），记作 $\det (A)$。注意公式里要拿这个数做分母，所以一旦它是 $0$，就没法除，矩阵自然没有逆。
举例：矩阵 $A=\left[\begin{array}{cc}1& 2\\ 2& 4\end{array}\right]$ 不可逆，因为 $ad-bc=1\times 4-2\times 2=0$。同时也能验证：取 $x=(2,-1)$，算一下 $Ax$：
$$Ax=\left[\begin{array}{cc}1& 2\\ 2& 4\end{array}\right]\left[\begin{array}{c}2\\ -1\end{array}\right]=\left[\begin{array}{c}1\times 2+2\times (-1)\\ 2\times 2+4\times (-1)\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]$$
正好对应方法 4：存在非零向量 $x=(2,-1)$ 使得 $Ax=0$，所以不可逆。三种判断方法（行列式、零向量、主元个数）在这个例子里完全一致地给出"不可逆"的结论。

方法 6：对角矩阵（Note 6）

对角矩阵（只有主对角线上有非零数字，其余全是 $0$）只要对角线上没有 $0$，就一定可逆，而且逆矩阵也很简单——把每个对角元素取倒数即可。

方法补充：对角占优矩阵一定可逆

如果一个矩阵每一行里，对角线上那个数的绝对值，比该行其余所有数的绝对值之和还大，即：
$$|a_{ii}|>\sum _{j\ne i}|a_{ij}|$$
那么这个矩阵叫"对角占优"（diagonally dominant），它一定可逆，不需要再去算行列式或做消元验证。
直观证明思路：假设有非零向量 $x$，设它绝对值最大的分量是 $x_i$。如果 $Ax=0$，那么 $Ax$ 的第 $i$ 行展开后，$a_{ii}{x}_i$ 这一项的绝对值，会比同一行其他所有项加起来的绝对值还大（因为对角占优），所以它们根本不可能正好抵消为 $0$。这就矛盾了，所以只能 $x=0$，从而 $A$ 可逆。

三、两个矩阵乘积的逆——顺序要反过来！

直观类比

还是穿鞋子的例子：先穿袜子（动作 $B$），再穿鞋子（动作 $A$），合起来的动作是"先 $B$ 后 $A$"，写成矩阵形式就是 $AB$（矩阵作用在最右边的向量上，所以最先生效的写在右边）。
撤销这个过程，必须先脱鞋（撤销 $A$，也就是 $A^{-1}$），再脱袜子（撤销 $B$，也就是 $B^{-1}$）。所以：
$$(AB{)}^{-1}=B^{-1}{A}^{-1}$$
注意顺序反过来了！ 不是 $A^{-1}{B}^{-1}$。

严格证明

要证明 $B^{-1}{A}^{-1}$ 确实是 $AB$ 的逆，只需要把它们乘起来，看是否等于 $I$：
$$(AB)(B^{-1}{A}^{-1})=A(B{B}^{-1})A^{-1}=A\cdot I\cdot A^{-1}=A{A}^{-1}=I$$
这里的关键技巧是"挪动括号"：先把中间的 $B{B}^{-1}$ 消成 $I$，外面的 $A$ 和 $A^{-1}$ 自然就挨在一起，消成 $I$。
同样可以验证反过来乘也成立：$(B^{-1}{A}^{-1})(AB)=I$。所以 $B^{-1}{A}^{-1}$ 确实是 $(AB{)}^{-1}$。
这个规律可以推广到任意多个矩阵相乘：
$$(ABC{)}^{-1}=C^{-1}{B}^{-1}{A}^{-1}$$

一个很微妙的例子：消元矩阵的顺序问题

设 $E$ 是"用第1行消去第2行"的消元矩阵（把第2行减去5倍第1行），$F$ 是"用第2行消去第3行"的消元矩阵（把第3行减去4倍第2行）：
$$E=\left[\begin{array}{ccc}1& 0& 0\\ -5& 1& 0\\ 0& 0& 1\end{array}\right],E^{-1}=\left[\begin{array}{ccc}1& 0& 0\\ 5& 1& 0\\ 0& 0& 1\end{array}\right]$$
$$F=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& -4& 1\end{array}\right],F^{-1}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 4& 1\end{array}\right]$$
现在分别计算 $FE$ 和 $E^{-1}{F}^{-1}$，会发现一个很有意思的现象：

• $FE$（先做 E 再做 F 的消元顺序）：第3行会受到第1行的"二次影响"——因为 $F$ 是在 $E$ 已经改变了第2行之后才作用的，所以 $FE$ 里会出现一个额外的数字 “20”（$5\times 4$）。
• $E^{-1}{F}^{-1}$（撤销的顺序，先撤 F 再撤 E）：因为 $F^{-1}$ 先把第3行加回去，这时候第2行还没被 $E^{-1}$ 改变，所以不会产生这个"20"，第3行不会感受到第1行的二次影响。
  这正是为什么后面学习 $A=LU$ 分解时，要用 $E^{-1}{F}^{-1}\cdots$ 这种"逆序乘起来"的方式——因为这样得到的下三角矩阵 $L$ 里，消元乘数会非常整齐地"原样落位"，没有交叉污染的多余项。

四、怎么算出逆矩阵？——高斯-若尔当消元法（Gauss-Jordan）

基本思路

我们想求解 $A{A}^{-1}=I$。把这件事拆开看：$A^{-1}$ 的每一列乘上 $A$，应该得到单位矩阵 $I$ 对应的那一列。
比如 $A$ 是 $3\times 3$ 矩阵，设 $A^{-1}$ 的三列分别是 $x_1,x_2,x_3$，那么：
$$A{x}_1=e_1=\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right),A{x}_2=e_2=\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right),A{x}_3=e_3=\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)$$
也就是说，求逆矩阵其实是同时解 $n$ 个方程组，而且它们的系数矩阵都是同一个 $A$。

具体做法：增广矩阵 $[A\mid I]$

把 $A$ 和单位矩阵 $I$ 拼在一起，写成增广矩阵 $[A\mid I]$，然后对整行做消元操作（高斯消元 + 继续往上消，即"若尔当"的贡献），目标是把左边的 $A$ 变成 $I$。等左边变成 $I$ 之后，右边自然就变成了 $A^{-1}$：
$$[A\mid I]\stackrel{\text{消元操作}}{\to }[I\mid A^{-1}]$$

完整例子：求矩阵 $K$ 的逆

书中用到的经典例子：
$$K=\left[\begin{array}{ccc}2& -1& 0\\ -1& 2& -1\\ 0& -1& 2\end{array}\right]$$
第一步：写出增广矩阵
$$[K\mid I]=\left[\begin{array}{cccccc}2& -1& 0& 1& 0& 0\\ -1& 2& -1& 0& 1& 0\\ 0& -1& 2& 0& 0& 1\end{array}\right]$$
第二步：正向消元（变成上三角，这是高斯的部分）
把第1行的 $\frac{1}{2}$ 倍加到第2行（消去第2行第1列的 $-1$）：
$$\left[\begin{array}{cccccc}2& -1& 0& 1& 0& 0\\ 0& \frac{3}{2}& -1& \frac{1}{2}& 1& 0\\ 0& -1& 2& 0& 0& 1\end{array}\right]$$
再把新第2行的 $\frac{2}{3}$ 倍加到第3行（消去第3行第2列的 $-1$）：
$$\left[\begin{array}{cccccc}2& -1& 0& 1& 0& 0\\ 0& \frac{3}{2}& -1& \frac{1}{2}& 1& 0\\ 0& 0& \frac{4}{3}& \frac{1}{3}& \frac{2}{3}& 1\end{array}\right]$$
主元是 $2,\frac{3}{2},\frac{4}{3}$，左边三列已经是上三角矩阵 $U$ 了。
第三步：反向消元，往上消零（这是若尔当多走的一步）
用第3行把第2行、第1行里第3列的非零项消掉，再用第2行把第1行里第2列的非零项消掉，最终结果（具体中间过程略，直接给出消元后的结果）：
$$\left[\begin{array}{cccccc}2& 0& 0& \frac{3}{4}& \frac{1}{2}& \frac{1}{4}\\ 0& \frac{3}{2}& 0& \frac{1}{2}& 1& \frac{1}{2}\\ 0& 0& \frac{4}{3}& \frac{1}{4}& \frac{1}{2}& \frac{3}{4}\end{array}\right]$$
第四步：每行除以自己的主元，让左边变成单位矩阵
第1行除以2，第2行除以 $\frac{3}{2}$，第3行除以 $\frac{4}{3}$：
$$[I\mid K^{-1}]=\left[\begin{array}{cccccc}1& 0& 0& \frac{3}{8}& \frac{1}{4}& \frac{1}{8}\\ 0& 1& 0& \frac{1}{4}& \frac{1}{2}& \frac{1}{4}\\ 0& 0& 1& \frac{1}{8}& \frac{1}{4}& \frac{3}{8}\end{array}\right]$$
所以：
$$K^{-1}=\frac{1}{4}\left[\begin{array}{ccc}3& 2& 1\\ 2& 4& 2\\ 1& 2& 3\end{array}\right]$$

关于 $K^{-1}$ 的三个观察

1. 对称性保留：$K$ 是对称矩阵（沿主对角线翻折后不变），它的逆 $K^{-1}$ 也是对称的。
2. 稀疏变稠密：$K$ 是"三对角矩阵"（只有主对角线及其紧邻两条线上有非零数），非常稀疏；但 $K^{-1}$ 算出来却是处处非零的"稠密矩阵"。这正是为什么实际工程中很少真的去算逆矩阵——原本省内存的稀疏结构，求逆之后就完全消失了。
3. 行列式 = 主元之积：三个主元 $2,\frac{3}{2},\frac{4}{3}$ 相乘：
   $$2\times \frac{3}{2}\times \frac{4}{3}=4$$
   这个 $4$ 正是 $K$ 的行列式。能看到 $K^{-1}$ 公式里恰好有"除以 4"，这不是巧合——逆矩阵的计算里必然要除以行列式，这就是为什么行列式一旦是 0，矩阵就没有逆（除以 0 是不允许的）。
   

通用流程图（mermaid）

下面用流程图总结高斯-若尔当法求逆矩阵的完整步骤：

五、求逆矩阵到底"值不值"？——计算成本分析

高斯-若尔当法求 $A^{-1}$，本质上是同时解 $n$ 个方程组 $A{x}_1=e_1,A{x}_2=e_2,\dots ,A{x}_n=e_n$，所以总计算量是：
$$\text{求 }{A}^{-1}\text{ 的乘除运算次数}\approx n^3$$
而如果只是想解一个具体的方程 $Ax=b$（不需要求出 $A^{-1}$，直接消元到底再代入回去），计算量大约是：
$$\text{直接解一个 }Ax=b\text{ 的运算次数}\approx \frac{n^3}{3}$$
结论很直接：如果你只是想解一两个具体的方程，直接消元解 $x$ 比先算出整个 $A^{-1}$ 再做 $A^{-1}b$ 要快得多（大约快 3 倍）。这就是为什么书里反复强调——“在大多数实际问题里，我们根本不需要显式地算出 $A^{-1}$”。

六、三角矩阵的逆——保持"形状不变"的好性质

上三角矩阵的逆仍是上三角矩阵

如果 $A$ 是可逆的上三角矩阵（对角线及其右上方有数字，左下方全是 $0$），那么 $A^{-1}$ 也一定是上三角矩阵。
直观理解：求 $A^{-1}$ 的第 $j$ 列，就是在解 $A{x}_j=e_j$，这里 $e_j$ 这个单位向量除了第 $j$ 个位置是 $1$，其余都是 $0$，特别是它后面（第 $j+1$ 到第 $n$ 个）位置全是 $0$。因为 $A$ 是上三角的，往回代入（back substitution）求解时，这些后面的 $0$ 会一直保持是 $0$，不会被"激活"成非零数。所以算出来的 $x_j$ 这一列，后面 $n-j$ 个位置必然也是 $0$——拼起来正好就是上三角的形状。

下三角矩阵（对角线全是1）的逆，仍是下三角（对角线全是1）

设：
$$L=\left[\begin{array}{ccc}1& 0& 0\\ 3& 1& 0\\ 4& 5& 1\end{array}\right]$$
通过 Gauss-Jordan 消元（具体是用消元矩阵 $E_{21},E_{31},E_{32}$ 依次把 $L$ 化成 $I$），可以算出：
$$L^{-1}=\left[\begin{array}{ccc}1& 0& 0\\ -3& 1& 0\\ 11& -5& 1\end{array}\right]$$
注意这里出现了一个"奇怪"的数字 $11$，它其实来自 $3\times 5-4=11$ 这样的交叉计算——这正是上面提到的"二次影响"问题在逆矩阵里的体现，提醒我们：逆矩阵的元素不是简单地把原矩阵元素取负号就行，中间往往有非线性的组合。

七、完整可运行的 C++ 代码：高斯-若尔当法求逆矩阵

下面这份代码完整实现了高斯-若尔当消元法，输入一个 $n\times n$ 矩阵，输出它的逆矩阵（如果存在）。代码不使用万能头文件，所有用到的标准库都单独列出，并且每一步都配有详细中文注释。

#include <iostream>   // 用于标准输入输出 std::cin, std::cout
#include <vector>     // 用于 std::vector 动态数组，存放矩阵
#include <cmath>      // 用于 std::fabs 取绝对值（判断主元是否接近0）
#include <iomanip>    // 用于 std::setprecision 控制输出的小数位数
// 使用 using 简化类型名，Matrix 就是"二维 vector"，表示一个矩阵
using Matrix = std::vector<std::vector<double>>;
// ---------------------------------------------------------
// 函数：打印矩阵，方便调试和查看中间结果
// ---------------------------------------------------------
void printMatrix(const Matrix& M) {
    // 设置输出精度为4位小数，固定小数格式
    std::cout << std::fixed << std::setprecision(4);
    for (const auto& row : M) {        // 遍历矩阵的每一行
        for (double val : row) {       // 遍历这一行的每一个数
            std::cout << val << "\t";  // 用 tab 分隔，方便对齐
        }
        std::cout << "\n";             // 一行结束后换行
    }
    std::cout << "\n";
}
// ---------------------------------------------------------
// 函数：高斯-若尔当消元法求逆矩阵
// 输入：n x n 的矩阵 A
// 输出：如果 A 可逆，返回 true，并把结果存入 invA；
//      如果 A 不可逆（出现全零主元列），返回 false。
// ---------------------------------------------------------
bool gaussJordanInverse(Matrix A, Matrix& invA) {
    int n = static_cast<int>(A.size());  // 矩阵的阶数（行数=列数）
    // ---------- 第一步：构造增广矩阵 [A | I] ----------
    // aug 是一个 n 行、2n 列的矩阵：左边 n 列放 A，右边 n 列放单位矩阵 I
    Matrix aug(n, std::vector<double>(2 * n, 0.0));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            aug[i][j] = A[i][j];   // 左半部分：复制原矩阵 A
        }
        aug[i][n + i] = 1.0;       // 右半部分：单位矩阵，只有对角线是1
    }
    // ---------- 第二步：对每一列做消元，化简到 [I | A^{-1}] ----------
    for (int col = 0; col < n; ++col) {
        // (a) 选主元：在当前列里，从第 col 行往下找绝对值最大的数，
        //     这样做（"部分主元消元"）可以减少数值误差，也能在
        //     当前行恰好是0时，通过换行避免除以0的问题。
        int pivotRow = col;
        double maxVal = std::fabs(aug[col][col]);
        for (int r = col + 1; r < n; ++r) {
            if (std::fabs(aug[r][col]) > maxVal) {
                maxVal = std::fabs(aug[r][col]);
                pivotRow = r;
            }
        }
        // (b) 如果最大主元都几乎是 0，说明这一列消不出非零主元，
        //     矩阵不可逆，直接返回 false。
        const double EPS = 1e-9;  // 数值容差，比这个还小就当作0
        if (maxVal < EPS) {
            return false;  // A 不可逆
        }
        // (c) 如果最大主元不在当前行，就交换两行（对应"行交换"操作）
        if (pivotRow != col) {
            std::swap(aug[pivotRow], aug[col]);
        }
        // (d) 把当前行除以主元，让主元变成1
        //     这一步对应书里"每行除以自己的主元"
        double pivotVal = aug[col][col];
        for (int j = 0; j < 2 * n; ++j) {
            aug[col][j] /= pivotVal;
        }
        // (e) 用当前这一行，把其他所有行在这一列的数消成0
        //     这一步同时完成了"正向消元"（消下面）和
        //     "反向消元"（消上面），这正是 Gauss-Jordan 比
        //     普通 Gauss 多做的部分。
        for (int r = 0; r < n; ++r) {
            if (r == col) continue;          // 自己不用消自己
            double factor = aug[r][col];     // 要消掉的那个数
            if (std::fabs(factor) < EPS) continue;  // 已经是0就跳过
            for (int j = 0; j < 2 * n; ++j) {
                aug[r][j] -= factor * aug[col][j];
            }
        }
    }
    // ---------- 第三步：把增广矩阵右半部分取出，就是 A^{-1} ----------
    invA.assign(n, std::vector<double>(n, 0.0));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            invA[i][j] = aug[i][n + j];
        }
    }
    return true;  // 成功求出逆矩阵
}
// ---------------------------------------------------------
// 函数：矩阵乘法 C = A * B，用来验证 A * A^{-1} 是否等于单位矩阵
// ---------------------------------------------------------
Matrix matMul(const Matrix& A, const Matrix& B) {
    int n = static_cast<int>(A.size());     // A 的行数
    int m = static_cast<int>(B[0].size());  // B 的列数
    int k = static_cast<int>(B.size());     // A 的列数 = B 的行数
    Matrix C(n, std::vector<double>(m, 0.0));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < m; ++j) {
            double sum = 0.0;
            for (int t = 0; t < k; ++t) {
                sum += A[i][t] * B[t][j];   // 点积法：第i行 点乘 第j列
            }
            C[i][j] = sum;
        }
    }
    return C;
}
// ---------------------------------------------------------
// 主函数：用书中的矩阵 K 作为示例，验证求逆是否正确
// ---------------------------------------------------------
int main() {
    // 书中的经典例子 K：三对角矩阵，对角线是2，旁边是-1
    Matrix K = {
        { 2.0, -1.0,  0.0},
        {-1.0,  2.0, -1.0},
        { 0.0, -1.0,  2.0}
    };
    std::cout << "原矩阵 K：\n";
    printMatrix(K);
    Matrix invK;  // 用来存放求出的逆矩阵
    // 调用高斯-若尔当消元函数求逆
    bool ok = gaussJordanInverse(K, invK);
    if (!ok) {
        // 如果矩阵不可逆，打印提示并结束程序
        std::cout << "矩阵 K 不可逆（出现了全零主元列）。\n";
        return 0;
    }
    std::cout << "求得的逆矩阵 K^{-1}：\n";
    printMatrix(invK);
    // 验证：K 乘以 K^{-1} 应该等于单位矩阵 I
    Matrix shouldBeI = matMul(K, invK);
    std::cout << "验证 K * K^{-1}（应该等于单位矩阵）：\n";
    printMatrix(shouldBeI);
    // 额外验证：构造一个不可逆的矩阵（行列式为0），测试程序能否正确识别
    Matrix singularA = {
        {1.0, 2.0},
        {2.0, 4.0}   // 第二行正好是第一行的2倍，ad-bc = 1*4-2*2 = 0
    };
    Matrix invSingular;
    bool okSingular = gaussJordanInverse(singularA, invSingular);
    std::cout << "测试不可逆矩阵 [[1,2],[2,4]]：\n";
    if (!okSingular) {
        std::cout << "程序正确识别：该矩阵不可逆。\n";
    } else {
        std::cout << "程序出错：本应不可逆，却返回了结果。\n";
    }
    return 0;
}

代码运行结果说明

这份代码会输出三部分内容：

1. 原矩阵 $K$。
2. 用高斯-若尔当法求出的 $K^{-1}$（应当与本文第四节手算的结果 $\frac{1}{4}\left[\begin{array}{ccc}3& 2& 1\\ 2& 4& 2\\ 1& 2& 3\end{array}\right]$ 一致）。
3. 验证 $K\times K^{-1}$ 是否等于单位矩阵 $I$（对角线是1，其余是0，允许有极小的浮点误差）。
4. 用矩阵 $\left[\begin{array}{cc}1& 2\\ 2& 4\end{array}\right]$（行列式为0）测试程序是否能正确识别"不可逆"的情况。
   代码中的部分主元选取（选当前列中绝对值最大的数作为主元）是工程实现中的标准做法，可以避免主元恰好是0或者数值非常小时带来的计算误差问题，这是对书中"理想化消元"的一个实用补充。

八、本节核心知识点总结