Scipy系列(二)：线性代数及积分

Scipy系列目录

---

---

一、 概述

scipy.linalg(scipy线性代数模块)包含numpy.linalg(numpy线性代数模块)中的所有函数，并且还增加了一些高级函数。

Numpy包括

1. 矩阵和向量乘
2. 范数，行列式，秩等
3. 矩阵分解
4. 矩阵特征值和特征向量
5. 逆矩阵和解线性方程组

scipy线性代数模块可以分为几个部分：

1. 基本功能
2. 解线性方程组
3. 矩阵分解
4. 特征值和特征向量
5. 矩阵函数
6. 矩阵方程求解器

scipy积分模块可以主要包括：

1. 对给定函数进行积分
2. 对给定样本进行离散积分

二、 Scipy线性代数

Scipy大约有200个线性代数函数，很多函数与Numpy相同，因此这里只简单介绍函数功能。
Numpy中线性代数模块详细说明见Numpy系列：线性代数

1 Scipy基本功能

• scipy.linalg.issymmetric(a)：判断方阵a是否对称矩阵。
• scipy.linalg.ishermitian(a)：判断方阵a是否hermitian矩阵。
• scipy.linalg.inv(a)：求方阵a的逆矩阵。和np.linalg.inv(a)结果相同，但是多两个可选参数，较少用，一般可以认为是一样的。
• scipy.linalg.pinv(b)：求矩阵a的伪逆矩阵。同样比np.linalg.pinv(a)多几个可选参数。
• scipy.linalg.pinvh(b)：Hermitian矩阵的伪逆。相当于np.linalg.pinv(a,hermitian=True)，但是有更多可选参数。
• scipy.linalg.norm©：求矩阵或向量范数。比np.linalg.norm(a)多一个是否检查a包含无穷大的参数。
• scipy.linalg.det(a)：求方阵的行列式。比np.linalg.det(a)多两个参数，一个参数表示结果覆盖a，一个参数表示检查是否包含无穷大。
• scipy.linalg.kron(a,b)：求张量积(克罗内克积)。同np.kron(a,b)
• scipy.linalg.khatri_rao(a,b):KR积，两个列数相等的矩阵对应列向量的kron积。

KR积示例

>> A=np.array([[1,2],[3,4],[5,6]])
>> B=np.array([[2,6],[3,4]])
>> scipy.linalg.khatri_rao(A,B)
 
array([[ 2, 12],
       [ 3,  8],
       [ 6, 24],
       [ 9, 16],
       [10, 36],
       [15, 24]])

2 解线性方程组

• scipy.linalg.solve(a,b)：类似np.linalg.solve(a,b)，但是支持更多参数。求解ax=b的根，a为方阵。
• scipy.linalg.lstsq(a,b)：类似np.linalg.lstsq(a,b)，最小二乘法求解线性方程组的系数。多两个参数，是否重写覆盖参数及检查是否包含无穷大。

同时还支持几个特殊的求解线性方程的函数：

• solve_banded(l_and_u, ab, b)：求解ax=b，假设a为带状矩阵。ab为带状矩阵，l_and_u为带状矩阵非零的上、下带数。
• solveh_banded(ab, b)：求解a为hermitan带状矩阵的线性方程ax=b。
• solve_circulant(c, b)：求解cx=b方程。c为循环矩阵。
• solve_triangular(a, b)：求解ax=b方程。a为三角方阵。
• solve_toeplitz(c_or_cr, b)：求解Tx=b方程。T为托普利兹(toeplitz)矩阵。T由c_or_cr确定。

3 矩阵分解

• scipy.linalg.cholesky(a)：Cholesky分解。返回L，a=L×L.T。类似np.linalg.cholesky(a)。
• scipy.linalg.qr(a)：QR分解。返回q,r，a=q×r。类似np.linalg.qr(a)。
• scipy.linalg.svd(a)：SVD奇异值分解。返回u,s,vh，a=u×s×vh。类似np.linalg.svd(a)。
• scipy.linalg.svdvals(a)：求解a的奇异值s。
• scipy.linalg.lu(a)：LU分解，返回p,l,u，a为(M,N)形状，p为(M,M)矩阵，l为(M,K)下三角矩阵，u为(K,N)上三角或trapezoidal矩阵。

同时scipy还支持更多的矩阵分解，如ldl，polar，cossin等。

4 特征值和特征向量

• scipy.linalg.eig(a):求解方阵的特征值和特征向量，类似np.linalg.eig(a)。
• scipy.linalg.eigvals(a)：求特征值，类似np.linalg.eigvals(a)。
• scipy.linalg.eigh(a)：求Hermitian矩阵的特征值和特征向量，类似np.linalg.eigh(a)。
• scipy.linalg.eigvalsh(a)：求Hermitian矩阵的特征值，类似np.linalg.eigvalsh(a)。

scipy求解特征值和特征向量的函数支持更多的参数。同时scipy还增加了几个函数：

• scipy.linalg.eig_banded(a_band):求解实对称或复厄米带矩阵特征值和特征向量。
• scipy.linalg.eigvals_banded(a_band):求解实对称或复厄米带矩阵特征值。
• scipy.linalg.eigh_tridiagonal(d,e):求解实对称三对角矩阵的特征值和特征向量。
• scipy.linalg.eigvalsh_tridiagonal(d,e):求解实对称三对角矩阵的特征值。

5 矩阵函数

矩阵函数与numpy函数不同，比如np.sqrt(a)表示对a的每一个元素求平方根。矩阵平方根则表示求得的值x与x进行矩阵乘等于a。

矩阵A的平方表示A与A的矩阵乘积。
矩阵A的n次幂表示n个A依次进行矩阵乘。
矩阵A的平方根x表示：x与x的矩阵乘积等于A。
矩阵A作为指数，求解e的A次幂，可以用级数分解为A的n的方的函数。即可以通过A的幂求解e的A次幂。
同样矩阵A的三角函数也可以用e的iA次幂求解。

• scipy.linalg.sqrtm(A)：求矩阵A的平方根。
• scipy.linalg.expm(A)：求e的矩阵A次幂。
• scipy.linalg.logm(A)：求A对e的对数。
• scipy.linalg.sinm(A)：求A的正弦。
• scipy.linalg.cosm(A)：求A的余弦。
• scipy.linalg.tanm(A)：求A的正切。
• scipy.linalg.sinhm(A)：求A的双曲正弦。
• scipy.linalg.coshm(A)：求A的双曲余弦。
• scipy.linalg.tanhm(A)：求A的双曲正切。
• scipy.linalg.signm(A)：求方阵符号。可以参考https://nhigham.com/2020/12/15/what-is-the-matrix-sign-function/
• scipy.linalg.funm(A,fun)：对A执行fun函数。比如funm(A,lambda x:x*x)相当于x与x的矩阵乘。

6 矩阵方程求解器(控制理论方程)

scipy线型模块中提供了5个控制理论中的矩阵方程求解器。

• scipy.linalg.solve_sylverster(a,b,q)：求解控制理论中的西尔维斯特方程AX+XB=Q。a为M阶方阵，b为N阶方阵，q为(M,N)矩阵。
• scipy.linalg.solve_continuous_are(a,b,q,r)：求解连续时间代数黎卡提方程（CARE）（最优控制的非线性方程）。
• scipy.linalg.solve_discrete_are(a,b,q,r)：求解离散时间代数黎卡提方程（CARE）。
• scipy.linalg.solve_continuous_lyapunov(a,q)：求解控制理论中的连续李雅普诺夫方程。
• scipy.linalg.solve_discrete_lyapunov(a,q)：求解控制理论中的离散李雅普诺夫方程。

三、 Scipy积分

1 对给定函数积分

1.1 通用积分

scipy.integrate.quad(func,a,b)：对函数func从a到b进行积分。

返回值为元组，元组第一个值为积分值，第二个值为绝对误差

>>> scipy.integrate.quad(np.sin,0,2*np.pi)
(2.221501482512777e-16, 4.3998892617845996e-14)

a,b可以取正负无穷大

>>> scipy.integrate.quad(lambda x:1/x**2,1,np.inf)
(1.0, 1.1102230246251565e-14)

func可以为多个参数的函数，用args参数为后边的参数指定值。

>>> scipy.integrate.quad(np.add,1,2,args=(1,))   #函数为x+y，y取值为1，即对x+1进行积分。
(2.5, 2.7755575615628914e-14)

发散积分并不会出错，而是会在命令行中返回一个警告消息。

>>> scipy.integrate.quad(lambda x:1/x,0, np.inf)
<input>:1: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
(48.720960971461565, 16.30167063049395)
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and