inMultidimensionalArrays本文将以图表的形式讲解如何在NumPy中对多维数据进行线性代数运算。多维数据的线性代数通常用于图像处理中的图像变换。本文将使用图像示例进行说明。图形加载及说明熟悉颜色的朋友应该知道，一种颜色可以用R、G、B来表示，再高级一点的话，还有一个表示透明度的A。通常我们使用一个包含四个属性的数组来表示。对于二维图像，其分辨率可以看作是一个X*Y矩阵，矩阵中每个点的颜色可以用(R,G,B)表示。有了以上知识，我们就可以对图像进行颜色分解了。首先，需要加载图像。我们使用imageio.imread方法加载本地图片，如下：查看img的类型。从运行结果可以看出img的类型是数组。类'imageio.core.util.Array'可以通过img.shape获得。img是一个(80,170,4)的三维数组，也就是说这张图片的分辨率是80*170，每个像素点是一个(R,B,G,A)数组。最后绘制图像如下：importmatplotlib.pyplotaspltplt.imshow(img)图形的灰度对于三维数组，我们可以得到三个颜色数组如下：red_array=img_array[:,:,0]green_array=img_array[:,:,1]blue_array=img_array[:,:,2]有了三种颜色后，我们可以使用下面的公式将它们转化为灰色：Y=0.2126R+0.7152G+0.0722BY在上图代表灰度。如何使用矩阵乘法？使用@可以解决问题：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只会延长或缩短向量，而不会改变其方向。特征分解，也称为谱分解，是一种将矩阵分解成由其特征值和特征向量表示的矩阵乘积的方法。若A为mn阶矩阵，q=min(m,n)，AA的q个非负特征值的算术平方根称为A的奇异值。特征值分解可以很容易地提取特征矩阵，但前提是矩阵是方阵。如果是非方阵，则需要使用奇异值分解。先看奇异值分解的定义：$A=UΣV^T$其中A为目标要分解的mn矩阵，U为mm方阵，Σ为mn矩阵，其非对角线元素均为0.$V^T$是V的转置，也是nn的矩阵。奇异值类似于特征值。它们在矩阵Σ中也是从大到小排列的，奇异值的减少特别快。很多时候，前10%甚至1%的奇异值之和就占了所有的奇异值。总和的99%以上。也就是说，我们也可以用第r个的奇异值来近似描述矩阵。r是一个比m和n小得多的数，因此矩阵可以被压缩。通过奇异值分解，我们可以用更小的数据量逼近原始矩阵。如果要使用奇异值分解svd，可以直接调用linalg.svd如下：U,s,Vt=linalg.svd(img_gray)其中U是一个m矩阵，Vt是一个nn矩阵。上图中，U为(80,80)矩阵，Vt为(170,170)矩阵。而s是一个80的数组，s包含了img中的奇异值。如果用图像表示s，我们可以看到大部分奇异值都集中在前面部分：这意味着我们可以拿s的前面部分来重构图像。要使用s重建图像，您需要将s恢复为80*170的矩阵：#rebuildimportnumpyasnpSigma=np.zeros((80,170))foriinrange(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=10approx=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))forjinrange(3):np.fill_diagonal(Sigma[j,:,:],s[j,:])reconstructed=U@Sigma@Vtprint(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/最流行的解读，最深刻的干货，最简洁的教程，很多你不知道的小技巧等你来发现！欢迎关注我的公众号：《程序那些事儿》，懂技术，更懂你！
