Python科学计算：数组💯数据生成💯数据交互💯微积分💯插值💯拟合💯FFT💯卷积💯滤波💯统计

求积

矩阵是线性代数的核心对象，是由$m$行$n$列的数组成的矩形数阵，从编程的角度理解，就是二维数组。在Numpy中，数组支持元素之间的各种运算，也支持与单个数值的各种运算，不足为奇。

另一方面，Numpy也为矩阵这种代数结构，提供了更加丰富的运算支持。由于矩阵和数组中均有维度的概念，为了不致混淆，下面提到矩阵维度是，特指$m\times n$中的$m$和$n$，而数组维度则指数组在排列时，其索引轴的个数。

numpy中提供了矩阵乘法【matmul】，内积【inner】、外积【outer】、直积【kron】等，其运算规则如下表所示。

矩阵乘法	内积	外积
$Cij=\sum a_{ik}{b}_{kj}$	$a\cdot b={\sum }_i{a}_i{b}_i$	$(AB{)}_{ij}=a_i{b}_j$

特别地，对于矩阵乘法，numpy重载了运算符@。但使用函数计算的好处是，可以指定用于计算的坐标轴，从而在数组维度大于2的时候，起到张亮乘法的作用。

所谓直积，又称克罗内克积，会对矩阵的维度进行扩张，定义如下。

$$A\otimes B=\left[\begin{array}{cccc}{a}_{11}B& a_{12}B& \cdots & a_{1n}B\\ a_{21}B& a_{22}B& \cdots & a_{2n}B\\ \cdots & \cdots & \cdots & \cdots \\ a_{m1}B& a_{m2}B& \cdots & a_{mn}B\end{array}\right]$$

求逆

对于两个方阵$A,B$，若$AB=E$，且$E$为单位阵，则$A,B$互逆，可记作$A=B^{-1},B=A^{-1}$。

numpy和scipy均提供了linalg模块，用于线性代数的相关计算，为了表述简洁，后面对这两个模块分别称为【nl】和【sl】。这两个模块均提供了求逆函数【inv】，但【sl】中的更快。二者在求逆时会有几乎可以忽略的误差，示例如下

import numpy as np
import numpy.linalg as nl
import scipy.linalg as sl

A = np.random.rand(3,3)
np.sum(np.abs(sl.inv(A)-nl.inv(A)))
# 6.106226635438361e-16

如果$A,B$不为方阵，那么久需要对逆的概念进行扩展，从而得到广义逆。对于复矩阵$A$而言，若存在复矩阵$G$，满足下列条件，其中${\ast }^H$表示共轭转置，则则称$G$为$A$的穆尔-彭罗斯广义逆，也叫加号逆，记作$G=A^{+}$。

1. $AGA=A$
2. $GAG=G$
3. $(AG{)}^H=AG$
4. $(GA{)}^H=GA$

【pinv】和【pinvh】均为【sl】中的加号逆函数，分别用于实矩阵和复矩阵。

最小二乘法

如果把线性方程组写作矩阵的形式，即$Ax=b$，若$A$是满秩方阵，则可解的$x=A^{-1}b$，否则可表示为$x=(A^{T}A{)}^{-1}{A}^{T}b$，这就是矩阵的逆在求解线性方程组时的应用。

所谓线性最小二乘法，可以理解为是解方程的延续，区别在于，当未知量远大于方程个数的时候，将得到一个无解的问题。最小二乘法的实质，是保证误差最小的情况下对未知数进行赋值。

【nl】和【sl】均实现了最小二乘法，函数名称均为【sltsq】，且必须输入的参数均为$a,b$，对应$Ax=b$中的$A,b$。二者返回值均有4个，分别是拟合得到的$x$、拟合误差、矩阵$a$的秩、以及矩阵$a$的单值形式。考虑到这两个函数在速度上几乎没有太大区别，所以有限推荐使用【nl】中的最小二乘法，示例如下。

np.random.seed(42)
M = np.random.rand(4,4)
x = np.arange(4)
y = M@x
xhat = np.linalg.lstsq(M,y)
print(xhat[0])
#[0. 1. 2. 3.]

特征值

能够描述矩阵特征的数值有很多，比如范数【norm】，迹【trace】，行列式【det】等，但被冠以特征值之名的，则只有$Ax=\lambda x$中的$\lambda$，相应地，每个$\lambda$对应一组特征向量$x$。

【sl】均提供了特征值计算函数，列表如下，其中被方括号括住的，表示【nl】中提供了同名函数。

适用情况	特征值+向量	只返回特征值
方阵	【eig】	【eigvals】
厄米矩阵	【eigh】	【eigvalsh】
厄米带状矩阵	eig_banded	eigvals_banded
对称三对角矩阵	eigh_tridiagonal	eigvalsh_tridiagonal

以方阵为例，示例如下

A = np.random.rand(3,3)

sl.eig(A)       # 返回特征值和特征向量
sl.eigvals(A)   # 只返回特征值