NumPy n 维数组上的线性代数

本篇教程适用于对 NumPy 中的线性代数和数组有基本了解并想了解 n 维 (n>=2) 数组被表示并且可以被操作。特别是，如果您不知道如何将常用函数应用于 n 维数组（不使用 for 循环），或者如果您想了解 n 维数组的轴和形状属性，本教程可能会有所帮助。

学习目标

完成本教程后，您应该能够：

• 理解NumPy中一维、二维和n维数组的区别；
• 了解如何在不使用 for 循环的情况下将一些线性代数运算应用于 n 维数组；
• 了解 n 维数组的轴和形状属性。

内容

在本教程中，我们将使用线性代数的矩阵分解，即奇异值分解，来生成图像的压缩近似值。我们将使用模块中的face图像scipy.misc：

>>> from scipy import misc
>>> img = misc.face()

笔记
如果您愿意，可以在学习本教程时使用自己的图像。为了将您的图像转换为可操作的 NumPy 数组，您可以使用子模块中的imread函数matplotlib.pyplot。另外，您也可以使用imageio.imread从函数 imageio库。请注意，如果您使用自己的图像，则可能需要调整以下步骤。有关如何，当转换成NumPy的阵列的图像处理的详细信息，请参阅上NumPy的甲速成为图像从所述scikit-image文档。

现在，img是一个 NumPy 数组，正如我们在使用type函数时所看到的：

>>> type(img)
<class 'numpy.ndarray'>

我们可以使用以下matplotlib.pyplot.imshow函数查看图像：

>>> import matplotlib.pyplot as plt
>>> plt.imshow(img)

笔记
如果您在 IPython shell 中执行上述命令，则可能需要使用该命令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 的每个值img[:,:,0]都是 0 到 255 之间的整数值，代表每个对应图像像素中的红色级别（请记住，如果您使用自己的图像而不是 ，这可能会有所不同scipy.misc.face）。

正如预期的那样，这是一个 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 模块将线性代数运算应用于 n 维数组。有关这方面的更多信息，请查看 scipy.linalg 参考。

要继续，请从 NumPy 导入线性代数子模块：

>>> from numpy import linalg

为了从给定的矩阵中提取信息，我们可以使用 SVD 获得 3 个数组，这些数组可以相乘以获得原始矩阵。从线性代数理论，给定一个矩阵A，可以计算以下乘积： UΣVT=A

在哪里 U 和 VT 是正方形和 A. Σ 是对角矩阵，并且包含 的奇异值的 A 从大到小排列。这些值总是非负的，可以用作矩阵表示的某些特征的“重要性”的指标 A. 让我们先看看这在实践中是如何使用一个矩阵的。请注意，根据色度学，如果我们应用公式，则可以获得我们彩色图像的相当合理的灰度版本。 Y=0.2126R+0.7152G+0.0722B 在哪里 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")

现在，将该linalg.svd函数应用于该矩阵，我们得到以下分解：

>>> U, s, Vt = linalg.svd(img_gray)

笔记
如果您使用自己的映像，则此命令可能需要一段时间才能运行，具体取决于您的映像大小和硬件。别担心，这是正常的！SVD 可能是一个非常密集的计算。如果您使用自己的映像，则此命令可能需要一段时间才能运行，具体取决于您的映像大小和硬件。别担心，这是正常的！SVD 可能是一个非常密集的计算。

让我们检查一下这是否是我们所期望的：

>>> U.shape, s.shape, Vt.shape
((768, 768), (768,), (1024, 1024))

请注意，s它具有特定的形状：它只有一个维度。这意味着某些需要二维数组的线性代数函数可能不起作用。例如，从理论来看，人们可能期望乘法s并Vt与之相容。然而，这是不正确的，因为s它没有第二个轴。执行

>>> s @ Vt
Traceback (most recent call last):
  ...
ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0,
with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 1024 is different from
768)

结果是ValueError. 发生这种情况是因为s在这种情况下，对于 的一维数组在实践中比使用相同的数据构建对角矩阵要经济得多。为了重建原始矩阵，我们可以重建对角矩阵Σs其对角线上的元素和适当的乘法维度：在我们的例子中，Σ应该U是 768x1024，因为是 768x768 和 1024x1024 Vt。

>>> import numpy as np
>>> Sigma = np.zeros((768, 1024))
>>> for i in range(768):
...     Sigma[i, i] = s[i]

现在，我们要检查重建的矩阵是否接近原始矩阵。U @ Sigma @ Vt``img_gray

近似值

该linalg模块包括一个norm函数，该函数计算以 NumPy 数组表示的向量或矩阵的范数。例如，从上面的 SVD 解释中，我们期望img_gray重构的 SVD 乘积之间的差异范数很小。正如预期的那样，你应该看到类似的东西

>>> linalg.norm(img_gray - U @ Sigma @ Vt)
1.3926466851808837e-12

（此操作的实际结果可能会因您的体系结构和线性代数设置而异。无论如何，您应该看到一个小数字。） 我们也可以使用该numpy.allclose函数来确保重构的乘积实际上接近我们的原始矩阵（两个数组之间的差异很小）：

>>> np.allclose(img_gray, U @ Sigma @ Vt)
True

要查看近似值是否合理，我们可以检查 中的值s：

>>> plt.plot(s)

在图中，我们可以看到，虽然我们在 中有 768 个奇异值 s，但大多数（在第 150 个条目之后）都非常小。因此，仅使用与第一个（例如 50 个）奇异值相关的信息来构建更经济的图像近似值可能是有意义的 。

这个想法是将除了第一个k奇异值之外的 所有Sigma（与 相同s）视为零，保持 U和Vt完整，并计算这些矩阵的乘积作为近似值。 例如，如果我们选择

>>> k = 10

我们可以通过做来建立近似

>>> approx = U @ Sigma[:, :k] @ Vt[:k, :]

请注意，我们必须仅使用 的第一k行Vt，因为所有其他行都将乘以对应于我们从该近似中消除的奇异值的零。

>>> plt.imshow(approx, cmap="gray")

现在，您可以继续使用k 的其他值重复此实验，并且您的每个实验都应该根据您选择的值为您提供稍好（或更差）的图像。

适用于所有颜色

现在我们想要做同样的操作，但是对所有三种颜色。我们的第一直觉可能是对每个颜色矩阵分别重复我们上面所做的相同操作。然而，NumPy 的广播为我们处理了这个问题。

如果我们的数组有两个以上的维度，那么 SVD 可以一次应用于所有轴。但是，NumPy 中的线性代数函数期望看到形式为 的数组，其中第一个轴表示矩阵的数量。(N, M, M)

在我们的例子中，

>>> 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.dotand numpy.matmul（或@运算符）。但是，对于 n 维数组，它们的工作方式非常不同。有关更多详细信息，请查看文档numpy.matmul。

现在，为了建立我们的近似值，我们首先需要确保我们的奇异值已准备好进行乘法运算，因此我们Sigma以与之前所做的类似的方式建立我们的矩阵。所述Sigma阵列必须具有尺寸 。为了将奇异值添加到 的对角线上 ，我们将使用NumPy 中的函数，使用 3 行中的每一行作为 中 3 个矩阵中的每一个的对角线：(3, 768, 1024)``Sigma``numpy.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)

和

>>> plt.imshow(np.transpose(reconstructed, (1, 2, 0)))

应该为您提供与原始图像无法区分的图像（尽管我们可能会为此重建引入浮点错误）。事实上，您可能会看到一条警告消息说“将输入数据裁剪到有效范围以使用 RGB 数据进行 imshow（[0..1] 表示浮点数或 [0..255] 表示整数）。” 这是我们刚刚对原始图像进行的操作所预期的。

现在，要进行近似，我们必须仅为k每个颜色通道选择第一个奇异值。这可以使用以下语法完成：

>>> approx_img = U @ Sigma[..., :k] @ Vt[..., :k, :]

您可以看到我们只选择k了最后一个轴的第一个组件Sigma（这意味着我们只使用k了堆栈中三个矩阵中每一个的第一列），并且我们只选择k了第二个中的第一个组件-to-last 轴Vt（这意味着我们只k从堆栈中的每个矩阵Vt和所有列中选择了第一行）。如果您不熟悉省略号语法，则它是其他轴的占位符。有关更多详细信息，请参阅有关Indexing的文档 。

现在，

>>> approx_img.shape
(3, 768, 1024)

这不是显示图像的正确形状。最后，将轴重新排序回我们的原始形状，我们可以看到我们的近似值：(768, 1024, 3)

>>> plt.imshow(np.transpose(approx_img, (1, 2, 0)))

即使图像不那么清晰，使用少量k奇异值（与原始的 768 个值集相比），我们可以从该图像中恢复许多区别特征。

总结

当然，这不是近似图像的最佳方法。然而，事实上，线性代数中的一个结果表明，我们上面建立的近似值是我们在差分范数方面可以得到的最好的原始矩阵。有关更多信息，请参阅GH Golub 和 CF Van Loan，矩阵计算，巴尔的摩，医学博士，约翰霍普金斯大学出版社，1985 年。