线性代数模块¶

scipy.linalg是SciPy中负责线性代数计算的模块。NumPy也有numpy.linalg模块，不过建议使用scipy.linalg，原因有二：

scipy.linalg包含numpy.linalg中的所有函数，同时还包含了很多numpy.linalg中没有的函数；
scipy.linalg默认使用线性代数库BLAS/LAPACK等对运算进行加速，而在numpy.linalg中，这些加速是需要自己配置，不是默认设置:

In [1]:

import numpy.linalg

In [2]:

import scipy.linalg

In [3]:

len(dir(numpy.linalg))

Out[3]:

In [4]:

len(dir(scipy.linalg))

Out[4]:

In [5]:

import numpy as np

In [6]:

from scipy import linalg

基本的矩阵操作¶

在NumPy中，矩阵有两种表示方法：矩阵类型和二维数组类型。

矩阵类型可以用np.mat()或者np.matrix()创建：

In [7]:

A = np.mat("[1, 2; 3, 4]")

In [8]:

Out[8]:

matrix([[1, 2],
        [3, 4]])

In [9]:

A = np.matrix("[1, 2; 3, 4]")

In [10]:

Out[10]:

matrix([[1, 2],
        [3, 4]])

转置矩阵为：

In [11]:

A.T

Out[11]:

matrix([[1, 3],
        [2, 4]])

矩阵的逆矩阵：

In [12]:

A.I

Out[12]:

matrix([[-2. ,  1. ],
        [ 1.5, -0.5]])

矩阵乘法：

In [13]:

A.I * A

Out[13]:

matrix([[1.00000000e+00, 0.00000000e+00],
        [1.11022302e-16, 1.00000000e+00]])

In [14]:

A * A.I

Out[14]:

matrix([[1.0000000e+00, 0.0000000e+00],
        [8.8817842e-16, 1.0000000e+00]])

数组类型也可以完成类似的操作：

In [15]:

A = np.array([[1, 2], [3, 4]])

In [16]:

Out[16]:

array([[1, 2],
       [3, 4]])

转置：

In [17]:

A.T

Out[17]:

array([[1, 3],
       [2, 4]])

其逆矩阵可以用linalg.inv()函数计算：

In [18]:

A.dot(linalg.inv(A))

Out[18]:

array([[1.0000000e+00, 0.0000000e+00],
       [8.8817842e-16, 1.0000000e+00]])

In [19]:

A @ linalg.inv(A)

Out[19]:

array([[1.0000000e+00, 0.0000000e+00],
       [8.8817842e-16, 1.0000000e+00]])

In [20]:

linalg.inv(A) @ A

Out[20]:

array([[1.00000000e+00, 0.00000000e+00],
       [1.11022302e-16, 1.00000000e+00]])

scipy.linalg也可以作用在矩阵类型上，返回的还是一个数组：

In [21]:

A = np.matrix([[1, 2], [3, 4]])

In [22]:

linalg.inv(A)

Out[22]:

array([[-2. ,  1. ],
       [ 1.5, -0.5]])

矩阵类型的使用比较少，通常都用数组类型替代矩阵类型进行计算。

线性方程组的求解¶

考虑下列方程组：

 x + 3y + 5z = 10
2x + 5y -  z = 8
2x + 3y + 8z = 3

矩阵形式为：

$\left[ \begin{array}{1} 1 & 3 & 5 \\ 2 & 5 & -1 \\ 2 & 3 & 8 \\ \end{array} \right] \left[ \begin{array}{1} x \\ y \\ z \\ \end{array} \right] = \left[ \begin{array}{1} 10 \\ 8 \\ 3 \\ \end{array} \right]$

其解为：

$\left[ \begin{array}{1} x \\ y \\ z \\ \end{array} \right] = \left[ \begin{array}{1} 1 & 3 & 5 \\ 2 & 5 & -1 \\ 2 & 3 & 8 \\ \end{array} \right]^{-1} \left[ \begin{array}{1} 10 \\ 8 \\ 3 \\ \end{array} \right]$

故：

In [23]:

A = np.array([[1, 3, 5],
              [2, 5, -1],
              [2, 3, 8]])

In [24]:

b = np.array([10, 8, 3])

解：

In [25]:

linalg.inv(A) @ b

Out[25]:

array([-8.83870968,  5.25806452,  0.61290323])

也可以直接使用linalg.solve()函数，直接求解这个方程：

In [26]:

linalg.solve(A, b)

Out[26]:

array([-8.83870968,  5.25806452,  0.61290323])

行列式的计算：¶

linalg.det()函数可以计算矩阵的行列式：

In [27]:

linalg.det(A)

Out[27]:

-31.000000000000007

行列式非零表示矩阵是可逆的。逆矩阵的行列式与矩阵的行列式满足相乘等于1的关系：

In [28]:

linalg.det(linalg.inv(A))

Out[28]:

-0.03225806451612903

In [29]:

linalg.det(linalg.inv(A)) * linalg.det(A)

Out[29]:

1.0000000000000002

范数的计算¶

In [30]:

A = np.array([[1, 2], [3, 4]])

默认情况下，linalg.norm()函数计算的是矩阵的Frobenius范数，该范数计算的是矩阵所有元素平方和的平方根：

In [31]:

linalg.norm(A)

Out[31]:

5.477225575051661

In [32]:

linalg.norm(A, "fro")

Out[32]:

5.477225575051661

可以通过linalg.norm()函数中的第二个参数，来修改矩阵范数的计算方法。例如，矩阵的1范数，计算的是矩阵中每列元素的和的最大值：

In [33]:

linalg.norm(A, 1)

Out[33]:

6.0

矩阵的-1范数，计算的是矩阵中每列元素的和的最小值：

In [34]:

linalg.norm(A, -1)

Out[34]:

4.0

矩阵的2范数，计算的是矩阵的最大奇异值：

In [35]:

linalg.norm(A, 2)

Out[35]:

5.464985704219043

矩阵的-2范数，计算的是矩阵的最小奇异值：

In [36]:

linalg.norm(A, -2)

Out[36]:

0.36596619062625746

矩阵的正无穷范数，计算的是矩阵中每行元素的和的最大值：

In [37]:

linalg.norm(A, np.inf)

Out[37]:

7.0

矩阵的负无穷范数，计算的是矩阵中每行元素的和的最小值：

In [38]:

linalg.norm(A, -np.inf)

Out[38]:

3.0

linalg.norm()函数也可以用来计算向量的范数（也称模）：

In [39]:

b = np.array([10, 8, 3])

对于向量来说，linalg.norm()函数默认是向量的2范数，即计算向量元素的平方和的平方根：

In [40]:

linalg.norm(b)

Out[40]:

13.152946437965905

In [41]:

linalg.norm(b, 2)

Out[41]:

13.152946437965905

特征值分解¶

对于方阵A，特征值λ和特征向量v满足：Av=λv。函数linalg.eig()可以用来求解特征值和特征向量：

In [42]:

A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9]])

In [43]:

l, v = linalg.eig(A)

In [44]:

np.allclose(A.dot(v), l * v)

Out[44]:

True

广义逆¶

对于形状为m×n的矩阵A，可以用linalg.pinv()求解其广义逆。A的广义逆B满足：A=ABA，B=ABA：

In [45]:

A = np.array([[1, 2, 3],
              [4, 5, 6]])

In [46]:

B = linalg.pinv(A)

In [47]:

np.allclose(A, A @ B @ A)

Out[47]:

True

In [48]:

np.allclose(B, B @ A @ B)

Out[48]:

True

奇异值分解¶

对于形状为m×n的矩阵A，其奇异值分解满足 $A=U\SigmaV^H$ ，矩阵Σ形状为M×N，主对角线上的元素被称为奇异值。U、V分别为M阶、N阶正交单位矩阵。

函数linalg.svd()可以对矩阵进行奇异值分解：

In [49]:

U, s, Vh = linalg.svd(A)

In [50]:

Out[50]:

array([[-0.3863177 , -0.92236578],
       [-0.92236578,  0.3863177 ]])

In [51]:

Out[51]:

array([9.508032  , 0.77286964])

In [52]:

Vh

Out[52]:

array([[-0.42866713, -0.56630692, -0.7039467 ],
       [ 0.80596391,  0.11238241, -0.58119908],
       [ 0.40824829, -0.81649658,  0.40824829]])

从奇异值恢复矩阵 $\Sigma$ ：

In [53]:

S = linalg.diagsvd(s, 2, 3)

In [54]:

np.allclose(A, U @ S @ Vh)

Out[54]:

True

正交矩阵：

In [55]:

U @ U.T

Out[55]:

array([[ 1.00000000e+00, -8.49433239e-18],
       [-8.49433239e-18,  1.00000000e+00]])

In [56]:

Vh @ Vh.T

Out[56]:

array([[ 1.00000000e+00,  4.13108370e-17,  2.21870637e-17],
       [ 4.13108370e-17,  1.00000000e+00, -1.77008767e-16],
       [ 2.21870637e-17, -1.77008767e-16,  1.00000000e+00]])

In [ ]: