python scipy.linalg线性代数计算

        当使用优化的 ATLAS LAPACK 和 BLAS 库构建 SciPy 时，它具有非常快速的线性代数功能。如果你深入挖掘，所有原始 LAPACK 和 BLAS 库都可供您使用，以提高速度。所有这些线性代数例程都需要一个可以转换为二维数组的对象。这些例程的输出也是一个二维数组。

1、scipy.linalg 与 numpy.linalg

        scipy.linalg包含 numpy.linalg中的所有函数，加上其他一些不包含在numpy.linalg。scipy.linalg总是使用 BLAS/LAPACK 支持编译，而对于 numpy，这是可选的。因此，scipy 版本可能会更快，具体取决于 numpy 的安装方式。因此，除非你不想将scipy依赖项添加到numpy程序中，否则请scipy.linalg使用numpy.linalg。

2、矩阵求逆（linalg.inv()）

        矩阵的逆A是矩阵B, 这样AB=I，在那里 I 是由沿主对角线的单位组成的单位矩阵。通常，B表示B=A-1，在 SciPy 中使用linalg.inv()求得。如：

以下示例演示了 SciPy 中的这种计算：

import numpy as np
from scipy import linalg
 
A = np.array([[1, 3, 5], [2, 5, 1], [2, 3, 8]])
print("A =", A)
 
A_1 = linalg.inv(A)
print("A的逆为：", A_1)

结果：

A = [[1 3 5] [2 5 1] [2 3 8]] 
A的逆为： [[-1.48 0.36 0.88] [ 0.56 0.08 -0.36] [ 0.16 -0.12 0.04]]

3、线性方程组求解（linalg.solve）

        使用 scipy 命令求解线性方程组很简单linalg.solve。该命令需要一个输入矩阵和一个右侧向量，然后计算解向量。提供了输入对称矩阵的选项，可在适用时加快处理速度。例如，假设需要求解以下联立方程：

        但是，最好使用 linalg.solve 命令，它可以更快且数值更稳定。在这种情况下，它给出了与以下示例中所示相同的答案：

import numpy as np
from scipy import linalg
 
A = np.array([[1, 3, 5], [2, 5, 1], [2, 3, 8]])
 
A_1 = linalg.inv(A)
print("A的逆为：", A_1)
 
I = A.dot(linalg.inv(A))  # np.dot表示点乘
print("A_1 * A =I", I)

结果：

A的逆 * B: [[-1.66887417]
     [ 2.04635762]
     [ 1.10596026]]
solve(A, b): [[-1.66887417]
     [ 2.04635762]
     [ 1.10596026]]

4、计算行列式的值（linalg.det）

        方阵A的行列式表示为|A|,计算方式为：

例如：

在 SciPy 中，其计算方式如下例所示：

A = np.array([[1, 3, 5], [2, 5, 1], [2, 3, 8]]) 
A_det = linalg.det(A) 
print("|A| =", A_det)

结果：

|A| = -25.00

5、范数（linalg.norm()）

        对于向量来说，范数在大多数情况下指的就是模长，scipy提供了函数linalg.norm来计算：

B=np.array([3,4,12]) # 
print(linalg.norm(B))
#返回值为 13.0

        矩阵和向量范数也可以用 SciPy 计算。使用 的 order 参数的不同参数可以使用范围广泛的规范定义linalg.norm。此函数采用 rank-1（向量）或 rank-2（矩阵）数组和可选的 order 参数（默认为 2）。基于这些输入，计算请求顺序的向量或矩阵范数。

        对于向量x，顺序参数可以是任何实数，包括 inf或-inf。计算的范数是：

        对于矩阵A, norm 的唯一有效值是+-2, +-1,+-in 和 'fro'（或 'f'）因此，

例子：

import numpy as np
>>> from scipy import linalg
>>> A=np.array([[1,2],[3,4]])
>>> A
array([[1, 2],
      [3, 4]])
>>> linalg.norm(A)
5.4772255750516612
>>> linalg.norm(A,'fro') # frobenius norm is the default
5.4772255750516612
>>> linalg.norm(A,1) # L1 norm (max column sum)
6
>>> linalg.norm(A,-1)
4
>>> linalg.norm(A,np.inf) # L inf norm (max row sum)
7

6、解决最小二乘问题和伪逆问题

        optimize函数含有实现最小二乘法的函数 leastsq，如下通过对正弦函数的拟合，求得最小二乘拟合参数。func三参数A,k,theta分别表示对应振幅，频率，相角。

例子：

        1、调用leastsq进行数据拟合

        2、residuals为计算误差的函数

        3、p0为拟合参数的初始值

        4、args为需要拟合的实验数据

import numpy as np
import pylab as pl
from pylab import mpl
from scipy.optimize import leastsq
 
mpl.rcParams['font.sans-serif'] = ['KaiTi']  # 解决中文乱码
mpl.rcParams['axes.unicode_minus'] = False  # 解决负号显示为方框的问题
 
def func(x, p):
    # 数据拟合所用函数： A*sin(2*pi*k*x + theta)
    A, k, theta = p
    return A * np.sin(2 * np.pi * k * x + theta)
 
def residuals(p, y, x):
    # 实验数据x,y和拟合函数之间的差，p为拟合需要找到的系数
    return y - func(x, p)
 
x = np.linspace(0, -2 * np.pi, 100)  # 创建等差数列，100表示数据点个数
A, k, theta = 10, 0.34, np.pi / 6  # 真实数据的函数参数
y0 = func(x, [A, k, theta])  # 真实数据
y1 = y0 + 2 * np.random.randn(len(x))  # 加入噪声后的实验数据
 
p0 = [7, 0.2, 0]  # 第一次猜测的函数拟合参数
plsq = leastsq(residuals, p0, args=(y1, x))  # 最小二乘拟合
print(u"真实参数：", [A, k, theta])
print(u"拟合参数：", plsq[0])  # 实验数据拟合后的参数
 
# 作图
pl.plot(x, y0, label=u'真实数据')
pl.plot(x, y1, label=u'带噪声的实验数据')
pl.plot(x, func(x, plsq[0]), label=u"拟合数据")
pl.legend()
pl.show()

结果：

广义逆：

        linalg.pinv使用命令或计算广义逆 linalg.pinv2。这两个命令在计算广义逆的方式上有所不同。第一个使用 linalg.lstsq 算法，而第二个使用奇异值分解。让A为M*N矩阵，那么如果M>N，广义逆是：

        而如果M：

        在这种情况下M=N，然后：

        只要A是可逆的。

参考：

1、SciPy 用户指南

SciPy User Guide — SciPy v1.8.0 Manual

2、Python scipy.linalg.norm用法及代码示例

Python scipy.linalg.norm用法及代码示例 - 纯净天空