教程：n维数组上的线性代数¶

先决条件¶

在阅读本教程之前，您应该了解一点Python。如果你想刷新你的记忆，看看 Python tutorial .

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

学习者简介¶

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

学习目标¶

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

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

内容¶

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

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

注解

如果您愿意，您可以在完成本教程时使用自己的图像。为了将图像转换为可以操纵的NumPy数组，可以使用 imread 函数来自 matplotlib.pyplot 子模块。或者，您可以使用 imageio.imread 函数来自 imageio 类库。请注意，如果您使用自己的图像，您可能需要调整以下步骤。有关转换为NumPy数组时如何处理图像的详细信息，请参见 A crash course on NumPy for images 从 scikit-image 文档。

现在， img 是一个NumPy数组，正如我们在使用 type 功能：

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

我们可以使用 matplotlib.pyplot.imshow 功能：

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

注解

如果在ipythonshell中执行上述命令，则可能需要使用该命令 plt.show() 显示图像窗口。

形状、轴和阵列属性¶

注意，在线性代数中，向量的维数是指数组中的项数。而不是NumPy定义轴数。例如，一维数组是一个向量，如 [1, 2, 3] ，二维数组是矩阵，依此类推。

首先，让我们检查数组中数据的形状。由于此图像是二维的（图像中的像素形成一个矩形），因此我们可能希望使用二维数组来表示它（矩阵）。但是，使用 shape 这个NumPy数组的属性给出了一个不同的结果：

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

输出是一个 tuple 有三个元素，这意味着这是一个三维数组。事实上，由于这是一个彩色图像，我们已经使用了 imread 函数读取它时，数据被组织在三个2D数组中，表示颜色通道（在本例中，为红色、绿色和蓝色-RGB）。您可以通过查看上面的形状看到这一点：它表示我们有一个由3个矩阵组成的数组，每个矩阵的形状为768x1024。

此外，使用 ndim 这个数组的属性，我们可以看到

>>> img.ndim
3

NumPy将每个维度称为 axis . 因为什么 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.misc.face ）

正如所料，这是一个768x1024矩阵：

>>> img[:, :, 0].shape
(768, 1024)

因为我们要对这些数据执行线性代数运算，所以在矩阵的每个条目中使用0到1之间的实数来表示RGB值可能更有趣。我们可以通过设置

>>> img_array = img / 255

此操作将数组除以标量，因为NumPy的 broadcasting rules ). （请注意，在实际应用中，最好使用 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 (Singular Value Decomposition) 尝试重建比原始图像使用较少奇异值信息的图像，同时仍保留其某些特征。

注解

我们将使用NumPy的线性代数模块， numpy.linalg ，以执行本教程中的操作。本模块中的大部分线性代数函数也可以在中找到 scipy.linalg ，并鼓励用户使用 scipy 实际应用程序模块。然而，目前不可能将线性代数运算应用于使用 scipy.linalg 模块。有关详细信息，请查看 scipy.linalg Reference .

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

>>> from numpy import linalg

为了从一个给定的矩阵中提取信息，我们可以使用奇异值分解得到3个数组，这些数组可以相乘得到原始矩阵。从线性代数理论出发，给出了一个矩阵 $A$ ，可计算以下乘积：

$U\Sigma V^T=A$

在哪里？ $U$ 和 $V^T$ 是方形的 $\Sigma$ 大小与 $A$ . $\Sigma$ 是对角矩阵，包含 singular values 属于 $A$ ，从大到小排列。这些值总是非负的，可以作为矩阵表示的某些特征的“重要性”的指标 $A$ .

让我们先来看看这是如何在一个矩阵的情况下工作的。注意，根据 colorimetry 如果我们应用这个公式，就有可能得到一个相当合理的彩色图像的灰度版本

System Message: WARNING/2 (Y=0.2126 R+0.7152克+0.0722硼)

latex exited with error [stdout] This is pdfTeX, Version 3.14159265-2.6-1.40.19 (TeX Live 2019/dev/Debian) (preloaded format=latex) restricted \write18 enabled. entering extended mode (./math.tex LaTeX2e <2018-12-01> (/usr/share/texlive/texmf-dist/tex/latex/base/article.cls Document Class: article 2018/09/03 v1.4i Standard LaTeX document class (/usr/share/texlive/texmf-dist/tex/latex/base/size12.clo)) (/usr/share/texlive/texmf-dist/tex/latex/base/inputenc.sty) (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsmath.sty For additional information on amsmath, use the `?' option. (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amstext.sty (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsgen.sty)) (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsbsy.sty) (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsopn.sty)) (/usr/share/texlive/texmf-dist/tex/latex/amscls/amsthm.sty) (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/amssymb.sty (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/amsfonts.sty)) (/usr/share/texlive/texmf-dist/tex/latex/anyfontsize/anyfontsize.sty) (/usr/share/texlive/texmf-dist/tex/latex/tools/bm.sty) (./math.aux) (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/umsa.fd) (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/umsb.fd) ! Package inputenc Error: Unicode character 克 (U+514B) (inputenc) not set up for use with LaTeX. See the inputenc package documentation for explanation. Type H <return> for immediate help. ... l.14 ...}Y=0.2126 R+0.7152克+0.0722硼\end{split} ! Package inputenc Error: Unicode character 硼 (U+787C) (inputenc) not set up for use with LaTeX. See the inputenc package documentation for explanation. Type H <return> for immediate help. ... l.14 ...}Y=0.2126 R+0.7152克+0.0722硼\end{split} ! Package inputenc Error: Unicode character 克 (U+514B) (inputenc) not set up for use with LaTeX. See the inputenc package documentation for explanation. Type H <return> for immediate help. ... l.14 ...}Y=0.2126 R+0.7152克+0.0722硼\end{split} ! Package inputenc Error: Unicode character 硼 (U+787C) (inputenc) not set up for use with LaTeX. See the inputenc package documentation for explanation. Type H <return> for immediate help. ... l.14 ...}Y=0.2126 R+0.7152克+0.0722硼\end{split} [1] (./math.aux) ) (see the transcript file for additional information) Output written on math.dvi (1 page, 328 bytes). Transcript written on math.log.

在哪里？ $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 将默认为与实际数据不对应的颜色映射）。

在我们的例子中，我们近似于图像的灰度部分，所以我们将使用colormap gray ：：

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

现在，应用 linalg.svd 函数，我们得到如下分解：

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

注解

如果您正在使用自己的映像，根据映像和硬件的大小，此命令可能需要一段时间才能运行。别担心，这很正常！奇异值分解是一个相当密集的计算。

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

>>> 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 在这种情况下，在实践中比用相同的数据建立对角矩阵要经济得多。为了重建原矩阵，我们可以重建对角矩阵 $\Sigma$ 具有 s 在对角线和适当的尺寸相乘：在我们的例子中， $\Sigma$ 应为768x1024，因为 U is 768x768和 Vt 为1024x1024。

>>> 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 重构后的奇异值分解产物较小。正如所料，您应该看到

>>> 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的 broadcasting 帮我们处理这个。

如果我们的阵列有两个以上的维度，那么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.dot 和 numpy.matmul （或 @ （操作员）可互换。然而，对于n维数组，它们的工作方式截然不同。有关更多详细信息，请查看文档 numpy.matmul .

现在，为了建立近似，我们首先需要确保我们的奇异值可以乘法，所以我们建立 Sigma 矩阵类似于我们以前做的。这个 Sigma 数组必须有维度 (3, 768, 1024) . 为了将奇异值加到 Sigma ，我们将使用 fill_diagonal 函数，使用 s 作为中3个矩阵的对角线 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)))

应该会给你一个与原始图像无法区分的图像（尽管我们可能会引入浮点错误进行重建）。事实上，你可能会看到一条警告信息说 "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 从第二个轴到最后一个轴的组件 Vt （这意味着我们只选择了第一个 k 堆栈中每个矩阵的行 Vt 和所有列）。如果不熟悉省略号语法，则它是其他轴的占位符。有关更多详细信息，请参阅上的文档 Indexing .

现在，

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

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

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

即使图像不是那么清晰，使用少量的 k 奇异值（与原始的一组768个值相比），我们可以从这幅图像中恢复许多显著特征。

最后一句话¶

当然，这不是解决问题的最佳方法近似一个形象。然而，事实上，在线性代数中，有一个结果表明，我们在上面建立的近似是我们能够得到的关于差模的原始矩阵的最佳近似。有关详细信息，请参阅 G、 H.Golub和C.F.Van Loan，《矩阵计算》，马里兰州巴尔的摩，约翰霍普金斯大学出版社，1985年 .

进一步阅读¶

NumPy教程教程：遮罩阵列