n 维数组上的线性代数#
前提条件#
在阅读本教程之前，你应该了解一些 Python 知识。如果你想复习一下，可以看看Python 教程。
如果你想能够运行本教程中的示例，你还应该在你的电脑上安装matplotlib 和SciPy。
学习者画像#
本教程适合那些对 NumPy 中的线性代数和数组有基本了解，并希望了解如何表示和操作 n 维 ( \(n>=2\) ) 数组的人。特别是，如果你不知道如何在不使用 for 循环的情况下将常用函数应用于 n 维数组，或者如果你想了解 n 维数组的轴和形状属性，本教程可能会对你有所帮助。
学习目标#
在本教程之后，你应该能够
理解 NumPy 中一维、二维和 n 维数组之间的区别；
理解如何在不使用 for 循环的情况下将一些线性代数运算应用于 n 维数组；
理解 n 维数组的轴和形状属性。
内容#
在本教程中，我们将使用线性代数中的矩阵分解，奇异值分解，来生成图像的压缩近似。我们将使用 scipy.datasets 模块中的 face 图像
# TODO: Rm try-except with scipy 1.10 is the minimum supported version
try:
from scipy.datasets import face
except ImportError: # Data was in scipy.misc prior to scipy v1.10
from scipy.misc import face
img = face()
Downloading file 'face.dat' from 'https://raw.githubusercontent.com/scipy/dataset-face/main/face.dat' to '/home/circleci/.cache/scipy-data'.
注意：如果你愿意，你可以在完成本教程时使用你自己的图像。为了将你的图像转换为可以操作的 NumPy 数组，你可以使用 matplotlib.pyplot 子模块中的 imread 函数。或者，你可以使用 imageio 库中的imageio.imread 函数。请注意，如果你使用你自己的图像，你可能需要调整下面的步骤。有关图像转换为 NumPy 数组时的处理方式的更多信息，请参阅 scikit-image 文档中的NumPy 图像速成课程。
现在，img 是一个 NumPy 数组，当我们使用 type 函数时可以看到这一点
type(img)
numpy.ndarray
我们可以使用matplotlib.pyplot.imshow 函数和特殊的 iPython 命令 %matplotlib inline 来内联显示图像。
import matplotlib.pyplot as plt
%matplotlib inline
plt.imshow(img)
plt.show()
形状、轴和数组属性#
请注意，在线性代数中，向量的维数指的是数组中条目的数量。在 NumPy 中，它定义的是轴的数量。例如，一维数组是一个向量，例如 [1, 2, 3]，二维数组是一个矩阵，依此类推。
首先，让我们检查数组中数据的形状。由于此图像是二维的（图像中的像素形成一个矩形），我们可能期望二维数组来表示它（一个矩阵）。但是，使用此 NumPy 数组的 shape 属性会给我们一个不同的结果
img.shape
(768, 1024, 3)
输出是一个包含三个元素的元组，这意味着这是一个三维数组。事实上，由于这是一个彩色图像，我们使用了 imread 函数来读取它，数据被组织在三个二维数组中，代表颜色通道（在本例中是红、绿和蓝 - RGB）。你可以通过查看上面的形状来看到这一点：它表明我们有一个包含 3 个矩阵的数组，每个矩阵的形状为 768x1024。
此外，使用此数组的 ndim 属性，我们可以看到
img.ndim
3
NumPy 将每个维度称为一个 *轴*。由于 imread 的工作方式，*第三个轴中的第一个索引*是图像的红色像素数据。我们可以使用以下语法访问它
img[:, :, 0]
array([[121, 138, 153, ..., 119, 131, 139],
[ 89, 110, 130, ..., 118, 134, 146],
[ 73, 94, 115, ..., 117, 133, 144],
...,
[ 87, 94, 107, ..., 120, 119, 119],
[ 85, 95, 112, ..., 121, 120, 120],
[ 85, 97, 111, ..., 120, 119, 118]], dtype=uint8)
从上面的输出中，我们可以看到 img[:, :, 0] 中的每个值都是 0 到 255 之间的整数值，代表每个对应图像像素中的红色水平（请记住，如果你使用你自己的图像而不是 scipy.datasets.face，这可能会有所不同）。
正如预期的那样，这是一个 768x1024 的矩阵
img[:, :, 0].shape
(768, 1024)
由于我们将对这些数据执行线性代数运算，因此在矩阵的每个条目中使用 0 到 1 之间的实数来表示 RGB 值可能会更有意义。我们可以通过设置来做到这一点
img_array = img / 255
此操作（将数组除以标量）由于 NumPy 的广播规则而有效。（请注意，在实际应用中，最好使用例如 scikit-image 中的img_as_float 实用程序函数）。
你可以通过进行一些测试来检查上述内容是否有效；例如，查询此数组的最大值和最小值
img_array.max(), img_array.min()
(np.float64(1.0), np.float64(0.0))
或检查数组中数据的类型
img_array.dtype
dtype('float64')
请注意，我们可以使用切片语法将每个颜色通道分配给单独的矩阵
red_array = img_array[:, :, 0]
green_array = img_array[:, :, 1]
blue_array = img_array[:, :, 2]
轴上的运算#
可以使用线性代数方法来逼近现有数据集。在这里，我们将使用SVD（奇异值分解）来尝试重建一个使用比原始图像更少的奇异值信息的图像，同时仍然保留其一些特征。
注意：我们将使用 NumPy 的线性代数模块 numpy.linalg 来执行本教程中的操作。此模块中的大多数线性代数函数也可以在scipy.linalg 中找到，并鼓励用户在实际应用中使用scipy 模块。但是，scipy.linalg 模块中的一些函数，例如 SVD 函数，仅支持二维数组。有关此方面的更多信息，请查看scipy.linalg 页面。
要继续，请从 NumPy 导入线性代数子模块
from numpy import linalg
为了从给定矩阵中提取信息，我们可以使用 SVD 来获得 3 个数组，这些数组可以相乘以获得原始矩阵。根据线性代数理论，给定一个矩阵 \(A\)，可以计算以下乘积
其中 \(U\) 和 \(V^T\) 是方阵，而 \(\Sigma\) 与 \(A\) 大小相同。\(\Sigma\) 是一个对角矩阵，包含 \(A\) 的奇异值，从大到小排列。这些值总是非负的，可以用作矩阵 \(A\) 表示的一些特征的“重要性”指标。
让我们首先看看这在一个矩阵中是如何工作的。请注意，根据色彩测量法，如果我们应用以下公式，可以获得彩色图像相当合理的灰度版本
其中 \(Y\) 是表示灰度图像的数组，\(R\)、\(G\) 和 \(B\) 是我们最初拥有的红、绿和蓝通道数组。注意我们可以为此使用 @ 运算符（NumPy 数组的矩阵乘法运算符，参见numpy.matmul）
img_gray = img_array @ [0.2126, 0.7152, 0.0722]
现在，img_gray 的形状为
img_gray.shape
(768, 1024)
要查看这在我们的图像中是否有意义，我们应该使用 matplotlib 中的与我们希望在图像中看到的颜色相对应的颜色图（否则，matplotlib 将默认为与真实数据不对应的颜色图）。
在我们的例子中，我们正在逼近图像的灰度部分，因此我们将使用颜色图 gray
plt.imshow(img_gray, cmap="gray")
plt.show()
现在，将linalg.svd 函数应用于此矩阵，我们得到以下分解
U, s, Vt = linalg.svd(img_gray)
注意 如果你使用的是你自己的图片，这个命令的运行时间可能会比较长，这取决于你的图片大小和硬件配置。不用担心，这是正常的！奇异值分解 (SVD) 计算可能比较密集。
让我们检查一下结果是否符合预期。
U.shape, s.shape, Vt.shape
((768, 768), (768,), (1024, 1024))
注意 s 具有特定的形状：它只有一维。这意味着某些期望二维数组的线性代数函数可能无法正常工作。例如，根据理论，人们可能期望 s 和 Vt 能够进行矩阵乘法。但是，事实并非如此，因为 s 没有第二维。
s @ Vt
执行上述操作会导致 ValueError 错误。出现这种情况是因为在这种情况下，使用一维数组表示 s 在实践中比构建具有相同数据的对角矩阵更经济。为了重建原始矩阵，我们可以使用 s 的元素重建对角矩阵 \(\Sigma\)，并使其具有合适的维度进行矩阵乘法：在本例中，\(\Sigma\) 应该为 768x1024，因为 U 是 768x768，而 Vt 是 1024x1024。为了将奇异值添加到 Sigma 的对角线上，我们将使用 NumPy 的 fill_diagonal 函数。
import numpy as np
Sigma = np.zeros((U.shape[1], Vt.shape[0]))
np.fill_diagonal(Sigma, s)
现在，我们要检查重建的 U @ Sigma @ Vt 是否接近原始 img_gray 矩阵。
近似值#
linalg 模块包含一个 norm 函数，该函数计算用 NumPy 数组表示的向量或矩阵的范数。例如，根据上面的 SVD 解释，我们预计 img_gray 与重建的 SVD 乘积之间的差的范数应该很小。正如预期的那样，你应该看到类似以下的结果：
linalg.norm(img_gray - U @ Sigma @ Vt)
np.float64(1.43712046073728e-12)
（此操作的实际结果可能因你的架构和线性代数设置而异。无论如何，你都应该看到一个很小的数字。）
我们还可以使用 numpy.allclose 函数来确保重建的乘积实际上与我们的原始矩阵“接近”（两个数组之间的差异很小）。
np.allclose(img_gray, U @ Sigma @ Vt)
True
要查看近似值是否合理，我们可以检查 s 中的值。
plt.plot(s)
plt.show()
在图表中，我们可以看到，尽管我们在 s 中有 768 个奇异值，但大多数奇异值（大约从第 150 个开始）都非常小。因此，只使用与前 (例如，50 个) 个奇异值相关的信息来构建图像的更经济的近似值可能是有意义的。
其思想是将 Sigma 中除前 k 个奇异值（与 s 中的奇异值相同）之外的所有奇异值都视为零，保持 U 和 Vt 不变，并将这些矩阵的乘积计算为近似值。
例如，如果我们选择
k = 10
我们可以通过执行以下操作来构建近似值：
approx = U @ Sigma[:, :k] @ Vt[:k, :]
请注意，我们只需要使用 Vt 的前 k 行，因为其他所有行都将与对应于我们从该近似值中消除的奇异值的零相乘。
plt.imshow(approx, cmap="gray")
plt.show()
现在，你可以继续使用其他 k 值重复此实验，并且根据你选择的 k 值，你的每次实验都应该得到略好（或略差）的图像。
应用于所有颜色#
现在，我们想要执行相同的操作，但要对所有三种颜色都进行操作。我们的第一直觉可能是将我们在上面对每个颜色矩阵分别执行的操作重复一遍。但是，NumPy 的广播功能可以为我们处理这个问题。
如果我们的数组具有两个以上的维度，则可以一次性将 SVD 应用于所有轴。但是，NumPy 中的线性代数函数期望看到一个 (n, M, N) 形式的数组，其中第一个轴 n 表示堆栈中 MxN 矩阵的数量。
在我们的例子中，
img_array.shape
(768, 1024, 3)
因此，我们需要对这个数组的轴进行置换，以获得类似 (3, 768, 1024) 的形状。幸运的是，numpy.transpose 函数可以为我们做到这一点。
np.transpose(x, axes=(i, j, k))
表示轴将被重新排序，以便转置数组的最终形状将根据索引 (i, j, k) 进行重新排序。
让我们看看这对我们的数组有什么影响。
img_array_transposed = np.transpose(img_array, (2, 0, 1))
img_array_transposed.shape
(3, 768, 1024)
现在，我们可以应用 SVD 了。
U, s, Vt = linalg.svd(img_array_transposed)
最后，为了获得完整的近似图像，我们需要将这些矩阵重新组合到近似值中。现在，请注意
U.shape, s.shape, Vt.shape
((3, 768, 768), (3, 768), (3, 1024, 1024))
为了构建最终的近似矩阵，我们必须了解跨不同轴的乘法是如何工作的。
多维数组的乘积#
如果你以前只在 NumPy 中使用过一维或二维数组，你可能会互换使用 numpy.dot 和 numpy.matmul（或 @ 运算符）。但是，对于多维数组，它们的工作方式大相径庭。有关更多详细信息，请查看有关 numpy.matmul 的文档。
现在，为了构建我们的近似值，我们首先需要确保我们的奇异值已准备好进行矩阵乘法，因此我们像之前一样构建 Sigma 矩阵。Sigma 数组的维度必须为 (3, 768, 1024)。为了将奇异值添加到 Sigma 的对角线上，我们将再次使用 fill_diagonal 函数，使用 s 中的 3 行中的每一行作为 Sigma 中 3 个矩阵中的每一个的对角线。
Sigma = np.zeros((3, 768, 1024))
for j in range(3):
np.fill_diagonal(Sigma[j, :, :], s[j, :])
现在，如果我们希望重建完整的 SVD（没有近似值），我们可以执行以下操作：
reconstructed = U @ Sigma @ Vt
请注意
reconstructed.shape
(3, 768, 1024)
重建后的图像应该与原始图像无法区分，除了由于重建过程中的浮点错误造成的差异。回想一下，我们的原始图像由范围为 [0., 1.] 的浮点值组成。重建过程中的浮点错误累积可能会导致值略微超出此原始范围。
reconstructed.min(), reconstructed.max()
(np.float64(-5.558487697898684e-15), np.float64(1.0000000000000053))
由于 imshow 期望范围内的值，我们可以使用 clip 来去除浮点误差。
reconstructed = np.clip(reconstructed, 0, 1)
plt.imshow(np.transpose(reconstructed, (1, 2, 0)))
plt.show()
事实上，imshow 在内部执行此剪切操作，因此，如果你跳过上一个代码单元格中的第一行，你可能会看到一条警告消息，提示 "Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers)."。
现在，为了进行近似，我们必须为每个颜色通道只选择前 k 个奇异值。这可以使用以下语法完成：
approx_img = U @ Sigma[..., :k] @ Vt[..., :k, :]
你可以看到，我们只选择了 Sigma 的最后一个轴的前 k 个分量（这意味着我们只使用了堆栈中三个矩阵中的每一个的前 k 列），并且我们只选择了 Vt 的倒数第二个轴的前 k 个分量（这意味着我们从堆栈 Vt 中的每个矩阵中选择了前 k 行和所有列）。如果你不熟悉省略号语法，它是一个用于其他轴的占位符。有关更多详细信息，请参阅有关 索引 的文档。
现在，
approx_img.shape
(3, 768, 1024)
这对于显示图像来说不是正确的形状。最后，将轴重新排序回我们原始的 (768, 1024, 3) 形状，我们可以看到我们的近似值。
plt.imshow(np.transpose(approx_img, (1, 2, 0)))
plt.show()
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-0.135671314724608..1.079453607915587].
即使图像不那么清晰，使用少量k奇异值（与原始的768个值相比），我们仍然可以恢复图像中的许多显著特征。
结语#
当然，这不是逼近图像的最佳方法。然而，线性代数中确实有一个结果表明，我们上面构建的逼近是根据差的范数而言，我们所能得到的与原始矩阵最接近的逼近。更多信息，请参见G. H. Golub和C. F. Van Loan，《矩阵计算》，马里兰州巴尔的摩，约翰·霍普金斯大学出版社，1985年。