Numpy的理解与应用（三）

（本文接上文，为numpy中的线性代数部分。）
三、 线性代数
NumPy 中包含了一个矩阵库 numpy.matlib，该模块中的函数返回的是一个矩阵，而不仅是 ndarray 对象。NumPy 提供了线性代数函数库 linalg，该库包含了线性代数所需的所有功能。本文做一些简要说明。
1、创建矩阵
Matlib中可以用简单的语句建立一些基本的矩阵，例如单位矩阵，对角矩阵，0值矩阵，1值矩阵，随机数矩阵

import numpy.matlib
print(np.matlib.eye(3,4))
print('eye矩阵对角线右移2位')
print(np.matlib.eye(3,4,2))
print('单位矩阵')
print(np.matlib.identity(5,dtype=float))

运行结果为：

[[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]]
eye矩阵对角线右移2位
[[0. 0. 1. 0.]
 [0. 0. 0. 1.]
 [0. 0. 0. 0.]]
单位矩阵
[[1. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0.]
 [0. 0. 1. 0. 0.]
 [0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 1.]]

2、矩阵的基本操作：矩阵的转置、迹、行列式、特征值等
1）转置：行列互换

import numpy as np
a=np.array([[1,2],[3,4]])
print('数组 a')
print(a)
print('a的转置')
print(a.transpose())

运行结果是：

数组 a
[[1 2]
 [3 4]]
a的转置
[[1 3]
 [2 4]]

2）、逆：设A是数域上的一个n阶矩阵，若在相同数域上存在另一个n阶矩阵B，使得： AB=BA=E ，则我们称B是A的逆矩阵，而A则被称为可逆矩阵。注：E为单位矩阵。

import numpy as np
import numpy.linalg
a=np.array([[1,2],[3,4]])
print('数组 a')
print(a)
print('a的逆矩阵')
print(np.linalg.inv(a))

运行结果是：

数组 a
[[1 2]
 [3 4]]
a的逆矩阵
[[-2.   1. ]
 [ 1.5 -0.5]]

3）、迹：在线性代数中，一个n×n矩阵A的主对角线（从左上方至右下方的对角线）上各个元素的总和被称为矩阵A的迹（或迹数），一般记作tr(A)。

import numpy as np
a=np.array([[1,2],[3,4]])
print('数组 a')
print(a)
print('a的迹')
print(np.trace(a))

运行结果是：

数组 a
[[1 2]
 [3 4]]
a的迹
5

4）、行列式：行列式在是一个函数，取值为一个标量，写作det(A)或 | A | 。行列式可以看做是有向面积或体积的概念在一般的欧几里得空间中的推广。或者说，在 n 维欧几里得空间中，行列式描述的是一个线性变换对“体积”所造成的影响。它从方阵的对角元素计算。 对于 2×2 矩阵，它是左上和右下元素的乘积与其他两个的乘积的差。换句话说，对于矩阵[[a，b]，[c，d]]，行列式计算为 ad-bc。 较大的方阵被认为是 2×2 矩阵的组合。

import numpy as np
import numpy.linalg
b=np.array([[6,1,1],[4,-2,5],[2,8,7]])
print(b)
print(np.linalg.det(b))
print(6*(-2*7-5*8)-1*(4*7-5*2)+1*(4*8-(-2)*2))

运行结果是：

[[ 6  1  1]
 [ 4 -2  5]
 [ 2  8  7]]
-306.0
-306

5）、奇异值分解：假设A是一个M* N的矩阵，那么通过矩阵分解将会得到U，Σ，V^T（V的转置）三个矩阵。

其中U是一个M * M的方阵，被称为左奇异向量，方阵里面的向量是正交的；Σ是一个M* N的对角矩阵，除了对角线的元素其他都是0，对角线上的值称为奇异值；V^T(V的转置)是一个N * N的矩阵，被称为右奇异向量，方阵里面的向量也都是正交的。
numpy.linalg模块中的svd函数可以对矩阵进行奇异值分解。该函数返回3个矩阵——U、Sigma和V^T，其中U和V的转置是正交矩阵， Sigma包含输入矩阵的奇异值。

import numpy as np
import numpy.matlib
import numpy.linalg
D=np.mat("4 11 14;8 7 -2")
U,Sigma,V=np.linalg.svd(D,full_matrices=False)
print('原数组')
print(D)
print("U:",U)
print("Sigma:",Sigma)
print("V:",V)
print("验证")
print(U*np.diag(Sigma)*V)

运行结果是：

原数组
[[ 4 11 14]
 [ 8  7 -2]]
U: [[-0.9486833  -0.31622777]
 [-0.31622777  0.9486833 ]]
Sigma: [18.97366596  9.48683298]
V: [[-0.33333333 -0.66666667 -0.66666667]
 [ 0.66666667  0.33333333 -0.66666667]]
验证
[[ 4. 11. 14.]
 [ 8.  7. -2.]]

6）、特征值：设A是n阶方阵，如果数λ和n维非零列向量x使关系式Ax=λx成立，那么这样的数λ称为矩阵A特征值，非零向量x称为A的对应于特征值λ的特征向量。式Ax=λx也可写成( A-λE)X=0。这是n个未知数n个方程的齐次线性方程组，它有非零解的充分必要条件是系数行列式| A-λE|=0。
可以调用eigvals函数求解特征值，使用eig函数求解特征值和特征向量 (该函数将返回一个元组，按列排放着特征值和对应的特征向量，其中第一列为特征值，第二列为特征向量)

import numpy as np
import numpy.matlib
import numpy.linalg
c=np.mat("3 -2;1 0")
print(c)
print('使用eigvals求解特征值')
c0=np.linalg.eigvals(c)
print(c0)
print('使用eig求解特征值和特征向量')
c1,c2=np.linalg.eig(c)
print(c1)
print(c2)

运行结果是：

[[ 3 -2]
 [ 1  0]]
使用eigvals求解特征值
[2. 1.]
使用eig求解特征值和特征向量
[2. 1.]
[[0.89442719 0.70710678]
 [0.4472136  0.70710678]]

3、矩阵运算。在numpy中，面向向量和矩阵的函数常常是不一样的，使用时需要留意。结合线性代数的基本运算，可以从以下维度理解：
1）、点积。返回一个数值，两个矩阵使用vdot，多个矩阵使用linalg.multi_dot。两个矩阵对应位置的元素乘积之和，维数必须匹配。此处，向量可以用同样的函数操作。

import numpy as np
import numpy.matlib
import numpy.linalg
print('点积')
a=np.array([[1,2],[3,4]])
b=np.array([[5,6],[7,8]])
print(a)
print(b)
print(np.vdot(a,b))

运行结果是：

点积
[[1 2]
 [3 4]]
[[5 6]
 [7 8]]
70

2）、矩阵内积。（2个矩阵的内积使用inner函数）

2个向量的内积也使用inner函数

import numpy as np
import numpy.matlib
import numpy.linalg
a=np.array([[1,2],[3,4]])
b=np.array([[5,6],[7,8]])
print(a)
print(b)
print('矩阵内积')
print(np.inner(a,b))

运行结果是：

[[1 2]
 [3 4]]
[[5 6]
 [7 8]]
矩阵内积
[[17 23]
 [39 53]]

3）、向量外积。（2个向量的外积使用outer函数）
向量外积的定义为：
|a ×b| = |a|•|b|•sin<a,b>.
方向根据右手法则确定，就是手掌立在a、b所在平面的向量a上，掌心由a转向b的过程中，大拇指的方向就是外积的方向。

import numpy as np
import numpy.matlib
import numpy.linalg
a=np.array([1,2,3])
b=np.array([4,5,6])
print("a:",a)
print("b:",b)
print("向量外积")
print(np.outer(a,b))

运行结果是：

a: [1 2 3]
b: [4 5 6]
向量外积
[[ 4  5  6]
 [ 8 10 12]
 [12 15 18]]

4）、矩阵相乘。（2个矩阵相乘使用matmul）
matmul计算的是线性代数中的矩阵相乘，即NM的矩阵必须与MK的矩阵（N、M、K可以为1。）。同时这个函数具备广播功能，并且可以处理高维的数组。

import numpy as np
import numpy.matlib
import numpy.linalg
a=np.array([[1,2],[3,4]])
b=np.array([[5,6],[7,8]])
print('a:',a)
print('b:',b)
print('矩阵a,b相乘')
print(np.matmul(a,b))

运行结果是：

a: [[1 2]
 [3 4]]
b: [[5 6]
 [7 8]]
矩阵a,b相乘
[[19 22]
 [43 50]]

*高维矩阵的乘法，会自动广播，并照矩阵乘法的顺序操作

import numpy as np
import numpy.matlib
import numpy.linalg
c=np.arange(1,9).reshape(2,2,2)
a=np.array([[1,2],[3,4]])
print(np.matmul(c,a))

运行结果是：

[[[ 7 10]
  [15 22]]

 [[23 34]
  [31 46]]]

5）、线性方程组求解：numpy.linalg中的函数solve可以求解形如 Ax = b 的线性方程组，其中 A 为矩阵，b 为一维或二维的数组，x 是未知变量。

import numpy as np
import numpy.matlib
import numpy.linalg
B=np.mat("1 -2 1;0 2 -8;-4 5 9")
b=np.array([0,8,-9])
print("系数矩阵：",B)
print("目标向量:",b)
x=np.linalg.solve(B,b)
print("x:",x)

运行结果是：

系数矩阵： [[ 1 -2  1]
 [ 0  2 -8]
 [-4  5  9]]
目标向量: [ 0  8 -9]
x: [29. 16.  3.]

6）、Numpy还可以执行张量积、克罗内克积、稀疏矩阵等运算。具体可参考：https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.linalg.html