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()
注意：如果你愿意，可以在学习本教程时使用自己的图片。为了将你的图片转换成可以操作的 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中，它定义的是轴的数量。例如，1D数组是一个向量，如[1, 2, 3]，2D数组是一个矩阵，依此类推。
首先，让我们检查数组中数据的形状。由于这张图像是二维的（图像中的像素形成一个矩形），我们可能会期望一个二维数组来表示它（一个矩阵）。然而，使用这个NumPy数组的shape属性给了我们一个不同的结果：
img.shape
(768, 1024, 3)
输出是一个包含三个元素的 tuple，这意味着这是一个三维数组。实际上，由于这是一张彩色图像，并且我们使用了 imread 函数来读取它，数据被组织成三个 2D 数组，分别代表颜色通道（在这种情况下，红色、绿色和蓝色 - RGB）。你可以通过查看上面的形状来看到这一点：它表明我们有一个包含 3 个矩阵的数组，每个矩阵的形状为 768x1024。
此外，使用这个数组的 ndim 属性，我们可以看到
img.ndim
3
NumPy 将每个维度称为一个 轴 。由于 imread 的工作方式，第3轴中的第一个索引 是我们图像的红色像素数据。我们可以通过使用语法来访问它。
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 函数，仅支持 2D 数组。有关更多信息，请查看 scipy.linalg 页面。
要继续，请从NumPy导入线性代数子模块：
from numpy import linalg
为了从一个给定的矩阵中提取信息，我们可以使用SVD来获得3个数组，这些数组可以相乘以获得原始矩阵。根据线性代数的理论，给定一个矩阵 \(A\)，可以计算以下乘积：
其中 \(U\) 和 \(V^T\) 是方阵，\(\Sigma\) 的大小与 \(A\) 相同。\(\Sigma\) 是一个对角矩阵，包含 \(A\) 的奇异值，从大到小排列。这些值总是非负的，可以作为矩阵 \(A\) 所表示的某些特征“重要性”的指标。
让我们先看看这在实践中如何使用一个矩阵。请注意，根据 colorimetry，如果我们应用公式，可以得到一个相当合理的彩色图像的灰度版本。
其中 \(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.6162976380358187e-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 值重复这个实验，并且每次实验应该会根据你选择的值给你一个稍微更好（或更差）的图像。
应用于所有颜色#
现在我们想对所有三种颜色进行同样的操作。我们的第一反应可能是对每个颜色矩阵单独重复上述操作。然而，NumPy 的 broadcasting 为我们处理了这个问题。
如果我们的数组有超过两个维度，那么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))
要构建最终的近似矩阵，我们必须理解如何在不同轴上进行乘法。
具有 n 维数组的产品#
如果你之前只使用过 NumPy 中的一维或二维数组，你可能会交替使用 numpy.dot 和 numpy.matmul（或 @ 运算符）。然而，对于 n 维数组，它们的工作方式非常不同。更多详情，请查看 numpy.matmul 的文档。
现在，为了构建我们的近似值，我们首先需要确保我们的奇异值已准备好进行乘法运算，因此我们以之前的方式构建我们的 Sigma 矩阵。Sigma 数组必须具有维度 (3, 768, 1024)。为了将奇异值添加到 Sigma 的对角线上，我们将再次使用 fill_diagonal 函数，使用 s 中的每一行作为 Sigma 中每个矩阵的对角线：
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(-4.330520317341602e-15), np.float64(1.0000000000000069))
由于 imshow 期望值在范围内，我们可以使用 clip 来切除浮点误差：
reconstructed = np.clip(reconstructed, 0, 1)
plt.imshow(np.transpose(reconstructed, (1, 2, 0)))
plt.show()
事实上，imshow 在幕后执行这种裁剪，所以如果你跳过前面代码单元中的第一行，你可能会看到一条警告信息，说 "裁剪输入数据到 imshow 的有效范围，对于 RGB 数据（浮点数为 [0..1]，整数为 [0..255]）。"
现在，为了进行近似，我们必须仅为每个颜色通道选择前 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.13567131472460894..1.0794536079155879].
尽管图像不如原来清晰，使用少量 k 奇异值（与原始的 768 个值相比），我们可以恢复图像中的许多显著特征。
结束语#
当然，这不是近似图像的最佳方法。然而，实际上，线性代数中有一个结果表明，我们上面构建的近似是我们在范数差异方面能得到的对原始矩阵的最佳近似。更多信息，请参见G. H. Golub 和 C. F. Van Loan，Matrix Computations，Baltimore，MD，Johns Hopkins University Press，1985。