线性代数基础 属于

关于 线性代数基础

什么是“线性代数”？

技术性的定义有很多，维基百科上有，教科书里有，像 3Blue1Brown 这样优秀的网站上也讲得很清楚。不过对我们来说，可以讲得更随意一些。

Julia（以及 Python）的维护者很喜欢这类东西。 Exercism 的管理层则_不太感冒_。

本概念只会涉及这门庞大领域中比较简单的部分，不过要有心理准备：这不可避免地是一个相当偏数学的概念。

什么是向量？

这取决于你问的是谁，也取决于你想怎样把这么抽象的东西可视化。

• 在 Julia 中，Vector{T}只是Array{T, 1}的别名，而eltype的类型参数T可以是任何类型。
• 在（大部分）物理学中，向量是 N 维空间里一支有length和direction的箭头，但没有固定的位置。
• 在线性代数中，物理学家的那支箭头被固定在空间中，尾端位于原点。这样一来，箭杆就是多余的，于是我们可以把这个向量表示为 N 维空间中的一个point，其中N个元素表示沿这N条轴到原点的距离。

就目前而言，先忽略由字符串或字符组成的向量。 本概念中我们处理的是数值类型：Int、Float 或 Complex。

什么是矩阵？

同样，你会找到各种各样的答案。

• 一个由数字组成的矩形二维数组_（没有参差不齐的边缘）_。方阵是一种常见的特殊情况。
• 由并排堆叠的列向量构成的linear combination。
• 一个可以作用于向量的transformation，就像function可以作用于其他输入类型一样。

有些方阵实在太常见，以至于有专门的名称。

对角矩阵： 所有非零元素都在main diagonal（从左上到右下）上。

julia> [1 0 0; 0 2 0; 0 0 3]
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

# a convenient shortcut
julia> diagm(1:3)
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

单位矩阵： 对角线上全是 1 的对角矩阵（原因之后会更清楚）。 通常简写为I。

julia> Matrix{Float64}(I, 3, 3)
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

上三角矩阵： 对角线及其上方的值非零，下方为零。 下三角矩阵就留给你自己猜了。

矩阵的转置

如果你确实想交换行和列，函数permutedims()可以做到，并会返回一个新矩阵。

这会涉及复制，而在处理大型矩阵时，复制既慢又费内存。

从线性代数的角度看，transpose()函数更实用，因为它会快速地在原矩阵外面创建一个惰性包装。

更实用的是adjoint()函数，它还会把复数的虚部符号取反（如果你想知道这为什么如此重要，一个简短的答案是量子力学）。 这个操作十分常见，所以只要在变量名后面加一个撇号'，就能得到伴随。

julia> m
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

# new, full copy
julia> permutedims(m)
3×2 Matrix{Int64}:
 1  4
 2  5
 3  6

# lazy version
julia> transpose(m)
3×2 transpose(::Matrix{Int64}) with eltype Int64:
 1  4
 2  5
 3  6

julia> mc = [1+2im 2+3im; 3+2im 1+2im]
2×2 Matrix{Complex{Int64}}:
 1+2im  2+3im
 3+2im  1+2im

# lazy conjugate transpose
julia> adjoint(mc)
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1-2im  3-2im
 2-3im  1-2im

# syntactic sugar for adjoint
julia> mc'
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1-2im  3-2im
 2-3im  1-2im

乘法

本文档余下的内容会用到LinearAlgebra模块中的功能。 加上一行using LinearAlgebra就能把它引入命名空间，但我们在示例里不会反复写上这一行（太影响观感）。

向量：逐元素乘积

这在向量运算概念中已经讲过。 运算符是.*，它成对地作用于输入向量，输出的尺寸和类型与输入相同。

julia> [1, 2] .* [3, 4]
2-element Vector{Int64}:
 3
 8

向量：点积

这个极其常见的运算在教科书里写作u ⋅ v，叫做“点积”。

点积等价于逐元素乘积的_和_。

两个向量可以用通常的*运算符相乘，但前提是左边的向量是adjoint：方便起见写作u' * v。 这一点在后面的矩阵乘法一节里会更清楚。

在 Julia 中，居中抬高的点（中间的点）作为dot()函数的语法糖提供（输入\cdot再按 Tab 键）。 不需要指定伴随，因为这个细节会自动处理。

julia> using LinearAlgebra

# with element-wise syntax
julia> sum([1, 2] .* [3, 4])
11

# with u' * v syntax
julia> [1, 2]' * [3, 4]
11

# with dot()
julia> dot([1, 2], [3, 4])
11

# with \cdot syntax
julia> [1, 2] ⋅ [3, 4]
11

对于复向量，左边的向量需要取共轭（虚部符号取反）。 dot函数会自动做这件事，u'语法则显式地做，但如果u和v是复数，sum(u .* v)会失败。

向量：叉积

对于以前上过电磁学课程的人来说，这是一点怀旧！ 对于熟悉计算力矩或角动量向量的工程师来说也是如此。

点积把两个向量转换成一个scalar，而叉积则把两个三维向量转换成第三个三维向量。

在向量空间的几何表示中，新向量_垂直_于包含这两个输入向量的平面。 如果两个输入平行（只差一个符号），它们就无法确定一个平面，因此输出会是[0, 0, 0]。

顺序很重要：u × v == -(v × u)。 这就是著名的右手定则，它让我们许多人盯着自己的大拇指和另外两根手指，在空间里扭来扭去（往往还一脸困惑）。

# with \times syntax
julia> [1, 2, 3] × [3, 4, 5]
3-element Vector{Int64}:
 -2
  4
 -2
# with cross()
julia> cross([1, 2, 3], [3, 4, 5])
3-element Vector{Int64}:
 -2
  4
 -2

注意，这个运算只限于长度为 3 的向量（相当于带有正交x, y, z轴的欧几里得空间）。 数学家们乐于把这一运算定义到其他维度，结果就是一堆让人看不懂的术语大杂烩（高阶张量、楔积、外积、多重向量……）。 到这一步，我们大多数人不是逃走，就是赶紧换个话题！

范数

一个向量有多“大”？

norm就是试图捕捉这一点，它把向量约简为一个合适的标量。

人们定义了整整一族范数，但迄今为止最常见的是 2-范数，也就是√(v ⋅ v)。

这种均方根运算就是（N 维空间中）到原点的_勾股距离_。 如果我们把向量想象成一支尾端在原点的箭头，那么 2-范数就是这支箭头的_长度_。

把p作为第二个参数传入，就能计算任意 p-范数。 1-范数有时很有用：它不过是各元素绝对值之和（所以算起来非常快、非常简单）。

# defaults to the 2-norm
julia> norm([1, 2, 3])
3.7416573867739413
# the 1-norm
julia> norm([1, -2, 3], 1)
6.0

矩阵乘法

有时感觉全世界的应用数学家在过去 80 年里，大部分时间都在把各种计算变成一连串的矩阵乘法。

这是一种计算机非常擅长的运算：

• 它高度重复。
• 它可以被高效并行化。
• 人们开发了大量专用硬件来让它更快，从 1970 年代价值 8000 万美元的 Cray-1，到你笔记本电脑里很可能集成的 GPU。

细节相当简单，尽管乍一看不太直观。

考虑矩阵A乘以向量v得到输出w（按照线性代数的惯例，我们用大写字母表示矩阵，小写字母表示向量）。

A的第一行与v做点积，得到w的第一个元素，第二行得到第二个元素，依此类推。

julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4
julia> v = [5, 6]
2-element Vector{Int64}:
 5
 6
julia> A * v
2-element Vector{Int64}:
 17  # equals [1, 2] ⋅ [5, 6]
 39  # equals [3, 4] ⋅ [5, 6]

矩阵 * 矩阵的乘法把这一做法扩展到右边矩阵的各列。

对于C = A * B，我们可以认为A分别乘以B的每一列，得到C中对应的一列：也就是一连串的矩阵-向量乘法。

等价地说，我们也可以说A的第一行分别与B的每一列做点积，得到C的第一行，第二行得到第二行，依此类推。

这样的思维模型有好几种，感兴趣的同学可以看一整堂麻省理工学院的讲座来了解它们。

julia> A
2×2 Matrix{Int64}:
 1  2
 3  4

julia> B = [5 7; 6 8]
2×2 Matrix{Int64}:
 5  7
 6  8

# left column is the same as A*v previously
julia> A * B
2×2 Matrix{Int64}:
 17  23
 39  53

julia> B * A
2×2 Matrix{Int64}:
 26  38
 30  44

如上面的例子所示，矩阵乘法_不满足交换律_：A*B和B*A之间没有简单的关系。

对于刚接触矩阵乘法的人来说，光靠读文字大概很难把它想象清楚。 YouTube 上有很多用图形演示它的视频，所以搜索“matrix multiplication”，然后挑一个你喜欢的风格、详细程度和语言的视频。

维度

两个向量的点积要求它们的长度相等。

推而广之，对于矩阵乘法，左边矩阵的_列数_必须等于右边矩阵的_行数_。

把尺寸表示为(nrows, ncols)这样的元组，也就是size(A)的输出，就有(a, b) * (b, c) -> (a, c)。 “内部”维度，这里指b和b，是可以做点积的。 “外部”维度，这里指a和c，决定输出的维度。

一个矩形矩阵的例子：

julia> D = reshape(1:6, 2, 3)
2×3 reshape(::UnitRange{Int64}, 2, 3) with eltype Int64:
 1  3  5
 2  4  6

julia> E = reshape(1:12, 3, 4)
3×4 reshape(::UnitRange{Int64}, 3, 4) with eltype Int64:
 1  4  7  10
 2  5  8  11
 3  6  9  12

julia> D * E
2×4 Matrix{Int64}:
 22  49   76  103
 28  64  100  136

# (3, 4) * (2, 3) not possible
julia> E * D
ERROR: DimensionMismatch: matrix A has axes (Base.OneTo(3),Base.OneTo(4)), matrix B has axes (Base.OneTo(2),Base.OneTo(3))

以矩阵的方式成对相乘向量，有两种可能。

按照惯例，我们会用u' * v来等价于点积。 它的维度是(1, 3) * (3, 1)，而 Julia 会把(1, 1)的输出简化为一个标量（这一点不像 R 等语言）。

另一种做法是用u * v'，它的维度是(3, 1) * (1, 3) -> (3, 3)，可以乘遍所有可能的元素对，把两个向量扩展成一个矩阵。

julia> u = [1, 2]
2-element Vector{Int64}:
 1
 2
julia> v = [3, 4]
2-element Vector{Int64}:
 3
 4
julia> u' * v
11
julia> u * v'
2×2 Matrix{Int64}:
 3  4
 6  8

u' * v有时被称为_内积_。

与之对应的，（不太常见的）u * v'则是_外积_。 它与张量积有关，不过那远远超出了我们的范围。

旋转

一种特别常见的矩阵乘法涉及旋转矩阵。

在二维中，有一个相对简单的矩阵，可以把向量逆时针旋转θ弧度。

julia> rot2d(θ, vec) = [cos(θ) -sin(θ); sin(θ) cos(θ)] * vec
rot2d (generic function with 1 method)

# unit vector in the x direction
julia> i_hat = [1, 0]
2-element Vector{Int64}:
 1
 0

# rotate 45 degrees
julia> rot2d(π/4, i_hat)
2-element Vector{Float64}:
 0.7071067811865476
 0.7071067811865475

# rotate 90 degrees -> unit vector in the y direction, j_hat
julia> rot2d(π/2, i_hat)
2-element Vector{Float64}:
 6.123233995736766e-17  # zero, within numerical error
 1.0

三维中的旋转需要一个更复杂的矩阵，涉及两个角度。 不在 Markdown 文档里嵌入大量 LaTeX 的话，这很难展示，所以公式请查看维基百科。

对于三维图形，比如 OpenGL 及其后续技术，惯例是使用齐次坐标。

每个顶点用一个四维向量（或者等价地，一个矩阵列）表示：对于位于(x, y, z)的点，就是[x, y, z, 1.0]。

这样就能构造更复杂的变换矩阵。 旋转仍然位于A[1:3, 1:3]，平移是A[1:3, 4]处的[Δx, Δy, Δz]，缩放位于对角线上，此外还有斜切、透视等其他可能。

性能

在世界各地的计算机系里，有许多博士论文用一章或几章的篇幅讨论矩阵乘法算法上小小的渐进改进。 这_真的_很重要，各种机构也愿意为了自身利益资助相关研究。

性能是一个很大的话题，我们没法深入展开，但为了说明问题，我们可以试着把各种大小的随机矩阵相乘。

下面的代码只是一个快速粗糙的估算（想要更好的做法，请使用 BenchmarkTools.jl）。

所用的机器是一台 430 美元的小型 PC：Ryzen 9 处理器、32GB 内存、Linux Mint 22.1、Julia 1.11.6。

julia> mmul(A) = A * A
mmul (generic function with 1 method)
julia> A = rand(Float64, 1_000, 1_000);
julia> @time mmul(A);
  0.011055 seconds (3 allocations: 7.629 MiB)
julia> A = rand(Float64, 10_000, 10_000);
julia> @time mmul(A);
  7.236284 seconds (1.61 k allocations: 763.027 MiB, 17.92% gc time, 0.15% compilation time)

大致说来，一对百万级元素的矩阵花了几毫秒，一亿级元素的矩阵花了几秒。 再往上加几个数量级，就需要更好的硬件了……