Python数值计算（4）

本次介绍一下NumPy中的向量和矩阵的创建与操作，以及线性方程组的求解，分两个章节篇幅介绍。

在numpy中，ndarray表示向量和矩阵，是N-dimensional array（N维数组）的缩写，由多个具有相同类型和尺寸的元素组成的多维容器。

1. 向量和矩阵的创建

使用numpy.array创建一个numpy.ndarray的实例对象，可以创建向量或者矩阵，例如：

v1=np.array([1,2,-1])
print(type(v1),v1) # <class 'numpy.ndarray'> [ 1  2 -1]
A = np.array([[1,1,1],[2,-3,1],[1,2,-3]])
print(A)
'''
[[ 1  1  1]
[ 2 -3  1]
[ 1  2 -3]]
'''

分别创建了一个向量v1和3×3的矩阵（也可以理解为是一个二维数组）。通过该函数，还可以创建更高维度的数组。

如果单纯是创建向量，其实前面介绍到的nupmy.linspace和numpy.arange，返回的都是ndarray类型：

print(type(np.linspace(1,5,10))) # <class 'numpy.ndarray'>
print(type(np.arange(1,5,0.3))) #<class 'numpy.ndarray'>

对于一些特殊的二维矩阵，也有一些特定的函数，例如，numpy.eye/numpy.identity numpy.diag,和numpy.vander分别创建单位矩阵、对角矩阵和范德蒙德矩阵：

E=np.eye(3)
D=np.diag([1,2,3])
M=np.vander([1,2,3],3)
print(M)
'''
[[1 1 1]
[4 2 1]
[9 3 1]]
'''

其中，eye可以创建行、列不相等的单位矩阵（只有主对角线元素为1），而identity则是创建一个方阵，另外可以看到范德蒙德矩阵和我们常见的有点不大一样。

另外两个函数是numpy.zeros和numpy.ones函数，用来构建全0或者全1矩阵，可以自行测试。

2. 向量和矩阵的属性

访问向量或者矩阵中的元素，可以通过以0开始的索引实现，例如：

v=np.linspace(1,2,5)
print(v) # [1.   1.25 1.5  1.75 2.  ]
print(v[3]) # 1.75
M=np.array([[1,2,3],[2,3,4],[7,6,5]])
print(M[1,2]) # 4

对于矩阵，还可以通过类似Python中列表切片的方式，实现单独的行或者列提取：

print(M[1,:]) # 第1行 [2 3 4]
print(M[:,0]) # 第0列 [1 2 7]

还可以提取子矩阵：

M1=M[0:2,1:3]
print(M1)
'''
[[2 3]
[3 4]]
'''

注意，此处和Python中一样，切片结尾的元素是不包括在内的，M[0:2,1:3]中，是指提取0～1行，1～2列的元素，如下所示：

更多高级的操作可以参考NumPy官方文档。

多维数组的一个属性是ndim，即维度，对于标量，维度为0，向量维度为1，矩阵的维度则是2，而shape属性，可以理解为各个维度上的长度，例如：

print(np.ndim(3)) # 0
M=np.array([[1,2,3],[2,3,4],[7,6,5]])
print(M.ndim) # 2
print(M.shape) #(3, 3)
d=np.ndim(M) # d=2
m,n=np.shape(M) # m=3,n=3

上例还演示了如下事实：使用ndarray实例的属性，或者使用numpy下的同名函数，获取这两个属性值。

通过T属性，可以获取矩阵的转置：

M=np.array([[1,2,3],[2,3,4],[7,6,5]])
print(M.T)
'''
[[1 2 7]
 [2 3 6]
 [3 4 5]]
'''

3. 向量和矩阵的计算

标量和向量/矩阵的运算

标量和与向量/矩阵的运算，表现出来的和标量之间运算一致，即向量/矩阵的每一个元素分别和标量进行运算。例如：

M=np.array([[1,2,3],[2,3,4],[7,6,5]])
print(M+1)
'''
[[2 3 4]
 [3 4 5]
 [8 7 6]]
    '''   
print(M*2)
'''
[[ 2  4  6]
 [ 4  6  8]
 [14 12 10]]
 '''

向量/矩阵之间的运算

向量/矩阵进行加减运算，也是对应位置元素之间的运算，唯一需要注意的是参与运算的两个N维数组具有相同的shape，例如：

M=np.array([[1,2,3],[2,3,4],[7,6,5]])
I=np.identity(3)
E=np.eye(2,3)
print(M+I)
'''
[[2. 2. 3.]
[2. 4. 4.]
[7. 6. 6.]]
'''
print(M-E)# ValueError

最后一个会产生异常：ValueError: operands could not be broadcast together with shapes (3,3) (2,3)。

加减运算没有什么特别需要注意的地方，但是涉及到矩阵的乘除运算时，需要特别留心。

例如，我们要计算矩阵的乘法，可以使用运算符“*”吗？测试一下：

M=np.array([[1,2,3],[2,3,4],[7,6,5]])
I=np.identity(3)
print(M*I)

按矩阵计算的方式去理解，矩阵与单位阵相乘，其结果应该为原矩阵，但是这个的结果是多少呢？

[[1. 0. 0.]
 [0. 3. 0.]
 [0. 0. 5.]]

看来，也是使用了对应位置元素的乘法，而不是标准的矩阵乘法。

真正实现矩阵乘法的方式，是使用numpy.dot函数：

M=np.array([[1,2,3],[2,3,4],[7,6,5]])
I=np.identity(3)
print(np.dot(M,I))

另外一种简便的方法是使用@运算符：

C=M@I

两种情况下的结果都是M。

总而言之，就是基本运算符是元素之间的运算，而矩阵计算要使用专门的函数和运算符。

Extra：NumPy中的行向量和列向量？

在Matlab中，行向量和列向量是分的很清楚的，不过在NumPy中，则没那么明显，例如，在数学上行向量的转置必然是一个列向量，但是在numpy中并非如此：

v1=np.array([1,1,-1])
v2=np.array([2,-1,1]).T
print(v1.shape,v2.shape) # (3,) (3,)

可见两者的shape都是3，而不是1×3和3×1的（向量的ndim都只有1）。

甚至，在和矩阵进行计算时，都可以进行左乘和右乘：

v1=np.array([1,1,-1])
M=np.array([[1,2,1],[3,-1,0],[0,3,-2]])

B=v1@M
print(B) # [ 4 -2  3]
print(B.shape) # (3,)
    
C=M@v1
print(C) # [2 2 5]
print(C.shape) # (3,)

B的计算从数学上很好理解，就是一个1×3的行向量和一个3×3的矩阵相乘，结果必然是一个1×3的行向量。

但是C的计算就有点费解了，按说是一个3×3的矩阵和一个1×3的向量相乘，这个在数学上是不允许的才对，但是依旧可以计算出结果来。

在numpy中没有严格的行列向量之分，但是结果是不同的，上面两个算式的说明如下：

从表现形式上看，向量位于矩阵乘法的右边时，会理解为列向量进行计算，但是计算后的结果还是表现为行向量。

不过，可以通过矩阵的方式定义向量：

K=np.array([[1],[1],[-1]])
print(K.shape) # (3, 1)
print(M@K)
'''
[[2]
[2]
[5]]
'''
print(K@M)  # ValueError: matmul

将K定义为一个3行1列的矩阵，这样，可以计算M*K，但是不能计算K*M，错误的详细信息是：

ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 3 is different from 1)