n 维数组上的线性代数#
先决条件#
在阅读本教程之前，您应该了解一些 Python。如果您想加深记忆，请查看Python 教程。
如果您希望能够运行本教程中的示例，您还应该在计算机上安装matplotlib和SciPy 。
学习者简介#
本教程面向对 NumPy 中的线性代数和数组有基本了解并想要了解 n 维 (\(n>=2\)) 数组被表示并且可以被操作。特别是，如果您不知道如何将常用函数应用于 n 维数组（不使用 for 循环），或者如果您想了解 n 维数组的轴和形状属性，本教程可能会有所帮助。
学习目标＃
学完本教程后，您应该能够：
了解 NumPy 中一维、二维和 n 维数组的区别；
了解如何在不使用 for 循环的情况下将一些线性代数运算应用于 n 维数组；
了解 n 维数组的轴和形状属性。
内容＃
在本教程中，我们将使用线性代数的矩阵分解（奇异值分解）来生成图像的压缩近似值。我们将使用scipy.datasetsface模块中的图像：
# 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.pyplotimread子模块中的函数。或者，您可以使用库中的imageio.imread函数。请注意，如果您使用自己的图像，则可能需要调整以下步骤。有关转换为 NumPy 数组时如何处理图像的更多信息，请参阅文档中的NumPy 图像速成课程。imageioscikit-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]
首先，让我们检查数组中数据的形状。由于该图像是二维的（图像中的像素形成一个矩形），因此我们可能期望用二维数组来表示它（矩阵）。然而，使用shape这个 NumPy 数组的属性给我们带来了不同的结果：
img.shape
(768, 1024, 3)
输出是一个包含三个元素的元组，这意味着这是一个三维数组。事实上，由于这是一个彩色图像，并且我们使用该imread函数来读取它，因此数据被组织在三个 2D 数组中，代表颜色通道（在本例中为红色、绿色和蓝色 - 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)
从上面的输出中，我们可以看到 in 中的每个值都是 0 到 255 之间的整数值，表示每个相应图像像素中的红色级别（请记住，如果您使用自己的图像而不是scipy.image，这可能会有所不同）。数据集.face）。img[:, :, 0]
正如预期的那样，这是一个 768x1024 矩阵：
img[:, :, 0].shape
(768, 1024)
由于我们要对此数据执行线性代数运算，因此在矩阵的每个条目中使用 0 到 1 之间的实数来表示 RGB 值可能会更有趣。我们可以通过设置来做到这一点
img_array = img / 255
此操作（将数组除以标量）之所以有效，是因为 NumPy 的广播规则。 （请注意，在实际应用中，最好使用例如中的img_as_float实用函数scikit-image）。
您可以通过做一些测试来检查上述内容是否有效；例如，查询该数组的最大值和最小值：
img_array.max(), img_array.min()
(1.0, 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，在实践中比使用相同数据构建对角矩阵要经济得多。为了重建原始矩阵，我们可以重建对角矩阵\(\Sigma\)其对角线中的元素s以及用于相乘的适当维度：在我们的例子中，\(\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 @ Vtimg_gray
近似值#
linalg模块包含一个norm函数，用于计算 NumPy 数组中表示的向量或矩阵的范数。例如，根据上面的 SVD 解释，我们预计img_graySVD 乘积与重构的 SVD 乘积之间的差异范数会很小。正如预期的那样，您应该看到类似的内容
linalg.norm(img_gray - U @ Sigma @ Vt)
1.43712046073728e-12
（此操作的实际结果可能会有所不同，具体取决于您的体系结构和线性代数设置。无论如何，您应该看到一个很小的数字。）
我们还可以使用numpy.allclose函数来确保重建的乘积实际上接近我们的原始矩阵（两个数组之间的差异很小）：
np.allclose(img_gray, U @ Sigma @ Vt)
True
要查看近似值是否合理，我们可以检查以下值s：
plt.plot(s)
plt.show()
在图中，我们可以看到，尽管 中有 768 个奇异值s，但其中大多数（在第 150 个条目之后）都非常小。因此，仅使用与第一个（例如 50 个）奇异值相关的信息来构建更经济的图像近似值可能是有意义的。
这个想法是将除了第一个k奇异值Sigma（与 中相同s）之外的所有奇异值视为零，保持U完整Vt，并计算这些矩阵的乘积作为近似值。
例如，如果我们选择
k = 10
我们可以通过这样做来建立近似值
approx = U @ Sigma[:, :k] @ Vt[:k, :]
请注意，我们必须仅使用k的第一行Vt，因为所有其他行都将乘以与我们从此近似中消除的奇异值相对应的零。
plt.imshow(approx, cmap="gray")
plt.show()
现在，您可以继续使用 的其他值重复此实验k，并且每个实验都应该根据您选择的值提供稍好（或更差）的图像。
适用于所有颜色#
现在我们想做同样的操作，但是对所有三种颜色。我们的第一直觉可能是对每个颜色矩阵分别重复上面所做的相同操作。然而，NumPy 的广播为我们解决了这个问题。
如果我们的数组超过二维，那么 SVD 可以同时应用于所有轴。但是，NumPy 中的线性代数函数期望看到 形式的数组，其中第一个轴表示堆栈中矩阵的数量。(n, M, N)nMxN
在我们的例子中，
img_array.shape
(768, 1024, 3)
所以我们需要排列这个数组上的轴以获得类似于 的形状。幸运的是，numpy.transpose函数可以为我们做到这一点：(3, 768, 1024)
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必须有维度。为了将奇异值添加到 的对角线上，我们将再次使用fill_diagonal函数，使用 中 3 行中的每一行作为 中 3 个矩阵中每一个的对角线：(3, 768, 1024)SigmasSigma
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()
(-5.558487697898684e-15, 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, :]
您可以看到，我们仅选择了k最后一个轴的第一个分量Sigma（这意味着我们仅使用了k堆栈中三个矩阵中每个矩阵的第一列），并且我们仅选择了k第二个轴中的第一个分量-to-last 轴（这意味着我们仅从堆栈中的每个矩阵和所有列中Vt选择了第一行）。如果您不熟悉省略号语法，它是其他轴的占位符。有关更多详细信息，请参阅有关索引的文档。kVt
现在，
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).
尽管图像不那么清晰，但使用少量的k奇异值（与原始的 768 个值相比），我们可以从该图像中恢复许多显着特征。
最后的话#
当然，这不是近似图像的最佳方法。然而，事实上，线性代数中的一个结果表明，我们上面构建的近似值是我们可以根据差值范数得到的原始矩阵的最佳近似值。有关详细信息，请参阅GH Golub 和 CF Van Loan，《矩阵计算》，马里兰州巴尔的摩，约翰霍普金斯大学出版社，1985 年。