Numpy和矩阵运算

March 2019 3 minute read

最近在复习线性代数,免不了要用到Numpy,所以在这里顺便总结一下Numpy里面常用的API.

基础

一个Numpy的ndarray对象其实就是一个多维的数组,这个数组里面的元素都是同一个类型的。

在这里我列举一些ndarray里面比较重要的属性。

生成ndarray

a = np.array([[[np.random.randint(1,10) for _ in range(3)] for _ in range(4)]for _ in range(5)])
>>c = np.array( [ [1,2], [3,4] ], dtype=complex )
>>c
array([[ 1.+0.j,  2.+0.j],
       [ 3.+0.j,  4.+0.j]])
>>np.ones((1,2),dtype=np.float64)
array([[1., 1.]])
>>np.zeros((4,3),dtype=np.float64)
array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])
>>np.arange(1,10,1)
array([1, 2, 3, 4, 5, 6, 7, 8, 9])
np.linspace(0.1,1,10)

这个技巧和可以和matplotlib 结合,可以很方便的画一些图

import matplotlib.pyplot as plt
x = np.linspace( 0, 2*pi, 100)
plt.plot(np.sin(x))

plot_sin.png

改变形状

>>a = np.floor(10*np.random.random((3,4)))
>>a
array([[5., 2., 5., 1.],
       [1., 3., 1., 1.],
       [0., 6., 5., 1.]])
>>a.ravel()
array([5., 2., 5., 1., 1., 3., 1., 1., 0., 6., 5., 1.])
>>a = np.ceil(np.random.random((3,4))*10)
>>a
array([[10.,  8.,  8.,  7.],
       [ 8.,  8.,  1., 10.],
       [ 3.,  7., 10.,  3.]])
>>a.reshape((6,2))
array([[10.,  8.],
       [ 8.,  7.],
       [ 8.,  8.],
       [ 1., 10.],
       [ 3.,  7.],
       [10.,  3.]])

最后一个维度可以省略,用-1 代替即可。

>>b = a.reshape(2,3,-1)
>>b
array([[[10.,  8.],
        [ 8.,  7.],
        [ 8.,  8.]],

       [[ 1., 10.],
        [ 3.,  7.],
        [10.,  3.]]])
>>b.shape
(2, 2, 3)

合并

>>a = np.ceil(10*np.random.random((2,2)))
>>a
array([[7., 7.],
       [4., 2.]])
>>b = np.ceil(10*np.random.random((2,2)))
>>b
array([[ 5.,  1.],
       [10.,  5.]])
>>np.vstack(a,b)
array([[ 7.,  7.],
       [ 4.,  2.],
       [ 5.,  1.],
       [10.,  5.]])
>>np.hstack((a,b))
array([[ 7.,  7.,  5.,  1.],
       [ 4.,  2., 10.,  5.]])

线性代数相关

矩阵的秩

$$ 在线性代数中,一个矩阵 {\displaystyle A}的列秩是 {\displaystyle A}的线性无关的纵列的极大数目。 $$ $$ 类似地,行秩是{\displaystyle A}的线性无关的横行的极大数目。 $$ $$ 矩阵的列秩和行秩总是相等的,因此它们可以简单地称作矩阵 {\displaystyle A}的秩。 $$ $$ 通常表示为 {\displaystyle \mathrm {r} (A)},{\displaystyle \mathrm {rank} (A)} 或 { \mathrm {rk} (A)}。 $$

在Numpy里面可以使用np.linalg.matrix_rank来计算矩阵的秩

>>np.linalg.matrix_rank(np.eye(3))
3
>>matrix_rank(np.ones((4,))) # 1 dimension - rank 1 unless all 0
1
>>matrix_rank(np.zeros((4,)))
0

对称矩阵

在线性代数中,对称矩阵(英语:symmetric matrix)是一个方形矩阵,其转置矩阵和自身相等。

在Numpy里的可以很方便地实现对称矩阵

def random_symmetric_matrix(size=4):
  matrix = np.ceil(10*np.random.random((size,size)))
  # only save the upper corner of matrix
  matrix = np.triu(matrix)
  # fold the down corner of matrix
  matrix += matrix.T - np.diag(matrix.diagonal())
  return matrix
>>a = random_symmetric_matrix(5)
>>a
array([[8., 5., 4., 9., 1.],
       [5., 2., 1., 2., 2.],
       [4., 1., 2., 8., 4.],
       [9., 2., 8., 2., 1.],
       [1., 2., 4., 1., 7.]])
>>np.all(a.T == a)
True

特征值和特征向量

$$ A\alpha=\lambda\alpha $$ $$ 左边是用矩阵A将向量\alpha 做了一个转换,右边是将向量\alpha 拉伸了\lambda倍 $$ $$ 说明A有这样一个功能:即对向量\alpha 变换后,长度拉伸\lambda倍,方向不变 $$ $$ \alpha就矩阵A的特征向量,\lambda就是矩阵A的特征值. $$

在Numpy里面,可以使用np.linalg.eig来计算特征值和特征向量

>>a = random_symmetric_matrix(3)
>>a 
array([[ 9.,  2.,  9.],
       [ 2.,  9.,  8.],
       [ 9.,  8., 10.]])
>>_lambda, alpha = np.linalg.eig(a)
>>_lambda
array([22.54663251,  7.01502574, -1.56165824])k
>>alpha
array([[-0.53197232, -0.66027727, -0.53013146],
       [-0.48743562,  0.7507226 , -0.44589471],
       [-0.69239582, -0.02120131,  0.72120631]])

如果矩阵是一个对称矩阵,则它的特征值为实数

奇异值分解

$$ 假设M是一个m×n阶矩阵,其中的元素全部属于域K,也就是实数域或复数域 $$ $$ 如此则存在一个分解使得 $$

$$ M=U\Sigma V^{*} $$

在Numpy里,可以使用numpy.linalg.svd对矩阵进行奇异值分解

>>a=random_symmetric_matrix(4)
>>U, s, V = np.linalg.svd(a, full_matrices=False)
>>np.allclose(a,U*np.dot())

用奇异值分解来实现图片压缩

import matplotlib.pyplot as plt
from sklearn.datasets import load_sample_images()

def get_image_data():
    data = datasets.load_sample_images()
    image_without_color = data.images[0][:,:,:1].reshape((data.images[0].shape[0],-1))
    return image_without_color

image_matrix = get_image_data()
plt.imshow(image_without_color, cmap=plt.cm.gray)

plot_nocolor.png

def compress_svb(image_matrix,k):
    U, s, V = np.linalg.svd(image_matrix, full_matrices=False)
    S = np.diag(s)
    np.allclose(a, np.dot(U, np.dot(S, V)))
    Uk = U[:, 0:k]
    Sk = np.diag(s[0:k])
    Vk = V[0:k, :]
    imgMat_new = np.dot(Uk, np.dot(Sk, Vk))
    return imgMat_new
plt.imshow(compress_svb(image_without_color,30), cmap=plt.cm.gray)

compress.png