线性方程求解 属于

关于 线性方程求解

还记得你第一次学代数的时候吗，那还是在高中？

通常，我们会拿到一个方程组，然后被要求解出x和y。

4x - 3y = -2
2x + 7y = 16

两个方程，两个未知数：代入消元一下，很快就能看出x = 1, y = 2。

再看一眼上一个例子的左边。它看起来就像矩阵[4 -3; 2 7]乘以一个未知向量[x, y]，得到向量[-2, 16]。

用线性代数记法写就是：A x = b，遵循大写字母表示矩阵、小写字母表示向量的惯例。

A x = 0与零空间

一个重要特例是，所有方程右边都为零。

4x - 3y = 0
2x + 7y = 0

在上面的例子里，唯一的解是x = y = 0，这通常没什么意思。

A x = 0存在非平凡解，当且仅当A的行列式为零。

julia> using LinearAlgebra

julia> A = [4 -3; 2 7]
2×2 Matrix{Int64}:
 4  -3
 2   7

julia> det(A)  # not zero!
34.0

换成下面这组方程试试：

4x - 3y = 0
8x + 6y = 0

第二个方程只是第一个方程的两倍，没有带来任何新信息。只要满足x = 0.75y，任何取值都是解。

这些(x, y)取值落在二维空间中的一条直线上。解有无穷多个，它们构成了这个矩阵的零空间。

这种情况下，矩阵的各行不是线性无关的，矩阵的秩也小于行数或列数。

julia> A = [4 -3; 8 -6]
2×2 Matrix{Int64}:
 4  -3
 8  -6

julia> det(A)
0.0

julia> size(A)  # total number of rows and columns
(2, 2)

julia> rank(A)  # number of linearly-independent rows (or columns)
1

julia> nullspace(A)
2×1 Matrix{Float64}:
 -0.6
 -0.8

julia> nullspace(A) |> norm  # Julia returns unit vectors when possible
1.0

nullspace() 函数返回一个点，它被归一化为单位向量，但这个向量的任意标量倍数也都在零空间中。

求解 A x = b

现在假设你有 1000 个方程、1000 个未知数？如果这听起来很荒唐，别忘了现在无线数字传感器既便宜又多功能（可以测量应变、风速、三轴加速度，什么都行……）。用 1000 个这样的传感器来监测一座现代建筑，比如悬索桥，是完全合理的。还需要有一个控制系统来解释这条数据流，一旦情况变得不妙就能发出警报。

同样是A x = b，其中我们已知A（来自工程设计）和b（传感器的测量值），但需要求出x。

对不熟悉这个的人来说，第一个念头往往是想办法“两边同除以A”，把它移到等式右边。

事实上，A有逆矩阵（大多数方阵都有，但并非全部：行列式必须非零），而且这种算法也确实能算（不过很慢！）。

julia> A = [4 -3; 2 7]
2×2 Matrix{Int64}:
 4  -3
 2   7

julia> b = [-2, 16]
2-element Vector{Int64}:
 -2
 16

# the inverse matrix
julia> inv(A)
2×2 Matrix{Float64}:
  0.205882   0.0882353
 -0.0588235  0.117647

julia> inv(A) * b
2-element Vector{Float64}:
 0.9999999999999999
 2.0

可惜，计算逆矩阵对小矩阵来说就慢，对大矩阵更是慢得离谱。要是你的桥在计算过程中塌了，那可就糟了！

幸好，其他算法要快得多（这里用的是高斯消元法）。

Julia（照搬了 Matlab）直接用反斜杠来求解（严格来说叫“左除”）。

julia> x = A \ b
2-element Vector{Float64}:
 1.0
 2.0

这只是个简单的例子，但有几个细节要注意。

• 各项系数需要按列对齐，缺少的项要补 0。
• 各个方程（也就是矩阵的各行）需要线性无关。如果第 3 行只是第 1 行和第 2 行之和，它就没有带来新信息，问题就是欠定的。

rank() 命令会给出线性无关的行数/列数，所以要让它等于未知变量的个数。

julia> rank(A)
2

长方形矩阵

大概在你还没到青春期的时候，就有人教过你：解N个未知数需要N个方程。

和数学里的大多数事情一样，现实要更微妙一些（不只是上一节提到的线性无关那个细节）。

行数“太少”

只有N-1个方程（也就是矩阵A里的N-1行）时，问题就是欠定的。我们要绝望地放弃吗？

那得看情况！

N-1个方程仍然包含很多信息。从几何上说，要把解确定成 N 维空间里的一个_点_，需要N个方程；但N-1个方程能告诉我们，解一定在某条_直线_上。

这样一来，解有无穷多个，但同时也有无穷多的参数空间是_没有_解的。

你的这条解直线，会不会穿过某个会让正经工程师担心的参数空间区域？或者穿过某个足以成为封闭建筑、禁止公众进入的理由的区域（悬索桥在强侧风里可是会活蹦乱跳的）？也许值得多做点计算来确定这一点！

julia> A3 = [4 -3 1; 2 2 2; 3 -4 -3]
3×3 Matrix{Int64}:
 4  -3   1
 2   2   2
 3  -4  -3

julia> b3 = [1, 12, -14]
3-element Vector{Int64}:
   1
  12
 -14

# fully-determined solution
julia> A3 \ b3
3-element Vector{Float64}:
 1.0
 2.0
 3.0

# remove row 3 from A3 and B3
julia> A2 = A3[1:2, :]
2×3 Matrix{Int64}:
 4  -3  1
 2   2  2

julia> b2 = b3[1:2]
2-element Vector{Int64}:
  1
 12

# an under-determined solution
julia> z1 = A2 \ b2
3-element Vector{Float64}:
 1.594594594594594
 2.445945945945947
 1.9594594594594592

反斜杠求解器给了我们一个解！据非官方说法，这似乎就是范数最小的那个解（从几何上看，就是离原点最近的那个）。

要得到其他解，我们需要A2的nullspace。

把零空间的任意标量倍数加到原解上，都会得到另一个有效解。

# get the nummspace of A2
julia> N = nullspace(A2)
3×1 Matrix{Float64}:
 -0.46499055497527725
 -0.3487429162314578
  0.813733471206735

# add a random multiple of the nullspace to the earlier solution
julia> z = z1 + rand() * N
3×1 Matrix{Float64}:
 1.274331860650258
 2.205748895487695
 2.5199192438620472

# our random z is a valid solution, recovering b2 = [1, 12]
julia> A2 * z
2×1 Matrix{Float64}:
  0.9999999999999947
 12.0

类似地，N-2个方程会把（无穷多个）解限制在一个平面上，原理是一样的。

行数“太多”

相反的情况是方程的个数_多于_未知数，但它们居然还是线性无关的。

在现实工程中，这完全正常，而且被视为一件好事！

传感器精度有限，b向量精度有限，计算中还有噪声。这时解就是对带噪数据的最小二乘拟合。

最简单的做法是用矩阵的pseudoinverse，Julia 把它实现为pinv()函数。

它的用法很像（非奇异）方阵的逆矩阵，给出的是变量的最小二乘估计，而不是精确解。

# create 5-row A and b for 2 variables
julia> A5 = [4 -3; 2 7; -1 2; 1 -1; 3 1]
5×2 Matrix{Int64}:
  4  -3
  2   7
 -1   2
  1  -1
  3   1

julia> rank(A5)
2

# add some random-normal noise to b5
julia> b5n = [-2, 16, 3, -1, 5] + randn(5) * 0.1
5-element Vector{Float64}:
 -1.9337349223581917
 16.024913008059457
  3.0087646177373992
 -0.8675748721176717
  4.914568047308805

# solve to get x ≈ 1, y ≈ 2, using the pseudoinverse
julia> pinv(A5) * b5n
2-element Vector{Float64}:
 1.006117942769543
 1.9962973764508272

特征值与特征向量

上一节讲的是A x = b的解，这等价于我们熟悉的代数。

这一节讲的是A x = λ x的解，其中 λ 是一个标量。

这是一个很有线性代数特色的概念，不过我们可以从几何上理解它：

对于合适的x和 λ 取值，方阵A会把非零向量x的长度按 λ 倍缩放，而不改变它的方向（除非 λ 为负，会把方向反过来）。

这听起来很冷门，但它其实_异常_有用！

术语： λ 的有效取值是A的_特征值_，对应的x取值是A的_特征向量_。很不幸，我们只能接受这些说到一半就换语言的词。

学生通常学的是手算 2×2 矩阵的特征值/特征向量，但用计算机要容易得多（不过还是相当慢，一般情况下，n×n矩阵的复杂度按O(n^3)增长）。

julia> A = rand(-9:9, 2, 2)
2×2 Matrix{Int64}:
 9  5
 3  2

julia> F = eigen(A)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
2-element Vector{Float64}:
  0.27984674554472466
 10.720153254455274
vectors:
2×2 Matrix{Float64}:
 -0.497417  0.945605
  0.867511  0.325317

# eigenvalues
julia> F.values
2-element Vector{Float64}:
  0.27984674554472466
 10.720153254455274

# each column is an eigenvector, normalized to a unit vector
julia> F.vectors
2×2 Matrix{Float64}:
 -0.497417  0.945605
  0.867511  0.325317

一般来说，n×n矩阵有n个特征值，不过这些值不一定是_互不相同_的。可以把它们看作某个n阶多项式的根（这个多项式叫特征多项式），这些根可以重复，而且即使矩阵是实值的，它们也常常是复数。

每个特征向量代表一个方向，它的任意标量倍数也都是有效的特征向量。为了方便后续计算，Julia 返回的是范数为 1 的单位向量。

应用

要讲清楚特征向量为什么重要，本该由一本 500 页的教科书来完成，而不是短短几段话。这个话题渗透在现代应用数学的_太多_领域里了。

从最高层面看，特征向量代表数据集里“最重要”的轴（方向）。特征值则表示每个轴的“相对重要程度”（这要取决于一些关于输入归一化的假设）。

这具体_意味着_什么，取决于应用场景。

主成分分析

在一个典型的数据科学问题里，我们可能有 100 个甚至更多的“特征”，以数据列的形式存放。几乎不可避免地会有噪声、冗余和不需要的相关性。

我们需要做降维，而PCA就是把混乱理出头绪的一种方法：

• 计算数据的协方差矩阵。
• 现在我们得到了一个方阵，接下来计算它的特征值和特征向量。
• 按特征值从大到小排列特征向量（更准确地说，是按绝对值|λ|）。
• 排在最前面的k个特征向量就是你的principal components，其中k远小于原来的特征数量。
• 把整个数据集投影到这k个轴上，然后开始寻找有趣的规律。

以这种形式，PCA 被各种各样的群体使用过：既有想在药物设计中改进分子性质的医学科学家，也有想从民调数据里读懂选民偏好的政治竞选团队。大概还有想向你推销东西的营销团队，不过任何技术都可以用来做好事，也可以用来做坏事。

PCA 不在 Julia 的最小安装里，不过（在 Exercism 之外）MultivariateStats 包里就有你需要的东西。

图像处理

把 PCA 的思路再推进一步：数字图像不过是由像素值构成的矩阵，我们也可以为它计算主成分。

几十年来，这被用在图像压缩里，PCA 可以指导我们在缩小文件大小时哪些内容最值得保留。

越来越多地，PCA 也成为人脸识别等图像分类器的重要组成部分。下次你走过机场时，老大哥可不只是在看着你，他还在用线性代数理解他看到的东西！

机械工程

相对没那么有争议的是，机械部件的主轴就是它惯性张量 I的特征向量（张量不过是换了个名字的矩阵）。

你车上的每个轮子大概都有一块或几块平衡配重，用来把这个张量 I的非对角元素调成零。修车师傅想必是不会去做这些数学计算的（跟飞机设计师和火箭工程师正好相反），但“抖动”不过是描述张量里交叉项的一个日常说法，在这个例子里，这些交叉项会让你的旅途没那么舒服，还会加剧轴承的机械磨损。