NumPy 特性：n维数组上的线性代数

先决条件

在阅读本教程之前，您应该对 Python 有一定的了解。如果您想恢复记忆，请参考 Python 教程。

如果您想要运行本教程中的示例，您还应该在计算机上安装 matplotlib 和 SciPy。

学习对象

本教程适用于对线性代数和 NumPy 中的数组有基本了解，并希望了解如何表示和操作 n 维数组的人。特别是，如果您不知道如何将常见函数应用于 n 维数组（而不使用 for 循环），或者如果您想了解 n 维数组的轴和形状属性，那么本教程可能会对您有所帮助。

学习目标

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

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

内容

在本教程中，我们将使用线性代数中的矩阵分解方法——奇异值分解（Singular Value Decomposition，简称 SVD），生成图像的压缩近似。我们将使用 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()

从 'https://raw.githubusercontent.com/scipy/dataset-face/main/face.dat' 下载文件 'face.dat' 到 '/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 之间的实数可能更有趣。我们可以通过设置

img_array = img / 255

来实现这一点。这个操作是将一个数组除以一个标量，它能够工作是因为 NumPy 的广播规则（请注意，在实际应用中，最好使用 scikit-image 中的 img_as_float 实用函数）。

您可以通过进行一些测试来验证上述操作是否有效；例如，查询该数组的最大值和最小值：

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 个数组，这些数组可以相乘以获得原始矩阵。根据线性代数的理论，给定一个矩阵 ，可以计算以下乘积：

其中和 是方阵，大小与 相同。是一个对角矩阵，包含按从大到小排列的 奇异值。这些值始终是非负的，并且可以用作表示矩阵所表示的某些特征的“重要性”的指标。

让我们先看看如何在实践中使用一个矩阵。请注意，根据 色度学，如果我们应用以下公式，可以获得我们彩色图像的相当合理的灰度版本：

其中 是表示灰度图像的数组，和 是我们最初的红色、绿色和蓝色通道数组。请注意，我们可以使用 @ 运算符（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()

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

注意 如果您使用自己的图像，则此命令可能需要一些时间才能运行，具体取决于图像的大小和硬件配置。不用担心，这是正常的！SVD 可能是一项相当密集的计算。

让我们检查一下这是否是我们预期的结果：

U.shape, s.shape, Vt.shape

((768, 768), (768,), (1024, 1024))

请注意，s 具有特定的形状：它只有一个维度。这意味着一些期望 2D 数组的线性代数函数可能无法正常工作。例如，从理论上讲，人们可能期望 s 和 Vt 可以相乘。然而，这是不正确的，因为 s 没有第二个轴。执行

s @ Vt

会导致 ValueError。这是因为在实践中，使用一个一维数组来表示 s 要比使用相同数据构建对角矩阵更经济。为了重建原始矩阵，我们可以使用 s 的元素重新构建对角矩阵，并具有适当的维度进行乘法运算：在我们的例子中，应该是 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)

1.43712046073728e-12

（此操作的实际结果可能因您的架构和线性代数设置而异。无论如何，您应该看到一个很小的数字。）

我们也可以使用 numpy.allclose 函数来确保重建的乘积实际上与我们的原始矩阵“接近”（两个数组之间的差异很小）：

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

True

为了判断近似是否合理，我们可以检查 s 中的值：

plt.plot(s)
plt.show()

从图中可以看出，尽管 s 中有 768 个奇异值，但大多数（大约从第 150 个开始）都非常小。因此，使用与前 50 个奇异值相关的信息来构建对图像的更经济的近似可能是有意义的。

我们的想法是将除了前 k 个奇异值之外的所有奇异值都视为零，保持 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()

(-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, :]

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

现在，

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 个值相比），我们可以恢复出许多图像的特征。文章来源地址https://www.toymoban.com/news/detail-841093.html

最后的话

进一步阅读

• Python 教程
• NumPy 参考手册
• SciPy 教程
• SciPy 讲义
• 一个 Matlab、R、IDL、NumPy/SciPy 字典

到了这里，关于NumPy 特性：n维数组上的线性代数的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！