簡介
本文將會以圖表的形式為大家講解怎麼在NumPy中進行多維資料的線性代數運算。
多維資料的線性代數通常被用在影像處理的圖形變換中，本文將會使用一個影像的例子進行說明。
圖形載入和說明
熟悉顏色的朋友應該都知道，一個顏色可以用R，G，B來表示，如果更高階一點，那麼還有一個A表示透明度。通常我們用一個四個屬性的陣列來表示。
對於一個二維的影像來說，其解析度可以看做是一個X*Y的矩陣，矩陣中的每個點的顏色都可以用（R，G，B）來表示。
有了上面的知識，我們就可以對影像的顏色進行分解了。
首先需要載入一個影像，我們使用imageio.imread方法來載入一個本地影像，如下所示：
import imageio
img=imageio.imread('img.png')
print(type(img))
上面的程式碼從本地讀取圖片到img物件中，使用type可以檢視img的型別，從執行結果，我們可以看到img的型別是一個陣列。
class 'imageio.core.util.Array'
通過img.shape可以得到img是一個(80, 170, 4)的三維陣列，也就是說這個影像的解析度是80*170，每個畫素是一個（R，B，G，A）的陣列。
最後將影像畫出來如下所示：
import matplotlib.pyplot as plt
plt.imshow(img)
圖形的灰度
對於三維陣列來說，我們可以分別得到三種顏色的陣列如下所示：
red_array = img_array[:, :, 0]
green_array = img_array[:, :, 1]
blue_array = img_array[:, :, 2]
有了三個顏色之後我們可以使用下面的公式對其進行灰度變換：
Y=0.2126R + 0.7152G + 0.0722B
上圖中Y表示的是灰度。
怎麼使用矩陣的乘法呢？使用 @ 就可以了：
img_gray = img_array @ [0.2126, 0.7152, 0.0722]
現在img是一個80 * 170的矩陣。
現在使用cmap="gray"作圖：
plt.imshow(img_gray, cmap="gray")
可以得到下面的灰度影像：
灰度影像的壓縮
灰度影像是對影像的顏色進行變換，如果要對影像進行壓縮該怎麼處理呢？
矩陣運算中有一個概念叫做奇異值和特徵值。
設A為n階矩陣，若存在常數λ及n維非零向量x，使得Ax=λx，則稱λ是矩陣A的特徵值，x是A屬於特徵值λ的特徵向量。
一個矩陣的一組特徵向量是一組正交向量。
即特徵向量被施以線性變換 A 只會使向量伸長或縮短而其方向不被改變。
特徵分解（Eigendecomposition），又稱譜分解（Spectral decomposition）是將矩陣分解為由其特徵值和特徵向量表示的矩陣之積的方法。
假如A是m * n階矩陣，q=min(m,n)，A*A的q個非負特徵值的算術平方根叫作A的奇異值。
特徵值分解可以方便的提取矩陣的特徵，但是前提是這個矩陣是一個方陣。如果是非方陣的情況下，就需要用到奇異值分解了。先看下奇異值分解的定義：
\(A=UΣV^T\)
其中A是目標要分解的m * n的矩陣，U是一個 m * m的方陣，Σ 是一個m * n 的矩陣，其非對角線上的元素都是0。\(V^T\)是V的轉置，也是一個n * n的矩陣。
奇異值跟特徵值類似，在矩陣Σ中也是從大到小排列，而且奇異值的減少特別的快，在很多情況下，前10%甚至1%的奇異值的和就佔了全部的奇異值之和的99%以上了。也就是說，我們也可以用前r大的奇異值來近似描述矩陣。r是一個遠小於m、n的數，這樣就可以進行壓縮矩陣。
通過奇異值分解，我們可以通過更加少量的資料來近似替代原矩陣。
要想使用奇異值分解svd可以直接呼叫linalg.svd 如下所示：
U, s, Vt = linalg.svd(img_gray)
其中U是一個m * m矩陣，Vt是一個n * n矩陣。
在上述的影像中，U是一個(80, 80)的矩陣，而Vt是一個(170, 170) 的矩陣。而s是一個80的陣列，s包含了img中的奇異值。
如果將s用影像來表示，我們可以看到大部分的奇異值都集中在前的部分：
這也就意味著，我們可以取s中前面的部分值來進行影像的重構。
使用s對影像進行重構，需要將s還原成80 * 170 的矩陣：
# 重建
import numpy as np
Sigma = np.zeros((80, 170))
for i in range(80):
Sigma[i, i] = s[i]
使用 U @ Sigma @ Vt 即可重建原來的矩陣，可以通過計算linalg.norm來比較一下原矩陣和重建的矩陣之間的差異。
linalg.norm(img_gray - U @ Sigma @ Vt)
或者使用np.allclose來比較兩個矩陣的不同：
np.allclose(img_gray, U @ Sigma @ Vt)
或者只取s陣列的前10個元素，進行重新繪圖，比較一下和原圖的區別：
k = 10
approx = U @ Sigma[:, :k] @ Vt[:k, :]
plt.imshow(approx, cmap="gray")
可以看到，差異並不是很大：
原始影像的壓縮
上一節我們講到了如何進行灰度影像的壓縮，那麼如何對原始影像進行壓縮呢？
同樣可以使用linalg.svd對矩陣進行分解。
但是在使用前需要進行一些處理，因為原始影像的img_array 是一個(80, 170, 3)的矩陣--這裡我們將透明度去掉了，只保留了R，B，G三個屬性。
在進行轉換之前，我們需要把不需要變換的軸放到最前面，也就是說將index=2，換到index=0的位置，然後進行svd操作：
img_array_transposed = np.transpose(img_array, (2, 0, 1))
print(img_array_transposed.shape)
U, s, Vt = linalg.svd(img_array_transposed)
print(U.shape, s.shape, Vt.shape)
同樣的，現在s是一個(3, 80)的矩陣，還是少了一維，如果重建影像，需要將其進行填充和處理，最後將重建的影像輸出：
Sigma = np.zeros((3, 80, 170))
for j in range(3):
np.fill_diagonal(Sigma[j, :, :], s[j, :])
reconstructed = U @ Sigma @ Vt
print(reconstructed.shape)
plt.imshow(np.transpose(reconstructed, (1, 2, 0)))
當然，也可以選擇前面的K個特徵值對影像進行壓縮：
approx_img = U @ Sigma[..., :k] @ Vt[..., :k, :]
print(approx_img.shape)
plt.imshow(np.transpose(approx_img, (1, 2, 0)))
重新構建的影像如下：
對比可以發現，雖然損失了部分精度，但是影像還是可以分辨的。
總結
影像的變化會涉及到很多線性運算，大家可以以此文為例，仔細研究。
本文已收錄於 http://www.flydean.com/08-python-numpy-linear-algebra/
最通俗的解讀，最深刻的乾貨，最簡潔的教程，眾多你不知道的小技巧等你來發現！
歡迎關注我的公眾號:「程式那些事」,懂技術，更懂你！