NumPy n 維數(shù)組上的線性代數(shù)

本篇教程適用于對 NumPy 中的線性代數(shù)和數(shù)組有基本了解并想了解 n 維 (n>=2) 數(shù)組被表示并且可以被操作。特別是，如果您不知道如何將常用函數(shù)應(yīng)用于 n 維數(shù)組（不使用 for 循環(huán)），或者如果您想了解 n 維數(shù)組的軸和形狀屬性，本教程可能會有所幫助。

學(xué)習(xí)目標(biāo)

完成本教程后，您應(yīng)該能夠：

• 理解NumPy中一維、二維和n維數(shù)組的區(qū)別；
• 了解如何在不使用 for 循環(huán)的情況下將一些線性代數(shù)運(yùn)算應(yīng)用于 n 維數(shù)組；
• 了解 n 維數(shù)組的軸和形狀屬性。

內(nèi)容

在本教程中，我們將使用線性代數(shù)的矩陣分解，即奇異值分解，來生成圖像的壓縮近似值。我們將使用模塊中的face圖像scipy.misc：

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

筆記
如果您愿意，可以在學(xué)習(xí)本教程時使用自己的圖像。為了將您的圖像轉(zhuǎn)換為可操作的 NumPy 數(shù)組，您可以使用子模塊中的imread函數(shù)matplotlib.pyplot。另外，您也可以使用imageio.imread從函數(shù)?imageio庫。請注意，如果您使用自己的圖像，則可能需要調(diào)整以下步驟。有關(guān)如何，當(dāng)轉(zhuǎn)換成NumPy的陣列的圖像處理的詳細(xì)信息，請參閱上NumPy的甲速成為圖像從所述scikit-image文檔。

現(xiàn)在，img是一個 NumPy 數(shù)組，正如我們在使用type函數(shù)時所看到的：

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

我們可以使用以下matplotlib.pyplot.imshow函數(shù)查看圖像：

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

筆記
如果您在 IPython shell 中執(zhí)行上述命令，則可能需要使用該命令plt.show()來顯示圖像窗口。

形狀、軸和數(shù)組屬性

請注意，在線性代數(shù)中，向量的維數(shù)是指數(shù)組中的條目數(shù)。在 NumPy 中，它改為定義軸的數(shù)量。例如，一維數(shù)組是向量，例如，二維數(shù)組是矩陣，等等。[1,?2,?3]

首先，讓我們檢查數(shù)組中數(shù)據(jù)的形狀。由于這個圖像是二維的（圖像中的像素形成一個矩形），我們可能期望一個二維數(shù)組來表示它（一個矩陣）。然而，使用shape?這個 NumPy 數(shù)組的屬性會給我們一個不同的結(jié)果：

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

輸出是一個包含三個元素的元組，這意味著這是一個三維數(shù)組。事實(shí)上，由于這是一張彩色圖像，我們已經(jīng)使用imread函數(shù)讀取它，數(shù)據(jù)被組織成三個 2D 數(shù)組，代表顏色通道（在這種情況下，紅色、綠色和藍(lán)色 - RGB）。您可以通過查看上面的形狀來看到這一點(diǎn)：它表明我們有一個由 3 個矩陣組成的數(shù)組，每個矩陣的形狀為 768x1024。

此外，使用ndim這個數(shù)組的屬性，我們可以看到

>>> img.ndim
3

NumPy 將每個維度稱為軸。由于imread?工作原理，第三軸中的第一個索引是我們圖像的紅色像素數(shù)據(jù)。我們可以使用語法訪問它

>>> 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 之間的整數(shù)值，代表每個對應(yīng)圖像像素中的紅色級別（請記住，如果您使用自己的圖像而不是 ，這可能會有所不同scipy.misc.face）。

正如預(yù)期的那樣，這是一個 768x1024 的矩陣：

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

由于我們將對該數(shù)據(jù)執(zhí)行線性代數(shù)運(yùn)算，因此在矩陣的每個條目中使用 0 到 1 之間的實(shí)數(shù)來表示 RGB 值可能更有趣。我們可以通過設(shè)置來做到這一點(diǎn)

>>> img_array = img / 255

由于 NumPy 的廣播規(guī)則，這個將數(shù)組除以標(biāo)量的操作有效?）。（請注意，在實(shí)際應(yīng)用程序中，最好使用例如img_as_float來自的?效用函數(shù)?scikit-image）。 您可以通過進(jìn)行一些測試來檢查上述內(nèi)容是否有效；例如，查詢這個數(shù)組的最大值和最小值：

>>> img_array.max(), img_array.min()
(1.0, 0.0)

或檢查數(shù)組中的數(shù)據(jù)類型：

>>> img_array.dtype
dtype('float64')

請注意，我們可以使用切片語法將每個顏色通道分配給一個單獨(dú)的矩陣：

>>> red_array = img_array[:, :, 0]
>>> green_array = img_array[:, :, 1]
>>> blue_array = img_array[:, :, 2]

軸上的操作

可以使用線性代數(shù)中的方法來逼近現(xiàn)有數(shù)據(jù)集。在這里，我們將使用SVD（奇異值分解）來嘗試重建一張圖像，該圖像使用比原始圖像更少的奇異值信息，同時仍保留其一些特征。

筆記
我們將使用 NumPy 的線性代數(shù)模塊numpy.linalg來執(zhí)行本教程中的操作。該模塊中的大多數(shù)線性代數(shù)函數(shù)也可以在 中找到scipy.linalg，鼓勵用戶將該scipy模塊用于實(shí)際應(yīng)用。但是，目前無法使用該scipy.linalg?模塊將線性代數(shù)運(yùn)算應(yīng)用于 n 維數(shù)組。有關(guān)這方面的更多信息，請查看?scipy.linalg 參考。

要繼續(xù)，請從 NumPy 導(dǎo)入線性代數(shù)子模塊：

>>> from numpy import linalg

為了從給定的矩陣中提取信息，我們可以使用 SVD 獲得 3 個數(shù)組，這些數(shù)組可以相乘以獲得原始矩陣。從線性代數(shù)理論，給定一個矩陣A，可以計算以下乘積： UΣVT=A

在哪里 U?和?VT 是正方形和?A.?Σ 是對角矩陣，并且包含?的奇異值的 A 從大到小排列。這些值總是非負(fù)的，可以用作矩陣表示的某些特征的“重要性”的指標(biāo) A. 讓我們先看看這在實(shí)踐中是如何使用一個矩陣的。請注意，根據(jù)色度學(xué)，如果我們應(yīng)用公式，則可以獲得我們彩色圖像的相當(dāng)合理的灰度版本。 Y=0.2126R+0.7152G+0.0722B 在哪里?Y 是表示灰度圖像的數(shù)組，和?R,G 和?B 是我們最初擁有的紅色、綠色和藍(lán)色通道陣列。請注意，我們可以@為此使用運(yùn)算符（NumPy 數(shù)組的矩陣乘法運(yùn)算符，請參閱 參考資料numpy.matmul）：

>>> img_gray = img_array @ [0.2126, 0.7152, 0.0722]

現(xiàn)在，img_gray有形狀

>>> img_gray.shape
(768, 1024)

要查看這在我們的圖像中是否有意義，我們應(yīng)該使用?matplotlib與我們希望在輸出圖像中看到的顏色相對應(yīng)的顏色圖（否則，matplotlib將默認(rèn)為與實(shí)際數(shù)據(jù)不對應(yīng)的顏色圖）。

在我們的例子中，我們正在逼近圖像的灰度部分，因此我們將使用顏色圖gray：

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

現(xiàn)在，將該linalg.svd函數(shù)應(yīng)用于該矩陣，我們得到以下分解：

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

筆記
如果您使用自己的映像，則此命令可能需要一段時間才能運(yùn)行，具體取決于您的映像大小和硬件。別擔(dān)心，這是正常的！SVD 可能是一個非常密集的計算。如果您使用自己的映像，則此命令可能需要一段時間才能運(yùn)行，具體取決于您的映像大小和硬件。別擔(dān)心，這是正常的！SVD 可能是一個非常密集的計算。

讓我們檢查一下這是否是我們所期望的：

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

請注意，s它具有特定的形狀：它只有一個維度。這意味著某些需要二維數(shù)組的線性代數(shù)函數(shù)可能不起作用。例如，從理論來看，人們可能期望乘法s并Vt與之相容。然而，這是不正確的，因?yàn)?code>s它沒有第二個軸。執(zhí)行

>>> 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)

結(jié)果是ValueError.?發(fā)生這種情況是因?yàn)?code>s在這種情況下，對于 的一維數(shù)組在實(shí)踐中比使用相同的數(shù)據(jù)構(gòu)建對角矩陣要經(jīng)濟(jì)得多。為了重建原始矩陣，我們可以重建對角矩陣Σs其對角線上的元素和適當(dāng)?shù)某朔ňS度：在我們的例子中，Σ應(yīng)該U是 768x1024，因?yàn)槭?768x768 和 1024x1024?Vt。

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

現(xiàn)在，我們要檢查重建的矩陣是否接近原始矩陣。U?@?Sigma?@?Vt``img_gray

近似值

該linalg模塊包括一個norm函數(shù)，該函數(shù)計算以 NumPy 數(shù)組表示的向量或矩陣的范數(shù)。例如，從上面的 SVD 解釋中，我們期望img_gray重構(gòu)的 SVD 乘積之間的差異范數(shù)很小。正如預(yù)期的那樣，你應(yīng)該看到類似的東西

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

（此操作的實(shí)際結(jié)果可能會因您的體系結(jié)構(gòu)和線性代數(shù)設(shè)置而異。無論如何，您應(yīng)該看到一個小數(shù)字。） 我們也可以使用該numpy.allclose函數(shù)來確保重構(gòu)的乘積實(shí)際上接近我們的原始矩陣（兩個數(shù)組之間的差異很小）：

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

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

>>> plt.plot(s)

在圖中，我們可以看到，雖然我們在 中有 768 個奇異值?s，但大多數(shù)（在第 150 個條目之后）都非常小。因此，僅使用與第一個（例如 50 個）奇異值相關(guān)的信息來構(gòu)建更經(jīng)濟(jì)的圖像近似值可能是有意義的?。

這個想法是將除了第一個k奇異值之外的?所有Sigma（與 相同s）視為零，保持?U和Vt完整，并計算這些矩陣的乘積作為近似值。 例如，如果我們選擇

>>> k = 10

我們可以通過做來建立近似

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

請注意，我們必須僅使用 的第一k行Vt，因?yàn)樗衅渌卸紝⒊艘詫?yīng)于我們從該近似中消除的奇異值的零。

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

現(xiàn)在，您可以繼續(xù)使用k 的其他值重復(fù)此實(shí)驗(yàn)，并且您的每個實(shí)驗(yàn)都應(yīng)該根據(jù)您選擇的值為您提供稍好（或更差）的圖像。

適用于所有顏色

現(xiàn)在我們想要做同樣的操作，但是對所有三種顏色。我們的第一直覺可能是對每個顏色矩陣分別重復(fù)我們上面所做的相同操作。然而，NumPy 的廣播為我們處理了這個問題。

如果我們的數(shù)組有兩個以上的維度，那么 SVD 可以一次應(yīng)用于所有軸。但是，NumPy 中的線性代數(shù)函數(shù)期望看到形式為 的數(shù)組，其中第一個軸表示矩陣的數(shù)量。(N,?M,?M)

在我們的例子中，

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

所以我們需要排列這個數(shù)組上的軸以獲得像?.?幸運(yùn)的是，該函數(shù)可以為我們做到這一點(diǎn)：(3,?768,?1024)``numpy.transpose

np.transpose(x, axes=(i, j, k))

表示將重新排序軸，以便轉(zhuǎn)置數(shù)組的最終形狀將根據(jù)索引重新排序。(i,?j,?k) 讓我們看看我們的數(shù)組如何：

>>> img_array_transposed = np.transpose(img_array, (2, 0, 1))
>>> img_array_transposed.shape
(3, 768, 1024)

現(xiàn)在我們準(zhǔn)備好應(yīng)用 SVD：

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

最后，為了獲得完整的近似圖像，我們需要將這些矩陣重新組合成近似值。現(xiàn)在，請注意

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

要構(gòu)建最終的近似矩陣，我們必須了解跨不同軸的乘法是如何工作的。

具有 n 維數(shù)組的產(chǎn)品

如果您之前只在 NumPy 中使用過一維或二維數(shù)組，則可以交替使用numpy.dotand?numpy.matmul（或@運(yùn)算符）。但是，對于 n 維數(shù)組，它們的工作方式非常不同。有關(guān)更多詳細(xì)信息，請查看文檔numpy.matmul。

現(xiàn)在，為了建立我們的近似值，我們首先需要確保我們的奇異值已準(zhǔn)備好進(jìn)行乘法運(yùn)算，因此我們Sigma以與之前所做的類似的方式建立我們的矩陣。所述Sigma陣列必須具有尺寸?。為了將奇異值添加到 的對角線上?，我們將使用NumPy 中的函數(shù)，使用 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, :])

現(xiàn)在，如果我們希望重建完整的 SVD（沒有近似值），我們可以這樣做

>>> reconstructed = U @ Sigma @ Vt

注意

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

和

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

應(yīng)該為您提供與原始圖像無法區(qū)分的圖像（盡管我們可能會為此重建引入浮點(diǎn)錯誤）。事實(shí)上，您可能會看到一條警告消息說“將輸入數(shù)據(jù)裁剪到有效范圍以使用 RGB 數(shù)據(jù)進(jìn)行 imshow（[0..1] 表示浮點(diǎn)數(shù)或 [0..255] 表示整數(shù)）。”?這是我們剛剛對原始圖像進(jìn)行的操作所預(yù)期的。

現(xiàn)在，要進(jìn)行近似，我們必須僅為k每個顏色通道選擇第一個奇異值。這可以使用以下語法完成：

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

您可以看到我們只選擇k了最后一個軸的第一個組件Sigma（這意味著我們只使用k了堆棧中三個矩陣中每一個的第一列），并且我們只選擇k了第二個中的第一個組件-to-last 軸Vt（這意味著我們只k從堆棧中的每個矩陣Vt和所有列中選擇了第一行）。如果您不熟悉省略號語法，則它是其他軸的占位符。有關(guān)更多詳細(xì)信息，請參閱有關(guān)Indexing的文檔?。

現(xiàn)在，

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

這不是顯示圖像的正確形狀。最后，將軸重新排序回我們的原始形狀，我們可以看到我們的近似值：(768,?1024,?3)

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

即使圖像不那么清晰，使用少量k奇異值（與原始的 768 個值集相比），我們可以從該圖像中恢復(fù)許多區(qū)別特征。

總結(jié)

當(dāng)然，這不是近似圖像的最佳方法。然而，事實(shí)上，線性代數(shù)中的一個結(jié)果表明，我們上面建立的近似值是我們在差分范數(shù)方面可以得到的最好的原始矩陣。有關(guān)更多信息，請參閱GH Golub 和 CF Van Loan，矩陣計算，巴爾的摩，醫(yī)學(xué)博士，約翰霍普金斯大學(xué)出版社，1985 年。