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.imread 函数，该函数来自 imageio 库。请注意，如果您使用自己的图像，您可能需要调整下面的步骤。有关图像在转换为 NumPy 数组时的处理方式的更多信息，请参阅 NumPy 图像速成课程，该课程来自 scikit-image 文档。
现在，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 值重复此实验，并且根据您选择的数值，每个实验都应该得到略微好（或差）的图像。
应用于所有颜色#
现在我们想做相同类型操作，但要对所有三种颜色进行操作。我们第一个直觉可能是对每个颜色矩阵分别重复我们上面所做的操作。但是，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))
为了构建最终的近似矩阵，我们必须了解跨不同轴的乘法是如何工作的。
与 n 维数组的乘积#
如果您之前在 NumPy 中只使用过一维或二维数组，您可能会交替使用 numpy.dot 和 numpy.matmul（或 @ 运算符）。但是，对于 n 维数组，它们以非常不同的方式工作。有关更多详细信息，请查看有关 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。