精通 SciPy（一）
零、前言
编写《精通 SciPy》的想法出现了，但是在发表《数值和科学计算学习》之后两个月。 在南卡罗来纳大学介绍这本书的过程中，我荣幸地向工程师，科学家和学生的不同受众介绍了本书的内容，他们每个人都有非常不同的研究问题和他们自己的一套首选计算资源 。 在演讲之后的几周中，我帮助了一些专业人士过渡到基于 SciPy 的环境。 在这些会议中，我们讨论了 SciPy 在幕后如何与他们已经使用的相同算法集（通常是相同代码）。 我们尝试了一些示例，并系统地获得了可比的性能。 我们立即看到基于健壮的脚本语言的通用环境的明显好处。 通过 SciPy 堆栈，我们发现了一种与同事，学生或雇主进行交流和共享结果的简便方法。 在所有情况下，切换到 SciPy 堆栈都为我们的小组提供了更快的设置，使新手可以快速掌握。
参与此过程的每个人都从新手变为高级用户，最终很快掌握了 SciPy 堆栈。 在大多数情况下，与我一起工作的个人的科学背景使过渡变得无缝。 当他们能够将他们的研究背后的理论与所提供的解决方案进行对比时，掌握过程就变成了现实。 啊哈时刻总是在复制过程中得到仔细的指导和解释的过程中发生。
这恰恰是这本书背后的哲学。 我邀请您参加类似的会议。 每章都被设想为与具有一定科学需求的个人进行对话，这些需求以数值计算表示。 我们在一起发现了相关的例子，这些例子是解决这些问题的不同可能方法，其背后的理论以及每种方法的利弊。
写作过程遵循类似的路径，以获取引人入胜的示例集。 我与几个不同领域的同事进行了交谈。 每个部分都清楚地反映了这些交流。 在编写最具挑战性的章节（最后四个章节）时，这一点至关重要。 为了确保整本书具有相同的质量，并且始终尝试遵循一组严格的标准，需要花更长的时间才能使这些章节满意。 特别要提到的是 NASA Langley 研究中心的 Aaron Dutle，他帮助塑造了计算几何学这一章的一部分； Facebook 的数据分析师 Parsa Bakhtary，激发了本章中关于统计计算在数据分析中的应用的许多技术。 。
这是一次了不起的旅程，有助于加深我对数值方法的理解，拓宽了我在解决问题上的视野，并增强了我的科学成熟度。 我希望它对您有相同的影响。
本书涵盖的内容
第 1 章和“数值线性代数”概述了矩阵在解决科学计算中的作用。 对于理解后续章节的大多数过程和思想而言，这是至关重要的章节。 您将学习如何在 Python 中有效地构造和存储大型矩阵。 然后，我们着手研究对其的基本操作和操作，然后进行因式分解，矩阵方程的解以及特征值/特征向量的计算。
第 2 章，“插值和近似”开发了近似函数的先进技术，并将其应用于科学计算。 这是接下来两章的 segway。
第 3 章，“微分和积分”探索了产生函数导数的不同技术，更重要的是，如何通过积分过程有效地计算面积和体积。 这是两章专门介绍科学计算中数值方法核心的第一章。 第二部分也是第 5 章和“常微分方程的初值问题”的介绍，其中提到了常微分方程。
第 4 章，“非线性方程式和优化”是一门技术性很强的章节，在其中我们讨论根据所涉及函数的种类获得函数系统的根和极值的最佳方法。
第 5 章和“常微分方程的初值问题”是有关实际问题的五个章节中的第一章。 通过示例，我们向您展示解决微分方程组的最流行技术以及一些应用。
第 6 章和“计算几何”浏览了计算机科学这一领域中最重要的算法。
第 7 章和“描述性统计”是关于统计计算及其在数据分析中的应用的两章中的第一章。 在本章中，我们重点介绍概率和数据探索。
第 8 章，“推断和数据分析”是有关数据分析的第二章。 我们专注于统计推断，机器学习和数据挖掘。
第 9 章，“数字图像处理”是本书的最后一章。 在其中，我们探索了图像压缩，编辑，还原和分析的技术。
这本书需要什么
要处理这些示例并尝试使用本书的代码，您所需要的只是带有 SciPy 堆栈的 Python 的最新版本（2.7 或更高版本）：NumPy，SciPy 库，matplotlib，IPython，pandas 和 SymPy。 尽管整本书提供了独立安装所有这些程序的方法，但是我们建议您通过科学的 Python 发行版（例如 Anaconda）执行全局安装。
这本书是给谁的
尽管这本书和技术最终旨在供应用数学家，工程师和计算机科学家使用，但此处介绍的材料面向的是更广泛的读者。 所需要做的只是精通 Python，熟悉 iPython，对科学计算中的数值方法有一定的了解，以及对在科学，工程学或数据分析中开发严肃应用程序的浓厚兴趣。
约定
在本书中，您将找到许多可以区分不同类型信息的文本样式。 以下是这些样式的一些示例，并解释了其含义。
文本中的代码字，数据库表名称，文件夹名称，文件名，文件扩展名，路径名，伪 URL，用户输入和 Twitter 句柄如下所示：“我们可以通过使用 include 指令来包含其他上下文。”
任何命令行输入或输出的编写方式如下：
In [7]: %time eigvals, v = spspla.eigsh(A, 5, which='SM')
CPU times: user 19.3 s, sys: 532 ms, total: 19.8 s
Wall time: 16.7 s
In [8]: print eigvals
[ 10.565523 10.663114 10.725135 10.752737 10.774503]
注意
警告或重要提示会出现在这样的框中。
提示
提示和技巧如下所示。
一、数值线性代数
术语“数值线性代数”是指使用矩阵来解决计算科学问题。 在本章中，我们首先学习如何在 Python 中有效地构造这些对象。 我们强调从在线存储库中导入大型稀疏矩阵。 然后，我们继续审查它们的基本操作和操作。 下一步是研究 SciPy 中实现的不同矩阵函数。 我们将继续探索矩阵方程解，特征值及其相应特征向量的计算的不同因式分解。
动机
下图显示了代表一系列网页（从 1 到 8）的图形：
从一个节点到另一个节点的箭头表示存在从发送节点代表的网页到接收节点代表页面的链接。 例如，从节点 2 到节点 1 的箭头表示网页 2 中有指向网页 1 的链接。请注意，网页 4 如何具有两个外部链接（至页面 2 和 8），并且有三个页面链接至网页 4（第 2，6 和 7 页）。 由节点 2，4 和 8 表示的页面乍一看似乎最受欢迎。
有没有一种数学方法可以真正表达网络在网络中的受欢迎程度？ Google 的研究人员提出了PageRank的想法，可以通过计算页面链接的数量和质量来大致估计该概念。 它是这样的：
-
我们以以下方式构造此图的过渡矩阵
T={a[i,j]}：如果存在从网页i到网页j的链接，则条目a[i,j]为1/k，并且总数目为 网页i中的外部链接等于k。 否则，该条目仅为零。N个网页的转换矩阵的大小始终为N × N。 在我们的例子中，矩阵的大小为8×8：0 1/2 0 0 0 0 0 0 1 0 1/2 1/2 0 0 0 0 0 0 0 0 0 0 1/3 0 0 1/2 0 0 0 1 1/3 0 0 0 1/2 0 0 0 0 0 0 0 0 0 0 0 0 1/2 0 0 0 0 1/2 0 0 1/2 0 0 0 1/2 1/2 0 1/3 0
让我们打开一个 iPython 会话并将此特定矩阵加载到内存中。
注意
请记住，在 Python 中，索引从零开始，而不是一个。
In [1]: import numpy as np, matplotlib.pyplot as plt, \
...: scipy.linalg as spla, scipy.sparse as spsp, \
...: scipy.sparse.linalg as spspla
In [2]: np.set_printoptions(suppress=True, precision=3)
In [3]: cols = np.array([0,1,1,2,2,3,3,4,4,5,6,6,6,7,7]); \
...: rows = np.array([1,0,3,1,4,1,7,6,7,3,2,3,7,5,6]); \
...: data = np.array([1., 0.5, 0.5, 0.5, 0.5, \
...: 0.5, 0.5, 0.5, 0.5, 1., \
...: 1./3, 1./3, 1./3, 0.5, 0.5])
In [4]: T = np.zeros((8,8)); \
...: T[rows,cols] = data
从过渡矩阵，我们通过在 0 和 1 之间固定一个正常数p并按照公式来创建PageRank矩阵G = (1-p) * T + p * B以获得合适的阻尼系数p。 在此，B是具有与T相同大小的矩阵，其所有条目均等于1/N。 例如，如果我们选择p = 0.15，我们将获得以下PageRank矩阵：
In [5]: G = (1-0.15) * T + 0.15/8; \
...: print G
[[ 0.019 0.444 0.019 0.019 0.019 0.019 0.019 0.019]
[ 0.869 0.019 0.444 0.444 0.019 0.019 0.019 0.019]
[ 0.019 0.019 0.019 0.019 0.019 0.019 0.302 0.019]
[ 0.019 0.444 0.019 0.019 0.019 0.869 0.302 0.019]
[ 0.019 0.019 0.444 0.019 0.019 0.019 0.019 0.019]
[ 0.019 0.019 0.019 0.019 0.019 0.019 0.019 0.444]
[ 0.019 0.019 0.019 0.019 0.444 0.019 0.019 0.444]
[ 0.019 0.019 0.019 0.444 0.444 0.019 0.302 0.019]]
PageRank矩阵具有一些有趣的属性：
- 1 是复数一的特征值。
- 1 实际上是最大的特征值； 所有其他特征值的模量均小于 1。
- 对应于特征值 1 的特征向量具有所有正项。 特别地，对于特征值 1，存在一个唯一的特征向量，其项的总和等于 1。这就是我们所说的
PageRank向量。
scipy.linalg.eig的快速计算为我们找到了特征向量：
In [6]: eigenvalues, eigenvectors = spla.eig(G); \
...: print eigenvalues
[ 1.000+0.j -0.655+0.j -0.333+0.313j -0.333-0.313j –0.171+0.372j -0.171-0.372j 0.544+0.j 0.268+0.j ]
In [7]: PageRank = eigenvectors[:,0]; \
...: PageRank /= sum(PageRank); \
...: print PageRank.real
[ 0.117 0.232 0.048 0.219 0.039 0.086 0.102 0.157]
这些值对应于图表上描述的八个网页中每个网页的PageRank。 如预期的那样，这些值的最大值与第二个网页（0.232）相关联，紧随其后的是第四个网页（0.219），然后是第八个网页（0.157）。 这些值向我们提供了我们正在寻找的信息：第二个网页最受欢迎，其次是第四个，然后是八个。
注意
请注意，此网页网络问题已被转换为数学对象，成为涉及矩阵，特征值和特征向量的等效问题，并已通过线性代数技术得以解决。
转换矩阵稀疏：其大多数条目为零。 在数值线性代数中，具有极大尺寸的稀疏矩阵特别重要，这不仅是因为它们编码具有挑战性的科学问题，而且还因为很难用基本算法来操纵它们。
与其将矩阵中的所有值存储到内存中，不如仅收集非零值，并使用利用这些智能存储方案的算法，这才有意义。 内存管理的好处是显而易见的。 这些方法通常对于此类矩阵更快，并且舍入误差较小，因为通常所涉及的运算要少得多。 这是 SciPy 的另一个优点，因为它包含许多程序来解决以这种方式存储数据的不同问题。 让我们通过另一个示例观察其力量：
佛罗里达大学稀疏矩阵集合是可在线访问的最大矩阵数据库。 截至 2014 年 1 月，它包含 157 组来自各种科学学科的矩阵。 矩阵的大小从非常小的（1×2）到非常大的（2800 万×2800 万）不等。 由于它们出现在不同的工程问题中，因此预计将不断添加更多的矩阵。
提示
有关此数据库的更多信息，请参见数学软件，第 1 卷，。 T.A.，38，Issue 1，2011，pp 1：1-1：25 戴维斯（Davis）和 Y.Hu，或在上在线访问 http://www.cise.ufl.edu/research/sparse/matrices/ 。
例如，数据库中矩阵最多的组是原始的 Harwell-Boeing Collection，具有 292 种不同的稀疏矩阵。 该组也可以在 Matrix Market 上在线访问： math.nist.gov/MatrixMarke… 。
数据库中的每个矩阵都有三种格式：
- Matrix Market Exchange 格式【Boisvert 等 1997】
- Rutherford-Boeing Exchange 格式【Duff 等 1997】
- 专有 Matlab 格式
.mat。
让我们从集合中以矩阵市场交易格式将两个矩阵导入到我们的 iPython 会话中，这两个矩阵打算用于最小二乘问题的解决方案。 这些矩阵位于 www.cise.ufl.edu/research/sp… 。这些数值对应于在 Sonata 1.5-T 扫描仪（Siemens，Erlangen）上获得的幻像数据。 ，德国）使用磁共振图像处理（ MRI ）设备。 被测物体是用几种金属物体制成的人体头部的模拟。 我们下载相应的 tar 捆绑包并将其解压缩以获取两个 ASCII 文件：
mri2.mtx（最小二乘问题中的主矩阵）mri2_b.mtx（方程式的右侧）
文件mri2.mtx的前 20 行内容如下：
前十六行是注释，为我们提供了有关矩阵生成的一些信息。
- 出现的计算机视觉问题：MRI 重建
- 作者信息：UCSD 的 Mark Bydder
- 应用于数据的过程：解决最小二乘问题
Ax - b，并对结果进行后可视化
第十七行指示矩阵的大小63240行×147456列，以及数据中的非零条目数569160。
文件的其余部分恰好包括 569160 行，每行包含两个整数和一个浮点数：这些是矩阵中非零元素的位置以及相应的值。
提示
我们需要考虑到这些文件使用 FORTRAN 约定，从 1 开始而不是从 0 开始数组。
将文件读入ndarray的一个好方法是通过 NumPy 中的函数loadtxt。 然后，我们可以使用scipy在模块scipy.sparse中使用函数coo_matrix将数组转换为稀疏矩阵（coo 代表坐标内部格式）。
In [8]: rows, cols, data = np.loadtxt("mri2.mtx", skiprows=17, \
...: unpack=True)
In [9]: rows -= 1; cols -= 1;
In [10]: MRI2 = spsp.coo_matrix((data, (rows, cols)), \
....: shape=(63240,147456))
可视化此矩阵稀疏性的最佳方法是借助模块matplotlib.pyplot中的例程spy。
In [11]: plt.spy(MRI2); \
....: plt.show()
我们获得以下图像。 每个像素对应于矩阵中的一个条目； 白色表示零值，非零值根据其大小（越大，越深）用不同的蓝色阴影表示：
这是第二个文件mri2_b.mtx的前十行，它不代表稀疏矩阵，而是列向量：
%% MatrixMarket matrix array complex general
%-------------------------------------------------------------
% UF Sparse Matrix Collection, Tim Davis
% http://www.cise.ufl.edu/research/sparse/matrices/Bydder/mri2
% name: Bydder/mri2 : b matrix
%-------------------------------------------------------------
63240 1
-.07214859127998352 .037707749754190445
-.0729086771607399 .03763720765709877
-.07373382151126862 .03766685724258423
这是六条带信息的带注释的行，另外一行表示向量的形状（63240行和1列），其余各行包含两列浮点值，即浮点值的实部和虚部 相应的数据。 我们继续将这个向量读取到内存中，解决了建议的最小二乘问题，并获得了以下表示模拟人类头部切片的重构：
In [12]: r_vals, i_vals = np.loadtxt("mri2_b.mtx", skiprows=7,
....: unpack=True)
In [13]: %time solution = spspla.lsqr(MRI2, r_vals + 1j*i_vals)
CPU times: user 4min 42s, sys: 1min 48s, total: 6min 30s
Wall time: 6min 30s
In [14]: from scipy.fftpack import fft2, fftshift
In [15]: img = solution[0].reshape(384,384); \
....: img = np.abs(fftshift(fft2(img)))
In [16]: plt.imshow(img); \
....: plt.show()
提示
如果对创建此矩阵背后的理论以及该问题的详细信息感兴趣，请阅读 H.Sedarat 和 D.G.Nishimura 撰写的文章关于网格重构算法的最优性，该文章发表于 IEEE Trans。 医学图像处理，第一卷。 19 号 4，第 306-317 页，2000 年。
对于结构良好的矩阵，这些矩阵将专门涉及矩阵乘法，通常可以用智能方式存储对象。 让我们考虑一个例子。
水平地震振荡会影响高层建筑的各个楼层，具体取决于楼层振荡的固有频率。 如果我们做出某些假设，可以通过竞争获得作为 N 微分方程的二阶系统的模型，该模型可以量化具有 N 层的建筑物的振动：牛顿第二力定律 设定为等于胡克力定律和地震波产生的外力之和。
这些是我们将需要的假设：
- 每个楼层都被视为位于其质量中心的质量点。 楼层有质量
m[1], m[2], ..., m[N]。 - 每个楼层都通过线性恢复力（胡克
-k * elongation）恢复到其平衡位置。 楼层的胡克常数为k[1], k[2], ..., k[N]。 - 表示楼层振动的质量的位置是
x[1], x[2], ..., x[N]。 我们假设它们都是时间函数，并且在平衡时，它们都等于零。 - 为了简化说明，我们将假设没有摩擦：楼层上的所有阻尼效果都将被忽略。
- 楼层的方程式仅取决于相邻楼层。
将质量矩阵M设置为对角矩阵，在其对角线上包含地板质量。 将K（胡克矩阵）设置为具有以下结构的三对角矩阵，对于j每行，除以下各项外，所有条目均为零：
j-1列（我们将其设置为k[j+1]）j列（我们设置为-k[j+1]-k[j+1]），以及- 列
j+1，我们将其设置为k[j+2]。
将H设置为包含地震引起的各层外力的列向量，将X设置为包含函数x[j]的列向量。
然后我们有了系统：M * X''= K * X + H。 该系统的齐次部分是M与K的倒数的乘积，我们将其表示为A。
为了解决齐次线性二阶系统X''= A * X，我们定义变量Y包含2 * N项：所有N函数x[j]，然后是它们的派生词x'[j]。 该二阶线性系统的任何解都具有一阶线性系统的相应解Y'= C * Y，其中C是大小为2 * N×2 * N的块矩阵 。 此矩阵C由大小为N × N的仅包含零的块组成，其后为单位（大小为N × N），并且在这两个下方，矩阵A在水平方向后跟另一个N×N零块。
不必将此矩阵C存储到内存或其任何因素或块中。 相反，我们将利用其结构，并使用线性运算符来表示它。 然后需要最少的数据来生成该运算符（仅质量值和胡克系数），远小于其任何矩阵表示形式。
让我们展示一个六层楼的具体例子。 我们首先指出它们的质量和胡克常数，然后继续构建A的表示形式为线性算子：
In [17]: m = np.array([56., 56., 56., 54., 54., 53.]); \
....: k = np.array([561., 562., 560., 541., 542., 530.])
In [18]: def Axv(v):
....: global k, m
....: w = v.copy()
....: w[0] = (k[1]*v[1] - (k[0]+k[1])*v[0])/m[0]
....: for j in range(1, len(v)-1):
....: w[j] = k[j]*v[j-1] + k[j+1]*v[j+1] - \
....: (k[j]+k[j+1])*v[j]
....: w[j] /= m[j]
....: w[-1] = k[-1]*(v[-2]-v[-1])/m[-1]
....: return w
....:
In [19]: A = spspla.LinearOperator((6,6), matvec=Axv, matmat=Axv,
....: dtype=np.float64)
现在C的结构非常简单（比其矩阵要简单得多！）：
In [20]: def Cxv(v):
....: n = len(v)/2
....: w = v.copy()
....: w[:n] = v[n:]
....: w[n:] = A * v[:n]
....: return w
....:
In [21]: C = spspla.LinearOperator((12,12), matvec=Cxv, matmat=Cxv,
....: dtype=np.float64)
这个齐次系统的解以C的指数作用形式出现：Y(t) = expm(C * t) * Y(0)，其中expm()表示矩阵指数函数。 在 SciPy 中，使用模块scipy.sparse.linalg中的例程expm_multiply执行此操作。
例如，在我们的情况下，给定包含值x[1](0)=0, ..., x[N](0)=0, x'[1](0)=1, ..., x'[N](0)=1的初始值，如果我们需要以大小0.1为步长在 0 和 1 之间的t值求解Y(t)，则可以发出以下命令：
提示
据报道，在某些安装中，在下一步中，必须给出 C 的矩阵，而不是实际的线性算子（因此与手册相矛盾）。 如果在您的系统中是这种情况，只需将下一行中的 C 更改为其矩阵表示即可。
In [22]: initial_condition = np.zeros(12); \
....: initial_condition[6:] = 1
In [23]: Y = spspla.exp_multiply(C, np.zeros(12), start=0,
....: stop=1, num=10)
然后可以计算出第一层中六层的振动，并绘制出曲线图。 例如，要查看一楼的振荡情况，我们可以发出以下命令：
In [24]: plt.plot(np.linspace(0,1,10), Y[:,0]); \
....: plt.xlabel('time (in seconds)'); \
....: plt.ylabel('oscillation')
我们获得以下图。 请注意，第一层如何在第一十分之一秒内上升，但从原来的高度从 0.1 秒下降到 0.9 秒，几乎下降到一米以下，然后开始缓慢上升：
提示
有关微分方程系统以及如何使用指数函数求解它们的更多详细信息，请阅读例如《 HTHT0]基本微分方程》第 10 版。 ，由 William E. Boyce 和 Richard C. DiPrima 撰写。 威利（Wiley），2012 年。
这三个示例说明了第一章“数值线性代数”的目标。 在 Python 中，这首先要通过以下任何类将数据存储为矩阵形式或作为相关的线性运算符来实现：
numpy.ndarray（确保它们是二维的）numpy.matrixscipy.sparse.bsr_matrix（块稀疏行矩阵）scipy.sparse.coo_matrix（坐标格式的稀疏矩阵）scipy.sparse.csc_matrix（压缩的稀疏列矩阵）scipy.sparse.csr_matrix（压缩的稀疏行矩阵）scipy.sparse.dia_matrix（带有对角存储的稀疏矩阵）scipy.sparse.dok_matrix（基于字典的稀疏矩阵）scipy.sparse.lil_matrix（基于链表的稀疏矩阵）scipy.sparse.linalg.LinearOperator
正如我们在示例中所看到的，不同类的选择主要遵循数据的稀疏性以及我们将应用于它们的算法。
提示
在以下各节中，我们将学习何时应用这些选择。
然后，该选择决定了我们用于不同算法的模块：scipy.linalg用于通用矩阵，而scipy.sparse和scipy.sparse.linalg都用于稀疏矩阵或线性运算符。 这三个 SciPy 模块是在高度优化的计算机库 BLAS（用 Fortran77 编写），LAPACK（在 Fortran90 中编写），ARPACK（在 Fortran77 中）和 SuperLU（在 C 中）的基础上编译的。
注意
为了更好地理解这些基础软件包，请阅读其创建者的描述和文档：
- BLAS： netlib.org/blas/faq.ht…
- LAPACK： netlib.org/lapack/lapa…
- ARPACK： www.caam.rice.edu/software/AR…
- SuperLU： crd-legacy.lbl.gov/~xiaoye/Sup…
这三个 SciPy 模块中的大多数例程都是上述库中函数的包装。 如果我们愿意，我们还可以直接调用基础函数。 在scipy.linalg模块中，我们具有以下内容：
scipy.linalg.get_blas_funcs从 BLAS 调用例程scipy.linalg.get_lapack_funcs从 LAPACK 调用例程
例如，如果我们要使用BLAS函数NRM2计算 Frobenius 范数：
In [25]: blas_norm = spla.get_blas_func('nrm2')
In [26]: blas_norm(np.float32([1e20]))
Out[26]: 1.0000000200408773e+20
创建矩阵和线性算子
在本章的第一部分中，我们将专注于有效创建矩阵。 我们首先回顾一些不同的方法来构造基本矩阵作为ndarray实例类，包括对 NumPy 和 SciPy 中已经包含的所有特殊矩阵的枚举。 我们继续研究从基本矩阵构造复杂矩阵的可能性。 我们在matrix实例类中回顾了相同的概念。 接下来，我们详细探讨输入稀疏矩阵的不同方法。 我们以线性算子的构造结束本节。
注意
我们假定您熟悉 NumPy 中的ndarray创建以及数据类型（dtype），索引，两个或多个数组组合的例程，数组操作或从这些对象中提取信息。 在本章中，我们将重点介绍仅对矩阵有意义的函数，方法和例程。 如果它们的输出没有转换为线性代数等价物，我们将不考虑其运算。 有关ndarray的入门知识，建议您浏览《数值科学计算学习》第二版的第二章。 为了快速回顾线性代数，我们推荐 Hoffman 和 Kunze，《线性代数第二版》，Pearson，1971 年。
在 ndarray 类中构造矩阵
我们可以通过三种不同的方式从作为ndarray实例的数据创建矩阵：从标准输入手动创建，通过为每个条目分配函数值或通过从外部文件检索数据来创建矩阵。
建设者
|
描述
|
| --- | --- |
| numpy.array(object) | 从object创建矩阵 |
| numpy.diag(arr, k) | 在对角线k上创建数组arr的条目的对角矩阵 |
| numpy.fromfunction(function, shape) | 通过在每个坐标上执行函数来创建矩阵 |
| numpy.fromfile(fname) | 从文本或二进制文件创建矩阵（基本） |
| numpy.loadtxt(fname) | 从文本文件创建矩阵（高级） |
让我们创建一些示例矩阵来说明上表中定义的一些函数。 和以前一样，我们开始一个 iPython 会话：
In [1]: import numpy as np, matplotlib.pyplot as plt, \
...: scipy.linalg as spla, scipy.sparse as spsp, \
...: scipy.sparse.linalg as spspla
In [2]: A = np.array([[1,2],[4,16]]); \...: A
Out[2]:
array([[ 1, 2],
[ 4, 16]])
In [3]: B = np.fromfunction(lambda i,j: (i-1)*(j+1),
...: (3,2), dtype=int); \
...: print B
...:
[[-1 -2]
[ 0 0]
[ 1 2]]
In [4]: np.diag((1j,4))
Out[4]:
array([[ 0.+1.j, 0.+0.j],
[ 0.+0.j, 4.+0.j]])
具有预定零和一的特殊矩阵可以使用以下函数构造：
|建设者
|
描述
|
| --- | --- |
| numpy.empty(shape) | 给定形状的数组，条目未初始化 |
| numpy.eye(N, M, k) | 二维数组，对角线为k - 1，其他位置为零 |
| numpy.identity(n) | 身份阵列 |
| numpy.ones(shape) | 所有条目等于 1 的数组 |
| numpy.zeros(shape) | 所有条目等于零的数组 |
| numpy.tri(N, M, k) | 在给定对角线处及以下的数组，否则为零 |
提示
除了numpy.tri以外，所有这些构造都具有一个伴随函数xxx_like，该函数创建具有所请求特征且具有与另一个源ndarray类相同的形状和数据类型的ndarray：
In [5]: np.empty_like(A)
Out[5]:
array([[140567774850560, 140567774850560],
[ 4411734640, 562954363882576]])
值得注意的是构造为数值范围的数组。
|建设者
|
描述
|
| --- | --- |
| numpy.arange(stop) | 间隔内的均等值 |
| numpy.linspace(start, stop) | 在一定间隔内均匀间隔的数字 |
| numpy.logspace(start, stop) | 对数刻度上等距的数字 |
| numpy.meshgrid | 来自两个或多个坐标向量的坐标矩阵 |
| numpy.mgrid | nd_grid实例返回密集多维meshgrid |
| numpy.ogrid | nd_grid实例返回开放的多维meshgrid |
可以在 NumPy 和模块scipy.linalg中轻松调用在线性代数中具有众多应用的特殊矩阵。
建设者
|
描述
|
| --- | --- |
| scipy.linalg.circulant(arr) | 一维数组arr生成的循环矩阵 |
| scipy.linalg.companion(arr) | arr编码的多项式的伴随矩阵 |
| scipy.linalg.hadamard(n) | Sylvester 构造大小为n × n的 Hadamard 矩阵。n必须为 2 的幂 |
| scipy.linalg.hankel(arr1, arr2) | 第一列为arr1，最后一列为arr2的汉克矩阵 |
| scipy.linalg.hilbert(n) | 大小为n × n的希尔伯特矩阵 |
| scipy.linalg.invhilbert(n) | n × n大小的希尔伯特矩阵的逆 |
| scipy.linalg.leslie(arr1, arr2) | 具有繁殖力数组arr1和生存系数arr2的莱斯利矩阵 |
| scipy.linalg.pascal(n) | n × n截断了二项式系数的 Pascal 矩阵 |
| scipy.linalg.toeplitz(arr1, arr2) | 第一列arr1和第一行arr2的 Toeplitz 数组 |
| numpy.vander(arr) | 数组arr的范德蒙德矩阵 |
例如，一种获得所有数量级的最大二项式系数（对应的帕斯卡三角形）的快速方法是借助精确的帕斯卡矩阵。 下面的示例显示如何计算这些系数，直到13为止：
In [6]: print spla.pascal(13, kind='lower')
除了这些基本的构造函数之外，我们始终可以以不同的方式堆叠数组：
|建设者
|
描述
|
| --- | --- |
| numpy.concatenate((A1, A2, ...)) | 结合矩阵 |
| numpy.hstack((A1, A2, ...)) | 水平堆叠矩阵 |
| numpy.vstack((A1, A2, ...)) | 垂直堆叠矩阵 |
| numpy.tile(A, reps) | 重复矩阵一定次数（由reps决定） |
| scipy.linalg.block_diag(A1,A2, ...) | 创建块对角线数组 |
让我们观察一下其中的一些构造函数：
In [7]: np.tile(A, (2,3)) # 2 rows, 3 columns
Out[7]:
array([[ 1, 2, 1, 2, 1, 2],
[ 4, 16, 4, 16, 4, 16],
[ 1, 2, 1, 2, 1, 2],
[ 4, 16, 4, 16, 4, 16]])
In [8]: spla.block_diag(A,B)
Out[8]:
array([[ 1, 2, 0, 0],
[ 4, 16, 0, 0],
[ 0, 0, -1, -2],
[ 0, 0, 0, 0],
[ 0, 0, 1, 2]])
在矩阵类中构造矩阵
对于matrix类，直接创建矩阵的通常方法是调用numpy.mat或numpy.matrix。 在创建类似于A的矩阵时，请观察numpy.matrix的语法比numpy.array的语法舒适得多。 使用此语法，用逗号分隔的不同值属于矩阵的同一行。 分号表示行的更改。 还要注意强制转换为matrix类！
In [9]: C = np.matrix('1,2;4,16'); \
...: C
Out[9]:
matrix([[ 1, 2],
[ 4, 16]])
这两个功能还将任何ndarray转换为matrix。 第三个功能可以完成此任务：numpy.asmatrix：
In [10]: np.asmatrix(A)
Out[10]:
matrix([[ 1, 2],
[ 4, 16]])
对于由块组成的矩阵的排列，除了前面描述的ndarray的常见堆栈操作外，我们还具有非常方便的功能numpy.bmat。 请注意与numpy.matrix的语法相似，尤其是使用逗号表示水平串联，使用分号表示垂直串联：
In [11]: np.bmat('A;B') In [12]: np.bmat('A,C;C,A')
Out[11]: Out[12]:
matrix([[ 1, 2], matrix([[ 1, 2, 1, 2],
[ 4, 16], [ 4, 16, 4, 16],
[-1, -2], [ 1, 2, 1, 2],
[ 0, 0], [ 4, 16, 4, 16]])
[ 1, 2]])
构造稀疏矩阵
有七种输入稀疏矩阵的方法。 每种格式旨在使特定问题或操作更有效。 让我们详细了解一下它们：
|方法
|
Name
|
最佳使用
|
| --- | --- | --- |
| BSR | 块稀疏行 | 如果矩阵包含块，则为高效算术。 |
| COO | 坐标 | 快速高效的施工格式。 转换为 CSC 和 CSR 格式的有效方法。 |
| CSC | 压缩稀疏柱 | 高效的矩阵算法和列切片。 矩阵向量乘积相对快。 |
| CSR | 压缩稀疏行 | 高效的矩阵算法和行切片。 执行矩阵向量乘积最快。 |
| DIA | 对角线存储 | 如果矩阵包含非零项的长对角线，则对于构造和存储有效。 |
| DOK | 钥匙字典 | 高效的增量构造和单个矩阵条目的访问。 |
| LIL | 基于行的链接列表 | 切片灵活。 有效更改矩阵稀疏度。 |
最多可以用五种方式填充它们，每种稀疏矩阵格式中共有三种：
-
它们可以强制转换为稀疏任何通用矩阵。
lil格式是使用此方法最有效的方法：In [13]: A_coo = spsp.coo_matrix(A); \ ....: A_lil = spsp.lil_matrix(A) -
它们可以将另一种稀疏格式的另一种稀疏矩阵转换为特定的稀疏格式：
In [14]: A_csr = spsp.csr_matrix(A_coo) -
可以通过指示形状和
dtype来构造任意形状的空稀疏矩阵：In [15]: M_bsr = spsp.bsr_matrix((100,100), dtype=int)
它们都有几种不同的额外输入法，每种输入法都特定于其存储格式。
-
花式索引：就像处理任何通用矩阵一样。 这仅适用于 LIL 或 DOK 格式：
In [16]: M_lil = spsp.lil_matrix((100,100), dtype=int) In [17]: M_lil[25:75, 25:75] = 1 In [18]: M_bsr[25:75, 25:75] = 1 NotImplementedError Traceback (most recent call last) <ipython-input-18-d9fa1001cab8> in <module>() ----> 1 M_bsr[25:75, 25:75] = 1 [...]/scipy/sparse/bsr.pyc in __setitem__(self, key, val) 297 298 def __setitem__(self,key,val): --> 299 raise NotImplementedError 300 301 ###################### NotImplementedError: -
键字典：当我们一次创建，更新或搜索每个元素时，此输入系统最有效。 它仅对 LIL 和 DOK 格式有效：
In [19]: M_dok = spsp.dok_matrix((100,100), dtype=int) In [20]: position = lambda i, j: ((i<j) & ((i+j)%10==0)) In [21]: for i in range(100): ....: for j in range(100): ....: M_dok[i,j] = position(i,j) ....: -
Data, rows, and columns: This is common to four formats: BSR, COO, CSC, and CSR. This is the method of choice to import sparse matrices from the Matrix Market Exchange format, as illustrated at the beginning of the chapter.
提示
使用数据，行和列输入法时，最好在构造中始终包含选项
shape。 如果未提供此值，将从行和列的最大坐标中推断出矩阵的大小，从而可能导致矩阵的大小小于所需的大小。 -
Data, indices, and pointers: This is common to three formats: BSR, CSC, and CSR. It is the method of choice to import sparse matrices from the Rutherford-Boeing Exchange format.
注意
Rutherford-Boeing Exchange 格式是 Harwell-Boeing 格式的更新版本。 它将矩阵存储为三个向量：
pointers_v，indices_v和data。j列的条目的行索引位于向量indices_v的位置pointers_v(j)至pointers_v(j+1)-1中。 矩阵的相应值在矢量数据中位于相同位置。
让我们以示例方式展示如何读取卢瑟福-波音矩阵交换格式Pajek/football中的有趣矩阵。 可以在 www.cise.ufl.edu/research/sp… 的集合中找到具有 118 个非零条目的 35×35 矩阵。
它是参加 1998 年法国举办的 FIFA 世界杯的所有国家橄榄球队的网络的邻接矩阵。网络中的每个节点代表一个国家（或国家橄榄球队），并且链接显示哪个国家向另一个国家输出了球员 国家。
这是football.rb文件的打印输出：
文件的头（前四行）包含重要信息：
- 第一行为我们提供了矩阵标题
Pajek/football; 1998; L. Krempel; ed: V. Batagelj和用于识别目的的数字键MTRXID=1474。 - 第二行包含四个整数值：
TOTCRD=12（在标头之后包含有效数据的行；请参见[24]），PTRCRD=2（包含指针数据的行数），INDCRD=5（包含索引数据的行数） 和VALCRD=2（包含矩阵非零值的行数）。 请注意，它必须为TOTCRD = PTRCRD + INDCRD + VALCRD。 - 第三行表示矩阵类型
MXTYPE=(iua)，在这种情况下，它代表整数矩阵，不对称的压缩列形式。 它还指示行数和列数（NROW=35，NCOL=35），以及非零条目的数量（NNZERO=118）。 对于压缩列格式，不使用最后一个条目，通常将其设置为零。 - 第四列包含以下各列中数据的 Fortran 格式。 指针为
PTRFMT=(20I4)，索引为INDFMT=(26I3)，非零值为VALFMT=(26I3)。
我们继续打开文件进行读取，将头文件后的每一行存储在 Python 列表中，并从文件的相关行中提取填充矢量indptr，indices和data所需的数据 。 我们首先使用data，indices和pointers方法创建对应的 CSR 格式的稀疏矩阵football：
In [22]: f = open("football.rb", 'r'); \
....: football_list = list(f); \
....: f.close()
In [23]: football_data = np.array([])
In [24]: for line in range(4, 4+12):
....: newdata = np.fromstring(football_list[line], sep=" ")
....: football_data = np.append(football_data, newdata)
....:
In [25]: indptr = football_data[:35+1] - 1; \
....: indices = football_data[35+1:35+1+118] - 1; \
....: data = football_data[35+1+118:]
In [26]: football = spsp.csr_matrix((data, indices, indptr),
....: shape=(35,35))
此时，可以借助名为networkx的 Python 模块将网络及其关联的图形可视化。 我们获得下图，描绘了不同国家的节点。 节点之间的每个箭头表示原始国家/地区已将玩家出口到接收国：
注意
networkx是用于处理复杂网络的 Python 模块。 有关更多信息，请访问其在 networkx.github.io 上的 Github 项目页面。
一种完成此任务的方法如下：
In [27]: import networkx
In [28]: G = networkx.DiGraph(football)
In [29]: f = open("football_nodename.txt"); \
....: m = list(f); \
....: f.close()
In [30]: def rename(x): return m[x]
In [31]: G = networkx.relabel_nodes(G, rename)
In [32]: pos = networkx.spring_layout(G)
In [33]: networkx.draw_networkx(G, pos, alpha=0.2, node_color='w',
....: edge_color='b')
模块scipy.sparse从 NumPy 借鉴了一些有趣的概念来创建构造函数和特殊矩阵：
建设者
|
描述
|
| --- | --- |
| scipy.sparse.diags(diagonals, offsets) | 对角线的稀疏矩阵 |
| scipy.sparse.rand(m, n, density) | 指定密度的随机稀疏矩阵 |
| scipy.sparse.eye(m) | 主对角线中带有 1 的稀疏矩阵 |
| scipy.sparse.identity(n) | n × n大小的身份稀疏矩阵 |
函数diags和rand都应通过示例说明其语法。 我们将从具有两个对角线的大小为 14×14 的稀疏矩阵开始：主对角线包含 1s，下面的对角线包含 2s。 我们还创建了一个具有函数scipy.sparse.rand的随机矩阵。 此矩阵的大小为 5×5，具有 25％的非零元素（density=0.25），并且采用 LIL 格式制作：
In [34]: diagonals = [[1]*14, [2]*13]
In [35]: print spsp.diags(diagonals, [0,-1]).todense()
[[ 1\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0.]
[ 2\. 1\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0.]
[ 0\. 2\. 1\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0.]
[ 0\. 0\. 2\. 1\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0.]
[ 0\. 0\. 0\. 2\. 1\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0.]
[ 0\. 0\. 0\. 0\. 2\. 1\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0.]
[ 0\. 0\. 0\. 0\. 0\. 2\. 1\. 0\. 0\. 0\. 0\. 0\. 0\. 0.]
[ 0\. 0\. 0\. 0\. 0\. 0\. 2\. 1\. 0\. 0\. 0\. 0\. 0\. 0.]
[ 0\. 0\. 0\. 0\. 0\. 0\. 0\. 2\. 1\. 0\. 0\. 0\. 0\. 0.]
[ 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 2\. 1\. 0\. 0\. 0\. 0.]
[ 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 2\. 1\. 0\. 0\. 0.]
[ 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 2\. 1\. 0\. 0.]
[ 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 2\. 1\. 0.]
[ 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 0\. 2\. 1.]]
In [36]: S_25_lil = spsp.rand(5, 5, density=0.25, format='lil')
In [37]: S_25_lil
Out[37]:
<5x5 sparse matrix of type '<type 'numpy.float64'>'
with 6 stored elements in LInked List format>
In [38]: print S_25_lil
(0, 0) 0.186663044982
(1, 0) 0.127636181284
(1, 4) 0.918284870518
(3, 2) 0.458768884701
(3, 3) 0.533573291684
(4, 3) 0.908751420065
In [39]: print S_25_lil.todense()
[[ 0.18666304 0\. 0\. 0\. 0\. ]
[ 0.12763618 0\. 0\. 0\. 0.91828487]
[ 0\. 0\. 0\. 0\. 0\. ]
[ 0\. 0\. 0.45876888 0.53357329 0\. ]
[ 0\. 0\. 0\. 0.90875142 0\. ]]
与我们组合ndarray实例的方式类似，我们有一些巧妙的方式来组合稀疏矩阵以构造更复杂的对象：
建设者
|
描述
|
| --- | --- |
| scipy.sparse.bmat(blocks) | 稀疏子块的稀疏矩阵 |
| scipy.sparse.hstack(blocks) | 水平堆叠稀疏矩阵 |
| scipy.sparse.vstack(blocks) | 垂直堆叠稀疏矩阵 |
线性算子
线性运算符基本上是一个函数，该函数通过将输入与矩阵进行左乘法运算，将一列向量作为输入并输出另一个列向量。 尽管从技术上讲，我们可以仅通过处理相应的矩阵来表示这些对象，但还有更好的方法可以做到这一点。
|建设者
|
描述
|
| --- | --- |
| scipy.sparse.linalg.LinearOperator(shape, matvec) | 执行矩阵向量乘积的通用接口 |
| scipy.sparse.linalg.aslinearoperator(A) | 将A返回为LinearOperator |
在scipy.sparse.linalg模块中，我们有一个处理这些对象的公共接口：LinearOperator类。 此类仅具有以下两个属性和三个方法：
shape：表示矩阵的形状dtype：矩阵的数据类型matvec：执行矩阵与向量的乘法rmatvec：通过矩阵与向量的共轭转置来进行乘法matmat：执行一个矩阵与另一个矩阵的乘法
通过示例可以最好地说明其用法。 考虑两个函数，它们分别使用大小为 3 的向量和大小为 4 的输出向量，并分别与两个各自大小为 4×3 的矩阵相乘。我们可以很好地用 lambda 谓词定义这些函数：
In [40]: H1 = np.matrix("1,3,5; 2,4,6; 6,4,2; 5,3,1"); \
....: H2 = np.matrix("1,2,3; 1,3,2; 2,1,3; 2,3,1")
In [41]: L1 = lambda x: H1.dot(x); \
....: L2 = lambda x: H2.dot(x)
In [42]: print L1(np.ones(3))
[[ 9\. 12\. 12\. 9.]]
In [43]: print L2(np.tri(3,3))
[[ 6\. 5\. 3.]
[ 6\. 5\. 2.]
[ 6\. 4\. 3.]
[ 6\. 4\. 1.]]
现在，当我们尝试对这两个函数进行加/减或将它们中的任何一个乘以标量时，就会出现一个问题。 从技术上讲，它应该像加/减相应的矩阵，或将它们乘以任意数字，然后再次执行所需的左乘法一样简单。 但事实并非如此。
例如，我们想写(L1+L2)(v)而不是L1(v) + L2(v)。 不幸的是，这样做会引发错误：
TypeError: unsupported operand type(s) for +: 'function' and
'function'
相反，我们可以实例化相应的线性运算符并按如下方式随意对其进行操作：
In [44]: Lo1 = spspla.aslinearoperator(H1); \
....: Lo2 = spspla.aslinearoperator(H2)
In [45]: Lo1 - 6 * Lo2
Out[45]: <4x3 _SumLinearOperator with dtype=float64>
In [46]: print Lo1 * np.ones(3)
[ 9\. 12\. 12\. 9.]
In [47]: print (Lo1-6*Lo2) * np.tri(3,3)
[[-27\. -22\. -13.]
[-24\. -20\. -6.]
[-24\. -18\. -16.]
[-27\. -20\. -5.]]
当描述具有相关矩阵的产品所需的信息量少于存储矩阵的非零元素所需的内存量时，线性运算符是一个很大的优势。
例如，置换矩阵是平方二进制矩阵（一和零），在每一行和每一列中只有一个条目。 考虑一个由大小为 512×512 的四个块组成的大排列矩阵，例如 1024×1024：一个零块，水平跟随一个标识块，在一个标识块的顶部，水平跟随一个另一个零块。 我们可以通过三种不同的方式存储此矩阵：
In [47]: P_sparse = spsp.diags([[1]*512, [1]*512], [512,-512], \
....: dtype=int)
In [48]: P_dense = P_sparse.todense()
In [49]: mv = lambda v: np.roll(v, len(v)/2)
In [50]: P_lo = spspla.LinearOperator((1024,1024), matvec=mv, \
....: matmat=mv, dtype=int)
在稀疏情况下P_sparse，我们可以将其视为仅存储 1024 个整数。 在密集情况下P_dense，我们在技术上存储 1048576 整数值。 对于线性算子，实际上看起来我们没有存储任何东西！ 指示如何执行乘法的函数mv的占用空间比任何相关矩阵小得多。 这也反映在执行与这些对象的乘法的时间上：
In [51]: %timeit P_sparse * np.ones(1024)
10000 loops, best of 3: 29.7 µs per loop
In [52]: %timeit P_dense.dot(np.ones(1024))
100 loops, best of 3: 6.07 ms per loop
In [53]: %timeit P_lo * np.ones(1024)
10000 loops, best of 3: 25.4 µs per loop
基本矩阵处理
本章第二部分的重点是掌握以下操作：
- 标量乘法，矩阵加法和矩阵乘法
- 痕迹和决定因素
- 转置和反转
- 规范和条件编号
标量乘法，矩阵加法和矩阵乘法
让我们从ndarray类存储的矩阵开始。 我们使用*运算符完成标量乘法，并使用+运算符完成矩阵加法。 但是对于矩阵乘法，我们将需要实例方法dot()或numpy.dot函数，因为运算符*保留用于逐元素乘法：
In [54]: 2*A
Out[54]:
array([[ 2, 4],
[ 8, 32]])
In [55]: A + 2*A
Out[55]:
array([[ 3, 6],
[12, 48]])
In [56]: A.dot(2*A) In [56]: np.dot(A, 2*A)
Out[56]: Out[56]:
array([[ 18, 68], array([[ 18, 68],
[136, 528]]) [136, 528]])
In [57]: A.dot(B)
ValueError: objects are not aligned
In [58]: B.dot(A) In [58]: np.dot(B, A)
Out[58]: Out[58]:
array([[ -9, -34], array([[ -9, -34],
[ 0, 0], [ 0, 0],
[ 9, 34]]) [ 9, 34]])
矩阵类使矩阵乘法更加直观：可以使用运算符*代替dot()方法。 还要注意不同实例类ndarray与矩阵之间的矩阵乘法如何始终强制转换为matrix实例类：
In [59]: C * B
ValueError: shapes (2,2) and (3,2) not aligned: 2 (dim 1) != 3 (dim 0)
In [60]: B * C
Out[60]:
matrix([[ -9, -34],
[ 0, 0],
[ 9, 34]])
对于稀疏矩阵，即使两个稀疏类都不相同，标量乘法和加法也可以与明显的运算符一起很好地工作。 注意每个操作后的结果类转换：
In [61]: S_10_coo = spsp.rand(5, 5, density=0.1, format='coo')
In [62]: S_25_lil + S_10_coo
Out[62]: <5x5 sparse matrix of type '<type 'numpy.float64'>'
with 8 stored elements in Compressed Sparse Row format>
In [63]: S_25_lil * S_10_coo
Out[63]: <5x5 sparse matrix of type '<type 'numpy.float64'>'
with 4 stored elements in Compressed Sparse Row format>
提示
numpy.dot不适用于稀疏矩阵与泛型的矩阵乘法。 我们必须改用运算符*。
In [64]: S_100_coo = spsp.rand(2, 2, density=1, format='coo')
In [65]: np.dot(A, S_100_coo)
Out[66]:
array([[ <2x2 sparse matrix of type '<type 'numpy.float64'>'
with 4 stored elements in COOrdinate format>,
<2x2 sparse matrix of type '<type 'numpy.float64'>'
with 4 stored elements in COOrdinate format>],
[ <2x2 sparse matrix of type '<type 'numpy.float64'>'
with 4 stored elements in COOrdinate format>,
<2x2 sparse matrix of type '<type 'numpy.float64'>'
with 4 stored elements in COOrdinate format>]], dtype=object)
In [67]: A * S_100_coo
Out[68]:
array([[ 1.81 , 1.555],
[ 11.438, 11.105]])
痕迹和决定因素
矩阵的迹线是对角线上元素的总和（假设两个维度上的索引总是增加）。 对于通用矩阵，我们使用实例方法trace()或函数numpy.trace进行计算：
In [69]: A.trace() In [71]: C.trace()
Out[69]: 17 Out[71]: matrix([[17]])
In [70]: B.trace() In [72]: np.trace(B, offset=-1)
Out[70]: -1 Out[72]: 2
为了计算通用平方矩阵的行列式，我们需要在模块scipy.linalg中使用函数det：
In [73]: spla.det(C)
Out[73]: 8.0
转置和反转
可以使用两种实例方法transpose()或T中的任何一种来计算转置，这是两类通用矩阵中的任何一种：
In [74]: B.transpose() In [75]: C.T
Out[74]: Out[75]:
array([[-1, 0, 1], matrix([[ 1, 4],
[-2, 0, 2]]) [ 2, 16]])
可以使用实例方法H为matrix类计算厄米转置：
In [76]: D = C * np.diag((1j,4)); print D In [77]: print D.H
[[ 0.+1.j 8.+0.j] [[ 0.-1.j 0.-4.j]
[ 0.+4.j 64.+0.j]] [ 8.-0.j 64.-0.j]]
使用模块scipy.linalg中的函数inv为ndarray类计算非奇异平方矩阵的逆。 对于matrix类，我们也可以使用实例方法I。 对于非奇异方形稀疏矩阵，我们可以在模块scipy.sparse.linalg中使用函数inv。
提示
稀疏矩阵的逆很少是稀疏的。 因此，不建议通过scipy.sparse.inv功能执行此操作。 解决此问题的一种可能方法是使用todense()实例方法将矩阵转换为通用矩阵，而改用scipy.linear.inv。
但是，由于难以对大矩阵求逆，因此通常最好对逆矩阵进行近似计算。 模块scipy.sparse.linalg中的函数spilu为我们提供了一种非常快速的算法，以CSC格式对平方稀疏矩阵执行此计算。 该算法基于 LU 分解，并在内部进行编码，作为库SuperLU中函数的包装。 它的使用相当复杂，我们将推迟其研究，直到我们探讨matrix分解。
In [78]: E = spsp.rand(512, 512, density=1).todense()
In [79]: S_100_csc = spsp.rand(512, 512, density=1, format='csc')
In [80]: %timeit E.I
10 loops, best of 3: 28.7 ms per loop
In [81]: %timeit spspla.inv(S_100_csc)
1 loops, best of 3: 1.99 s per loop
注意
在执行稀疏逆时，如果输入矩阵不是 CSC 或 CSR 格式，则会收到警告：
/scipy/sparse/linalg/dsolve/linsolve.py:88: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format
warn('spsolve requires A be CSC or CSR matrix format', SparseEfficiencyWarning)
/scipy/sparse/linalg/dsolve/linsolve.py:103: SparseEfficiencyWarning: solve requires b be CSC or CSR matrix format
可以使用模块scipy.linalg中的pinv或pinv2例程，为任何类型的矩阵（不一定为正方形）计算 Moore-Penrose 伪逆。 第一种方法pinv求助于解决最小二乘问题以计算伪逆。 函数pinv2通过基于奇异值分解的方法来计算伪逆。 对于 Hermitian 矩阵或没有复数系数对称的矩阵，我们还有一个名为pinvh的第三函数，它基于特征值分解。
已知在正方形非奇异矩阵的情况下，逆和伪逆是相同的。 这个简单的例子显示了使用以下五种方法计算大型通用对称矩阵的逆的时间：
In [82]: F = E + E.T # F is symmetric
In [83]: %timeit F.I
1 loops, best of 3: 24 ms per loop
In [84]: %timeit spla.inv(F)
10 loops, best of 3: 28 ms per loop
In [85]: %timeit spla.pinvh(E)
1 loops, best of 3: 120 ms per loop
In [86]: %timeit spla.pinv2(E)
1 loops, best of 3: 252 ms per loop
In [87]: %timeit spla.pinv(F)
1 loops, best of 3: 2.21 s per loop
规范和条件编号
对于通用矩阵，我们在scipy.linalg中有七个不同的标准规范。 我们可以在下表中总结它们：
建设者
|
描述
|
| --- | --- |
| norm(A,numpy.inf) | 每行中条目绝对值的总和。 选择最大的价值。 |
| norm(A,-numpy.inf) | 每行中条目绝对值的总和。 选择最小值。 |
| norm(A,1) | 每列中条目的绝对值总和。 选择最大的价值。 |
| norm(A,-1) | 每列中条目的绝对值总和。 选择最小值。 |
| norm(A,2) | 矩阵的最大特征值。 |
| norm(A,-2) | 矩阵的最小特征值。 |
| norm(A,'fro')或norm(A,'f') | Frobenius 范数：乘积A.H * A的迹线的平方根。 |
In [88]: [spla.norm(A,s) for s in (np.inf,-np.inf,-1,1,-2,2,'fro')]
Out[88]: [20, 3, 5, 18, 0.48087417361008861, 16.636368595013604, 16.643316977093239]
提示
对于稀疏矩阵，我们始终可以通过在计算之前应用todense()实例方法来计算范数。 但是，当矩阵的大小太大时，这是不切实际的。 在这些情况下，由于模块scipy.sparse.linalg中的函数onenormest，我们对于 1-范数可以获得的最佳下界是：
In [89]: spla.norm(S_100_csc.todense(), 1) - \
....: spspla.onenormest(S_100_csc)
Out[89]: 0.0
对于 2 范数，我们可能会找到最小和最大特征值的值，但仅适用于平方矩阵。 我们在模块scipy.sparse.linalg中有两种算法可以执行此任务：eigs（用于通用平方矩阵）和 eigsh（用于实对称矩阵）。 在下一节中讨论矩阵分解和因式分解时，我们将详细探讨它们。
注意 SciPy 和 NumPy 的范数计算之间的细微差异。 例如，在 Frobenius 范数的情况下，scipy.linalg.norm直接基于称为 NRM2 的 BLAS 函数，而numpy.linalg.norm等效于形式为sqrt(add.reduce((x.conj() * x).real))的纯粹直接的计算。 当单精度算术中某些数据太大或太小时，基于 BLAS 的代码的优点是速度更快，而且很明显。 在下面的示例中显示：
In [89]: a = np.float64([1e20]); \
....: b = np.float32([1e20])
In [90]: [np.linalg.norm(a), spla.norm(a)]
Out[90]: [1e+20, 1e+20]
In [91]: np.linalg.norm(b)
[...]/numpy/linalg/linalg.py:2056: RuntimeWarning: overflow encountered in multiply
return sqrt(add.reduce((x.conj() * x).real, axis=None))
Out[91]: inf
In [92]: spla.norm(b)
Out[92]: 1.0000000200408773e+20
这不可避免地使我们对非奇异方阵A的条件数的计算进行了讨论。 当我们对输入自变量b进行较小的更改时，此值测量线性方程A * x = b的解的输出将更改多少。 如果该值接近于 1，我们可以放心，解决方案的变化将很小（我们说系统状况良好）。 如果条件数很大，我们就知道系统的计算解可能有问题（然后我们说它是病态的）。
通过将A的范数与其反范数相乘来执行此条件数的计算。 请注意，有不同的条件编号，具体取决于我们为计算选择的标准。 尽管我们需要意识到其明显的局限性，但也可以使用函数numpy.linalg.cond为每个预定义的规范计算这些值。
In [93]: np.linalg.cond(C, -np.inf)
Out[93]: 1.875
矩阵函数
矩阵函数是通过幂级数将平方矩阵映射到另一个平方矩阵的函数。 这些不应与向量化混淆：将一个变量的任何给定函数应用于矩阵的每个元素。 例如，计算方阵A.dot(A)（例如In [8]）的平方与将A的所有元素都平方的矩阵不同（示例In [5]至In [] ）。
注意
为了在符号上进行适当的区分，我们将写A^2表示方阵的实际平方，写A^n表示随后的幂（对于所有正整数n）。
建设者
|
描述
|
| --- | --- |
| scipy.linalg.funm(A, func, disp) | 将名为func的标量值函数扩展到矩阵 |
| scipy.linalg.fractional_matrix_power(A, t) | 小数矩阵幂 |
| scipy.linalg.expm(A) or scipy.sparse.linalg.expm(A) | 矩阵指数 |
| scipy.sparse.linalg.expm_multiply(A,B) | A的矩阵指数对B的作用 |
| scipy.linalg.expm_frechet(A, E) | E方向上矩阵指数的 Frechet 导数 |
| scipy.linalg.cosm(A) | 矩阵余弦 |
| scipy.linalg.sinm(A) | 矩阵 |
| scipy.linalg.tanm(A) | 矩阵切线 |
| scipy.linalg.coshm(A) | 双曲矩阵余弦 |
| scipy.linalg.sinhm(A) | 双曲矩阵正弦 |
| scipy.linalg.tanhm(A) | 双曲矩阵切线 |
| scipy.linalg.signm(A) | 矩阵符号功能 |
| scipy.linalg.sqrtm(A, disp, blocksize) | 矩阵平方根 |
| scipy.linalg.logm(A, disp) | 矩阵对数 |
In [1]: import numpy as np, scipy as sp; \
...: import scipy.linalg as spla
In [2]: np.set_printoptions(suppress=True, precision=3)
In [3]: def square(x): return x**2
In [4]: A = spla.hilbert(4); print A
[[ 1\. 0.5 0.333 0.25 ]
[ 0.5 0.333 0.25 0.2 ]
[ 0.333 0.25 0.2 0.167]
[ 0.25 0.2 0.167 0.143]]
In [5]: print square(A)
[[ 1\. 0.25 0.111 0.062]
[ 0.5 0.333 0.25 0.2 ]
[ 0.333 0.25 0.2 0.167]
[ 0.25 0.2 0.167 0.143]]
In [6]: print A*A
[[ 1\. 0.25 0.111 0.062]
[ 0.25 0.111 0.062 0.04 ]
[ 0.111 0.062 0.04 0.028]
[ 0.062 0.04 0.028 0.02 ]]
In [7]: print A**2
[[ 1\. 0.25 0.111 0.062]
[ 0.25 0.111 0.062 0.04 ]
[ 0.111 0.062 0.04 0.028]
[ 0.062 0.04 0.028 0.02 ]]
In [8]: print A.dot(A)
[[ 1.424 0.8 0.567 0.441]
[ 0.8 0.464 0.333 0.262]
[ 0.567 0.333 0.241 0.19 ]
[ 0.441 0.262 0.19 0.151]]
矩阵的实际功率A^n是定义任何矩阵函数的起点。 在模块numpy.linalg中，我们具有例程matrix_power来执行此操作。 我们还可以使用通用功能funm或功能fractional_matrix_power来实现此结果，它们都在模块scipy.linalg中。
In [9]: print np.linalg.matrix_power(A, 2)
[[ 1.424 0.8 0.567 0.441]
[ 0.8 0.464 0.333 0.262]
[ 0.567 0.333 0.241 0.19 ]
[ 0.441 0.262 0.19 0.151]]
In [10]: print spla.fractional_matrix_power(A, 2)
[[ 1.424 0.8 0.567 0.441]
[ 0.8 0.464 0.333 0.262]
[ 0.567 0.333 0.241 0.19 ]
[ 0.441 0.262 0.19 0.151]]
In [11]: print spla.funm(A, square)
[[ 1.424 0.8 0.567 0.441]
[ 0.8 0.464 0.333 0.262]
[ 0.567 0.333 0.241 0.19 ]
[ 0.441 0.262 0.19 0.151]]
为了计算任何矩阵函数，理论上，我们首先通过泰勒展开将函数表示为幂级数。 然后，我们将输入矩阵应用于该扩展的近似值（因为不可能无限地添加矩阵）。 因此，大多数矩阵函数必然会带来计算错误。 在scipy.linalg模块中，矩阵功能按照该原理编码。
- 注意，三个函数带有可选的布尔参数
disp。 要了解此参数的用法，我们必须记住，大多数矩阵函数都计算近似值，但存在计算误差。 默认情况下，参数disp设置为True，如果近似误差较大，则会发出警告。 如果将disp设置为False，而不是警告，我们将获得估计误差的 1-范数。 - 函数
expm，矩阵上的指数作用expm_multiply和指数expm_frechet的 Frechet 导数背后的算法使用 Pade 近似代替泰勒展开。 这允许更健壮和准确的计算。 所有三角函数和双曲线三角函数均基于涉及expm的简单计算建立其算法。 - 称为
funm的通用矩阵函数和称为sqrtm的平方根函数应用了巧妙的算法，可与输入矩阵的 Schur 分解配合使用，并具有相应的特征值进行适当的代数运算。 它们仍然容易出现舍入错误，但是比任何基于泰勒展开的算法都更快，更准确。 - 称为
signm的矩阵符号函数最初是具有适当功能的funm的应用程序，但是如果该方法失败，则该算法将基于收敛到该解的近似值的迭代采用另一种方法。 - 函数
logm和fractional_matrix_power（将后者应用于非整数幂时）使用 Pade 近似和 Schur 分解的非常复杂的组合（并且进行了改进！）。
提示
当我们处理与特征值有关的矩阵分解时，我们将探讨 Schur 分解。 同时，如果您有兴趣学习这些聪明算法的细节，请阅读 Golub 和 Van Loan 的描述，《矩阵计算 4 版》，约翰霍普金斯大学数学科学学院，第一卷。 3。
有关 Schur-Pade 算法以及指数的 Frechet 导数背后的算法的改进的详细信息，请参阅：
- Nicholas J. Higham 和 Lijing Lin 的一种改进的 Schur-Pade 算法，用于求解矩阵及其分数函数的分数幂
- Awad H. Al-Mohy 和 Nicholas J. Higham 用于矩阵对数的改进的逆缩放和平方算法，《SIAM 科学计算杂志》，34（4）
与求解矩阵方程有关的矩阵分解
矩阵分解的概念使数值线性代数成为科学计算中的有效工具。 如果表示问题的矩阵足够简单，则任何基本的通用算法都可以最佳地找到解决方案（即，快速，最少的数据存储并且没有明显的舍入误差）。 但是，在现实生活中，这种情况很少发生。 在一般情况下，我们要做的是找到合适的矩阵分解并定制在每个因子上均最佳的算法，从而在每个步骤上都获得明显的优势。 在本节中，我们将探讨模块scipy.linalg和scipy.sparse.linalg中包含的不同分解因子，这些因子有助于我们实现矩阵方程的稳健解。
相关因式分解
在此类别中，我们具有以下分解式：
透视 LU 分解
始终可以对平方矩阵A进行乘积分解A = P ● L ● U作为置换矩阵P （执行A行的置换），下三角矩阵L和上三角矩阵U：
建设者
|
描述
|
| --- | --- |
| scipy.linalg.lu(A) | 透视 LU 分解 |
| scipy.linalg.lu_factor(A) | 透视 LU 分解 |
| scipy.sparse.linalg.splu(A) | 透视 LU 分解 |
| scipy.sparse.linalg.spilu(A) | 枢轴 LU 分解不完全 |
胆固醇分解
对于正方形，对称和正定矩阵A，我们可以将矩阵实现为上三角三角形的乘积A = U^T ● U矩阵U与其转置的矩阵，或者作为下三角矩阵L与其转置的乘积A = L^T ● L。 U或L的所有对角线条目均严格为正数：
建设者
|
描述
|
| --- | --- |
| scipy.linalg.cholesky(A) | 胆固醇分解 |
| scipy.linalg.cholesky_banded(AB) | Hermitian 正定带状矩阵的 Cholesky 分解 |
QR 分解
我们可以将大小为 m×n 的任何矩阵实现为大小为 m×m 的平方正交矩阵Q的乘积A = Q ● R与上三角矩阵[ 与A大小相同的R。
建设者
|
描述
|
| --- | --- |
| scipy.linalg.qr(A) | 矩阵的 QR 分解 |
奇异值分解
我们可以实现任何矩阵A作为乘积A = U ● D ● V^H带有对角矩阵D的 unit 矩阵U（对角线中的所有条目均为正数），以及另一个 unit 矩阵V的 Hermitian 转置。 D的对角线上的值称为A的奇异值。
建设者
|
描述
|
| --- | --- |
| scipy.linalg.svd(A) | 奇异值分解 |
| scipy.linalg.svdvals(A) | 奇异值 |
| scipy.linalg.diagsvd(s, m, n) | SVD 的对角矩阵，由奇异值 s 和指定大小 |
| scipy.sparse.linalg.svds(A) | 稀疏矩阵的最大k奇异值/向量 |
矩阵方程
在 SciPy 中，我们具有基于以下情况的强大算法来求解任何矩阵方程：
- 给定一个正方形矩阵
A和一个右侧b（它可以是一维矢量或与A相同行数的另一个矩阵），基本系统如下：A ● x = bA^T ● x = bA^H ● x = b
- 给定任何矩阵
A（不一定是正方形）和适当大小的右侧矢量/矩阵b，方程A ● x = b的最小二乘解。 也就是说，找到了使表达A ● x - b的 Frobenius 范数最小的向量x。 - 对于与之前相同的情况，以及一个额外的阻尼系数
d，方程A ● x = b的正则化最小二乘解使函数norm(A * x - b, 'f')**2 + d**2 * norm(x, 'f')**2最小。 - 给定方阵
A和B，以及具有适当大小的右侧矩阵Q，Sylvester 系统为A ● X + X ● B = Q。 - 对于适当大小的方阵
A和矩阵Q，连续的 Lyapunov 方程为A ● X + X ● A^H = Q。 - 对于矩阵
A和Q，与之前的情况一样，离散 Lyapunov 方程为X - A ● X ● A^H = Q。 - 给定平方矩阵
A，Q和R，以及具有适当大小的另一个矩阵B，连续代数 Riccati 方程为A^T ● X + X ● A - X ● B ● R^(-1) ● B^T ● X + Q = 0。 - 对于与前述情况相同的矩阵，离散代数 Riccati 方程为
X = A^T ● X ● A - (A^T ● X ● B) ● (R + B^T ● X ● B)^(-1) ● (B^T ● X ● A) + Q。
无论如何，用 SciPy 掌握矩阵方程式基本上意味着确定所涉及的矩阵并在库中选择最合适的算法来执行所请求的运算。 除了能够以最小的舍入误差量来计算解决方案外，我们还需要以尽可能最快的方式并使用尽可能少的内存资源来进行计算。
向前和向后替代
让我们从最简单的情况开始：线性方程组的基本系统A ● x = b（或其他两个变体），其中A是通用的下限或上限 三角方阵。 从理论上讲，这些系统可以通过正向替换（对于较低的三角矩阵）或向后替换（对于较高的三角矩阵）轻松解决。 在 SciPy 中，我们使用模块scipy.linalg中的功能solve_triangular完成此任务。
在此初始示例中，我们将A构造为大小为 1024×1024 的下三角 Pascal 矩阵，其中已过滤了非零值：奇数变为 1，而偶数变为 0。 右侧b是带有1024的矢量。
In [1]: import numpy as np, \
...: scipy.linalg as spla, scipy.sparse as spsp, \
...: scipy.sparse.linalg as spspla
In [2]: A = (spla.pascal(1024, kind='lower')%2 != 0)
In [3]: %timeit spla.solve_triangular(A, np.ones(1024))
10 loops, best of 3: 6.64 ms per loop
为了解决涉及矩阵A的其他相关系统，我们采用了可选参数trans（默认设置为0或N，为基本系统A ● x = b）。 如果将trans设置为T或1，我们将求解系统A^T ● x = b。 如果将trans设置为C或2，则我们求解A^H ● x = b。
注意
函数solve_triangular是LAPACK函数trtrs的包装器。
基本系统：带状矩阵
就算法简单性而言，接下来的情况是基本系统A ● x = b，其中A是方带矩阵。 我们使用例程solve_banded（用于通用带状矩阵）或solveh_banded（用于复杂的 Hermitian 带状矩阵的通用实对称）。 它们都属于模块scipy.linalg。
注意
函数solve_banded和solveh_banded分别是LAPACK函数GBSV和PBSV的包装。
这两个函数都不接受通常格式的矩阵。 例如，由于solveh_banded需要一个对称的带状矩阵，因此该函数仅需要输入从上到下依次存储的主对角线上，下和上的对角线元素。
通过一个具体的例子可以最好地解释这种输入方法。 取以下对称带状矩阵：
2 -1 0 0 0 0
-1 2 -1 0 0 0
0 -1 2 -1 0 0
0 0 -1 2 -1 0
0 0 0 -1 2 -1
0 0 0 0 -1 2
矩阵的大小为 6×6，只有三个非零对角线，由于对称性，其中两个是相同的。 我们通过以下两种方式之一收集大小为 2×6 的ndarray中两个相关的非零对角线：
-
如果我们决定从上三角矩阵输入条目，则首先收集从上到下的对角线（在主对角线结束），右对齐：
* -1 -1 -1 -1 -1 2 2 2 2 2 2 -
如果我们决定从较低的三角矩阵输入条目，则从顶部到底部（从主对角线开始）收集对角线，左对齐：
2 2 2 2 2 2 -1 -1 -1 -1 -1 * In [4]: B_banded = np.zeros((2,6)); \ ...: B_banded[0,1:] = -1; \ ...: B_banded[1,:] = 2 In [5]: spla.solveh_banded(B_banded, np.ones(6)) Out[5]: array([ 3., 5., 6., 6., 5., 3.])
对于非对称带状方矩阵，我们改用solve_banded，并且输入矩阵也需要以这种特殊方式存储：
- 计算主对角线下的非零对角线数（将其设置为
l）。 计算主对角线上非零对角线的数量（将其设置为u）。 设置r = l + u + 1。 - 如果矩阵的大小为
n × n，则用n列和r行创建ndarray。 我们简称为AB形式的矩阵或AB矩阵。 - 从上到下，仅按顺序将相关的非零对角线存储在 AB 矩阵中。 主对角线上的对角线是右对齐的； 主对角线下方的对角线保持左对齐。
让我们用另一个例子来说明这个过程。 我们输入以下矩阵：
2 -1 0 0 0 0
-1 2 -1 0 0 0
3 -1 2 -1 0 0
0 3 -1 2 -1 0
0 0 3 -1 2 -1
0 0 0 3 -1 2
In [6]: C_banded = np.zeros((4,6)); \
...: C_banded[0,1:] = -1; \
...: C_banded[1,:] = 2; \
...: C_banded[2,:-1] = -1; \
...: C_banded[3,:-2] = 3; \
...: print C_banded
[[ 0\. -1\. -1\. -1\. -1\. -1.]
[ 2\. 2\. 2\. 2\. 2\. 2.]
[-1\. -1\. -1\. -1\. -1\. 0.]
[ 3\. 3\. 3\. 3\. 0\. 0.]]
要调用求解器，我们需要手动输入对角线上方和下方的对角线数量，以及AB矩阵和系统的右侧：
In [7]: spla.solve_banded((2,1), C_banded, np.ones(6))
Out[7]:
array([ 0.86842105, 0.73684211, -0.39473684, 0.07894737,
1.76315789, 1.26315789])
让我们检查一下可以包含在这两个函数的调用中的可选参数：
|参数
|
默认值
|
描述
|
| --- | --- | --- |
| l_and_u | (int, int) | 非零上下对角线数 |
| ab | AB格式的矩阵 | 带状方矩阵 |
| b | ndarray | 右侧 |
| overwrite_ab | 布尔型 | 舍弃ab中的数据 |
| overwrite_b | 布尔型 | 舍弃b中的数据 |
| check_finite | 布尔型 | 是否检查输入矩阵是否包含有限数 |
提示
scipy.linalg模块中需要矩阵作为输入和输出的所有函数，或者方程组的解或因式分解的函数，都有两个我们需要熟悉的可选参数：overwrite_x（对于矩阵中的每个矩阵/矢量 输入）和check_finite。 它们都是布尔值。
默认情况下，overwrite选项设置为False。 如果我们不关心保留输入矩阵的值，则可以使用内存中的同一对象执行操作，而不是在内存中创建另一个具有相同大小的对象。 在这种情况下，我们可以提高速度并减少资源消耗。
默认情况下，check_finite选项设置为True。 在存在数据的算法中，存在可选的数据完整性检查。 如果在任何给定时刻，任何值为(+/-)numpy.inf或NaN，则过程将暂停，并引发异常。 我们可能会关闭此选项，从而导致更快的解决方案，但是如果数据在计算中的任何时候被破坏，代码可能会崩溃。
函数solveh_banded具有一个额外的可选布尔参数lower，该参数最初设置为False。 如果设置为True，则必须输入目标 AB 矩阵的下三角矩阵，而不是上三角矩阵（输入约定与以前相同）。
基本系统：通用平方矩阵
对于A是通用方矩阵的基本系统的解决方案，将A分解为一个好主意，以使某些（或全部）因子为三角形，然后在适当时应用来回替换。 这是枢轴LU和 Cholesky 分解背后的想法。
如果矩阵A是实对称（或复 Hermitian）并且是正定的，则最佳策略是通过应用两个可能的 Cholesky 分解中的任何一个A = U^H ● U或A = L ● L^H，带有U和L上/下三角矩阵。
例如，如果我们使用带有上三角矩阵的形式，则基本方程组A ● x = b的解变成 U H ●U ●x = b 。 设置 y = U ●x 并求解系统 U H ●y = b 对于y 通过正向替换。 现在我们有了一个新的三角系统 U ●x = y ，我们可以通过逆向替换来求解x。
为了用这种技术执行这样的系统的解决方案，我们首先通过使用函数cholesky，cho_factor或cholesky_banded来计算因式分解。 然后将输出用于求解器cho_solve。
对于 Cholesky 分解，称为cholesky，cho_factor和cholesky_banded的三个相关函数具有与solveh_banded相似的一组选项。 他们接受一个较低的附加布尔选项（默认设置为False），该选项决定是输出较低的三角形分解还是较高的三角形分解。 函数cholesky_banded需要AB格式的矩阵作为输入。
现在让我们使用所有三种方法测试矩阵B的 Cholesky 分解：
In [8]: B = spsp.diags([[-1]*5, [2]*6, [-1]*5], [-1,0,1]).todense()
...: print B
[[ 2\. -1\. 0\. 0\. 0\. 0.]
[-1\. 2\. -1\. 0\. 0\. 0.]
[ 0\. -1\. 2\. -1\. 0\. 0.]
[ 0\. 0\. -1\. 2\. -1\. 0.]
[ 0\. 0\. 0\. -1\. 2\. -1.]
[ 0\. 0\. 0\. 0\. -1\. 2.]]
In [9]: np.set_printoptions(suppress=True, precision=3)
In [10]: print spla.cholesky(B)
[[ 1.414 -0.707 0\. 0\. 0\. 0\. ]
[ 0\. 1.225 -0.816 0\. 0\. 0\. ]
[ 0\. 0\. 1.155 -0.866 0\. 0\. ]
[ 0\. 0\. 0\. 1.118 -0.894 0\. ]
[ 0\. 0\. 0\. 0\. 1.095 -0.913]
[ 0\. 0\. 0\. 0\. 0\. 1.08 ]]
In [11]: print spla.cho_factor(B)[0]
[[ 1.414 -0.707 0\. 0\. 0\. 0\. ]
[-1\. 1.225 -0.816 0\. 0\. 0\. ]
[ 0\. -1\. 1.155 -0.866 0\. 0\. ]
[ 0\. 0\. -1\. 1.118 -0.894 0\. ]
[ 0\. 0\. 0\. -1\. 1.095 -0.913]
[ 0\. 0\. 0\. 0\. -1\. 1.08 ]]
In [12]: print spla.cholesky_banded(B_banded)
[[ 0\. -0.707 -0.816 -0.866 -0.894 -0.913]
[ 1.414 1.225 1.155 1.118 1.095 1.08 ]]
cho_factor的输出是一个元组：第二个元素是布尔低位。 第一个元素是ndarray，代表方矩阵。 如果将lower设置为True，则在A的 Cholesky 分解中，此ndarray的下三角子矩阵为L。 如果将lower设置为False，则在A的因式分解中，上三角子矩阵为U。 矩阵中的其余元素是随机的，而不是零，因为cho_solve不使用它们。 以类似的方式，我们可以通过cho_banded的输出调用cho_solve_banded来解决适当的系统。
注意
cholesky和cho_factor都是同一个名为potrf的 LAPACK 函数的包装，具有不同的输出选项。 cholesky_banded调用pbtrf。 cho_solve函数是potrs的包装器，cho_solve_banded调用pbtrs。
然后，我们准备使用两个选项之一来解决系统问题：
In [13]: spla.cho_solve((spla.cholesky(B), False), np.ones(6))
Out[13]: array([ 3., 5., 6., 6., 5., 3.])
In [13]: spla.cho_solve(spla.cho_factor(B), np.ones(6))
Out[13]: array([ 3., 5., 6., 6., 5., 3.])
对于任何其他种类的通用方阵A，求解基本系统的下一个最佳方法 A ●x = b 被枢轴化LU分解。 这等效于找到置换矩阵P以及三角矩阵U（上部）和L（下部），以便 P ●A = L ●U 。 在这种情况下，根据P对系统中的行进行置换可以得到等式（ P ●A） ●x = P ●b 。 设置c = P ● b和y = U ● x，并使用正向替换在系统 L 中求解y。 然后，在系统 U ●x = y 中用逆向替换求解x。
执行此操作的相关功能是模块scipy.linalg中的lu，lu_factor（用于分解）和lu_solve（用于解决方案）。 对于稀疏矩阵，在模块scipy.sparse.linalg中具有splu和spilu。
让我们首先开始进行因式分解实验。 在此示例中，我们使用大循环矩阵（非对称）：
In [14]: D = spla.circulant(np.arange(4096))
In [15]: %timeit spla.lu(D)
1 loops, best of 3: 7.04 s per loop
In [16]: %timeit spla.lu_factor(D)
1 loops, best of 3: 5.48 s per loop
注意
lu_factor函数是 LAPACK 中所有*getrf例程的包装。 lu_solve函数是getrs的包装器。
函数lu具有一个额外的布尔选项：permute_l（默认设置为False）。 如果设置为True，则该功能仅输出两个矩阵 PL = P ●L （正确排列的下三角矩阵）和U。 否则，输出是该顺序的三元组P，L和U。
In [17]: P, L, U = spla.lu(D)
In [17]: PL, U = spla.lu(D, permute_l=True)
函数lu_factor的输出节省资源。 我们获得矩阵LU，其中上三角U和下三角L。 我们还获得了整数dtype，piv的一维ndarray类，表示表示置换矩阵P的枢轴索引。
In [18]: LU, piv = spla.lu_factor(D)
求解器lu_solve将lu_factor的两个输出，一个右侧矩阵b和可选指示器trans带到要解决的基本系统类型：
In [19]: spla.lu_solve(spla.lu_factor(D), np.ones(4096))
Out[19]: array([ 0., 0., 0., ..., 0., 0., 0.])
提示
此时，我们必须对模块scipy.linalg中的通用功能solve进行注释。 它是LAPACK功能POSV和GESV的包装。 它允许我们输入矩阵A和右侧矩阵b，并指示A是否对称且为正定。 无论如何，例程都会在内部决定要使用的两个分解因子（Cholesky 或透视 LU）中的哪一个，并据此计算解决方案。
对于大的稀疏矩阵，如果它们以 CSC 格式存储，则可以使用模块scipy.sparse.linalg中的函数splu或spilu更有效地执行枢轴LU分解。 这两个函数直接使用SuperLU库。 它们的输出不是一组矩阵，而是一个名为scipy.sparse.linalg.dsolve._superlu.SciPyLUType的 Python 对象。 该对象具有四个属性和一个实例方法：
shape：包含矩阵A形状的 2 元组nnz：矩阵A中非零条目的数量perm_c, perm_r：分别应用于矩阵A的列和行的置换，以获得计算出的LU分解solve：将对象转换为接受ndarray b的函数object.solve(b,trans)和可选描述字符串trans的实例方法。
一个大想法是，处理大量数据时，LU分解中的实际矩阵不像分解背后的主要应用（系统的解决方案）那么重要。 所有执行此操作的相关信息都以最佳方式存储在对象的方法solve中。
splu和spilu之间的主要区别在于，后者计算不完全分解。 有了它，我们可以获得矩阵A的逆的非常好的近似值，并使用矩阵乘法来计算大型系统的解，而这只花费了计算实际解的时间的一小部分。
注意
这两个功能的用法相当复杂。 目的是使用对角矩阵Dr和Dc以及置换矩阵Pr和Pc来计算形式为 Pr * Dr * A * Dc * Pc = L * U 的因式分解。 想法是手动平衡基质A，以便使乘积 B = Dr * A * Dc 比A更好。 如果有可能在并行体系结构中解决此问题，我们可以通过最佳地重新排列行和列来提供帮助。 然后，手动输入置换矩阵Pr和Pc以对B的行和列进行预排序。 所有这些选项都可以提供给splu或spilu。
该算法利用放宽超级节点的想法来减少无效的间接寻址和符号时间（允许使用更高级别的 BLAS 操作）。 我们可以选择确定这些对象的程度，以使算法适合手头的矩阵。
有关算法和所有不同选项的完整说明，最好的参考是《 SuperLU 用户指南》，可在 crd-legacy.lbl.gov/~xiaoye/Sup… 在线找到。
让我们用一个简单的示例来说明这一点，其中不需要对行或列进行置换。 在一个较大的下三角 Pascal 矩阵中，所有偶值条目都变为零，所有奇值条目都变为 1。 将此用作矩阵A。 对于右侧，使用一个 1 的向量：
In [20]: A_csc = spsp.csc_matrix(A, dtype=np.float64)
In [21]: invA = spspla.splu(A_csc)
In [22]: %time invA.solve(np.ones(1024))
CPU times: user: 4.32 ms, sys: 105 µs, total: 4.42 ms
Wall time: 4.44 ms
Out[22]: array([ 1., -0., 0., ..., -0., 0., 0.])
In [23]: invA = spspla.spilu(A_csc)
In [24]: %time invA.solve(np.ones(1024))
CPU times: user 656 µs, sys: 22 µs, total: 678 µs
Wall time: 678 µs
Out[24]: array([ 1., 0., 0., ..., 0., 0., 0.])
注意
将在稀疏矩阵上执行过程的时间与本节开始处在相应矩阵A上的初始solve_triangular过程进行比较。 哪个过程更快？
但是，通常，如果必须解决一个基本系统并且矩阵A很大且稀疏，我们更喜欢使用迭代方法快速收敛到实际解。 当它们收敛时，它们对舍入误差始终不那么敏感，因此在计算量非常大时更适合。
在scipy.sparse.linalg模块中，我们有八种不同的迭代方法，所有这些方法都接受以下参数：
- 任意格式的矩阵
A（矩阵，ndarray，稀疏矩阵，甚至线性运算符！），右侧矢量/矩阵b为ndarray。 - 初步猜测为
x0，为ndarray。 - 公差
l，一个浮点数。 如果连续迭代的差小于此值，则代码停止，最后计算出的值作为解决方案输出。 - 允许的最大迭代次数，maxiter，一个整数。
- 预调节器稀疏矩阵
M应当近似A的逆。 - 每次迭代后调用的当前解矢量
xk的callback函数。
建设者
|
描述
|
| --- | --- |
| bicg | 双共轭梯度迭代 |
| bicgstab | 双共轭梯度稳定迭代 |
| cg | 共轭梯度迭代 |
| cgs | 共轭梯度平方迭代 |
| gmres | 广义最小残差迭代 |
| lgmres | LGMRES 迭代 |
| minres | 最小残差迭代 |
| qmr | 拟最小残差迭代 |
选择正确的迭代方法，良好的初步猜测，尤其是成功的 Preconditioner，本身就是一门艺术。 它涉及学习诸如泛函分析中的运算符或 Krylov 子空间方法之类的主题，这些内容远远超出了本书的范围。 在这一点上，为了满足比较，我们满意地展示了一些简单的示例：
In [25]: spspla.cg(A_csc, np.ones(1024), x0=np.zeros(1024))
Out[25]: (array([ nan, nan, nan, ..., nan, nan, nan]), 1)
In [26]: %time spspla.gmres(A_csc, np.ones(1024), x0=np.zeros(1024))
CPU times: user 4.26 ms, sys: 712 µs, total: 4.97 ms
Wall time: 4.45 ms
Out[26]: (array([ 1., 0., 0., ..., -0., -0., 0.]), 0)
In [27]: Nsteps = 1
....: def callbackF(xk):
....: global Nsteps
....: print'{0:4d} {1:3.6f} {2:3.6f}'.format(Nsteps, \
....: xk[0],xk[1])
....: Nsteps += 1
....:
In [28]: print '{0:4s} {1:9s} {1:9s}'.format('Iter', \
....: 'X[0]','X[1]'); \
....: spspla.bicg(A_csc, np.ones(1024), x0=np.zeros(1024),
....: callback=callbackF)
....:
Iter X[0] X[1]
1 0.017342 0.017342
2 0.094680 0.090065
3 0.258063 0.217858
4 0.482973 0.328061
5 0.705223 0.337023
6 0.867614 0.242590
7 0.955244 0.121250
8 0.989338 0.040278
9 0.998409 0.008022
10 0.999888 0.000727
11 1.000000 -0.000000
12 1.000000 -0.000000
13 1.000000 -0.000000
14 1.000000 -0.000000
15 1.000000 -0.000000
16 1.000000 0.000000
17 1.000000 0.000000
Out[28]: (array([ 1., 0., 0., ..., 0., 0., -0.]), 0)
最小二乘
给定通用矩阵A（不一定是正方形）和右侧矢量/矩阵b，我们寻找矢量/矩阵x，使得表达式 A 的 Frobenius 范数 ●x-b 被最小化。
在scipy中考虑了用数字方式解决此问题的三种主要方法：
- 正态方程
- QR 分解
- 奇异值分解
正态方程
正规方程将最小二乘问题简化为使用对称（不一定是正定）矩阵求解线性方程组的基本系统。 它的速度非常快，但是由于存在舍入错误，因此可能不准确。 基本上，它等于求解系统*（A* H ●A） ●x = A H ●b 。 这等效于求解 x =（A H ●A） -1 ●A H ●b = pinv（A） ●b 。
让我们以示例显示：
In [29]: E = D[:512,:256]; b = np.ones(512)
In [30]: sol1 = np.dot(spla.pinv2(E), b)
In [31]: sol2 = spla.solve(np.dot(F.T, F), np.dot(F.T, b))
QR 分解
QR 分解将任何矩阵转换为正交/ unit 矩阵Q与正方形上三角矩阵R的乘积 A = Q ●R 。 这使我们无需对任何矩阵求逆即可求解系统（因为 Q H = Q -1 ），因此 A ●x = b 变成 R ●x = Q H ●b ，可通过反取代轻松解决。 请注意，以下两种方法是等效的，因为模式economic报告了最大秩的子矩阵：
In [32]: Q, R = spla.qr(E); \
....: RR = R[:256, :256]; BB = np.dot(Q.T, b)[:256]; \
....: sol3 = spla.solve_triangular(RR, BB)
In [32]: Q, R = spla.qr(E, mode='economic'); \
....: sol3 = spla.solve_triangular(R, np.dot(Q.T, b))
奇异值分解
正规方程和 QR 因式分解这两种方法都可以快速工作，并且只有在A的等级满时才可靠。 如果不是这种情况，则必须使用奇异值分解 A = U ●D ●V H 以及单一矩阵U和V和对角矩阵D，其中对角线上的所有条目均为正值。 这允许快速解决方案 x = V ●D -1 ●U H ● b 。
请注意，下面讨论的两种方法是等效的，因为设置为False的选项full_matrices报告了可能的最小大小的子矩阵：
In [33]: U, s, Vh = spla.svd(E); \
....: Uh = U.T; \
....: Si = spla.diagsvd(1./s, 256, 256); \
....: V = Vh.T; \
....: sol4 = np.dot(V, Si).dot(np.dot(Uh, b)[:256])
In [33]: U, s, Vh = spla.svd(E, full_matrices=False); \
....: Uh = U.T; \
....: Si = spla.diagsvd(1./s, 256, 256); \
....: V = Vh.T; \
....: sol4 = np.dot(V, Si).dot(np.dot(Uh, b))
模块scipy.linalg具有一个实际上使用 SVD 方法执行最小二乘的功能：lstsq。 无需手动转置，求逆和乘法所有需要的矩阵。 它是LAPACK函数GELSS的包装。 它输出所需的解以及计算的残差，有效秩和输入矩阵A的奇异值。
In [34]: sol5, residue, rank, s = spla.lstsq(E, b)
请注意，我们执行的所有计算如何提供彼此非常接近的解决方案（如果不相等！）：
In [35]: map(lambda x: np.allclose(sol5,x), [sol1, sol2, sol3, sol4])
Out[35]: [True, True, True, True]
正则化最小二乘
在大稀疏矩阵的情况下，模块scipy.sparse.linalg具有两种用于最小二乘的迭代方法，即lsqr和lsmr，这允许使用更通用的阻尼因子d进行迭代。 我们试图使功能norm(A * x - b, 'f')**2 + d**2 * norm(x, 'f')**2最小化。 用法和参数与我们之前研究的迭代函数非常相似。
其他矩阵方程求解器
下表总结了其余的矩阵方程求解器。 这些例程都没有任何参数可用于性能或内存管理或检查数据的完整性：
|建设者
|
描述
|
| --- | --- |
| solve_sylvester(A, B, Q) | 西尔维斯特方程 |
| solve_continuous_are(A, B, Q, R) | 连续代数 Riccati 方程 |
| solve_discrete_are(A, B, Q, R) | 离散代数 Riccati 方程 |
| solve_lyapunov(A, Q) | 连续李雅普诺夫方程 |
| solve_discrete_lyapunov(A, Q) | 离散 Lyapunov 方程 |
基于特征值的矩阵分解
在此类中，我们对平方矩阵有两种因式分解：频谱分解和 Schur 分解（尽管从技术上讲，频谱分解是 Schur 分解的一种特殊情况）。 两者的目的最初是同时显示一个或几个矩阵的特征值，尽管它们的应用程序有很大不同。
光谱分解
我们考虑以下四种情况：
- 给定一个平方矩阵
A，我们求出某些实数或复数值m（相应的特征值）满足 A●v = m●v 的所有向量v（右特征向量）。 如果所有特征向量都不相同，我们将它们收集为矩阵V的列（恰好是可逆的）。 它们对应的特征值以与对角矩阵D的对角条目相同的顺序存储。 然后我们可以将A理解为乘积 A = V●D●V -1 。 我们将此分解称为普通特征值问题。 - 给定一个平方矩阵
A，我们寻找特征值m满足 v●A = m●v 的所有向量v（左特征向量）。 如前所述，如果所有特征向量都不同，则将它们收集在矩阵V中。 它们的对应特征值收集在对角矩阵D中。 然后可以将矩阵A分解为乘积 A = V●D●V -1 。 我们也将此分解称为普通特征值问题。 特征值与前面的情况相同。 - 在给定具有相同大小的平方矩阵
A和B的情况下，我们寻求满足 m●A●v = n●B●v 的所有向量v（广义右特征向量） 复数值m和n。 当比率 r = n / m 是可计算的时，称为广义特征值。 特征向量被收集为矩阵V的列，其对应的广义特征值r被收集在对角矩阵D中。 然后我们可以通过标识 A = B●V●D●V -1 来实现A和B之间的关系。 我们将此身份称为广义特征值问题。 - 对于与之前相同的情况，如果我们求向量
v（广义左特征向量）以及满足 m●v●A = n●v●B 的值m和n另一个类似的分解。 我们再次将这种分解称为广义特征值问题。
模块scipy.linalg和scipy.sparse.linalg中的以下功能可帮助我们计算特征值和特征向量：
建设者
|
描述
|
| --- | --- |
| scipy.linear.eig(A[, B]) | 普通/广义特征值问题 |
| scipy.linalg.eigvals(A[, B]) | 普通/广义特征值问题的特征值 |
| scipy.linalg.eigh(A[, B]) | 普通/广义特征值问题。 厄米/对称矩阵 |
| scipy.linalg.eigvalsh(A[, B]) | 普通/广义特征值问题的特征值； 厄米/对称矩阵 |
| scipy.linalg.eig_banded(AB) | 普通特征值问题； 厄米/对称带矩阵 |
| scipy.linalg.eigvals_banded(AB) | 普通特征值问题的特征值； 厄米/对称带矩阵 |
| scipy.sparse.linalg.eigs(A, k) | 查找 k 个特征值和特征向量 |
| scipy.sparse.linalg.eigsh(A, k) | 查找k特征值和特征向量； 实对称矩阵 |
| scipy.sparse.linalg.lobpcg(A, X) | 具有可选预处理A对称的普通/广义特征值问题 |
对于矩阵不对称或不带状的任何一种特征值问题，我们使用函数eig，该函数是LAPACK例程GEEV和GGEV（后者用于广义特征值问题）的包装器。 对于仅输出特征值而不输出特征向量的情况，函数eigvals是语法糖。 为了报告我们是否需要左右特征向量，我们使用可选的布尔参数left和right。 默认情况下，left设置为False，right设置为True，因此提供了正确的特征向量。
对于具有非带实对称或 Hermitian 矩阵的特征值问题，我们使用函数eigh，该函数是形式为*EVR，*GVD和*GV的 LAPACK 例程的包装。 我们可以选择使用可选参数eigvals输出任意数量的特征值。 这是一个整数元组，指示所需的最低和最高特征值的索引。 如果省略，则返回所有特征值。 在这种情况下，可以使用基于分而治之技术的更快算法进行计算。 我们可以使用可选的布尔参数turbo（默认设置为False）来指示此选择。
如果仅希望报告特征值，则可以将可选参数eigvals_only设置为True，或使用相应的语法糖eighvals。
我们在scipy.linalg模块中考虑的最后一种情况是带状实对称或 Hermitian 矩阵的特征值问题。 我们使用函数eig_banded，确保输入矩阵为 AB 格式。 该函数是LAPACK例程*EVX的包装。
对于非常大的矩阵，特征值的计算通常在计算上是不可能的。 如果这些大型矩阵稀疏，则可以使用两种迭代算法来计算一些特征值，即隐式重启 Arnoldi 和隐式重启 Lanczos 方法（后者用于对称或 Hermitian 方法） 矩阵）。 模块scipy.sparse.linalg具有两个功能eigs和eigsh，它们是执行它们的ARPACK例程*EUPD的包装。 我们还具有执行另一种迭代算法的函数lobpcg，即局部最优块预处理共轭梯度方法。 该函数接受预处理器，因此有可能更快地收敛到所需的特征值。
我们将通过一个有趣的矩阵Andrews来说明所有这些函数的用法。 它创建于 2003 年，旨在对本征值问题的内存高效算法进行基准测试。 它是大小为 60,000×60,000 和 760,154 个非零条目的实对称稀疏矩阵。 可以从的稀疏矩阵集合中下载 www.cise.ufl.edu/research/sparse/matrices/Andrews/Andrews.html 。
对于此示例，我们以 Matrix Market 格式Andrews.mtx下载了该矩阵。 请注意，矩阵是对称的，并且文件仅提供主对角线上或下方的数据。 收集所有这些信息之后，我们确保也填充上三角形：
In [1]: import numpy as np, scipy.sparse as spsp, \
...: scipy.sparse.linalg as spspla
In [2]: np.set_printoptions(suppress=True, precision=6)
In [3]: rows, cols, data = np.loadtxt("Andrews.mtx", skiprows=14,
...: unpack=True); \
...: rows-=1; \
...: cols-=1
In [4]: A = spsp.csc_matrix((data, (rows, cols)), \
...: shape=(60000,60000)); \
...: A = A + spsp.tril(A, k=1).transpose()
我们首先计算绝对值中最大的前五个特征值。 我们使用选项which='LM'调用函数eigsh。
In [5]: %time eigvals, v = spspla.eigsh(A, 5, which='LM')
CPU times: user 3.59 s, sys: 104 ms, total: 3.69 s
Wall time: 3.13 s
In [6]: print eigvals
[ 69.202683 69.645958 70.801108 70.815224 70.830983]
通过切换到选项which='SM'，我们也可以根据绝对值来计算最小的特征值：
In [7]: %time eigvals, v = spspla.eigsh(A, 5, which='SM')
CPU times: user 19.3 s, sys: 532 ms, total: 19.8 s
Wall time: 16.7 s
In [8]: print eigvals
[ 10.565523 10.663114 10.725135 10.752737 10.774503]
提示
ARPACK中的例程在查找较小特征值时效率不高。 在这种情况下，通常最好采用移位反转模式以获得更好的性能。 有关此过程的信息，请阅读 www.caam.rice.edu/software/AR… 中的描述，或 RB Lehoucq，DC Sorensen 和 C.Yang，ARPACK 的文章。 用户指南：通过隐式重新启动的 Arnoldi 方法解决大规模特征值问题。 SIAM，宾夕法尼亚州费城，1998 年。
注意
函数eigsh允许我们通过指示接近所需特征值的值来执行移位反转模式。 如果我们有很好的猜测（如上一步所述），则可以将此过程与选项sigma结合使用，并将策略与选项模式结合使用。 在这种情况下，我们还需要提供线性运算符而不是矩阵。 执行时间要慢得多，但结果通常要精确得多（尽管给出的示例并不建议这样做！）。
In [9]: A = spspla.aslinearoperator(A)
In [10]: %time spspla.eigsh(A, 5, sigma=10.0, mode='cayley')
CPU times: user 2min 5s, sys: 916 ms, total: 2min 6s
Wall time: 2min 6s
In [11]: print eigvals
[ 10.565523 10.663114 10.725135 10.752737 10.774503]
舒尔分解
有四种情况：
- 具有复系数的方阵
A的复数 Schur 分解。 我们可以将A理解为 the 矩阵U与上三角矩阵T的乘积 A = U●T●U H 和U的 Hermitian 转置。 我们称T为A的复杂 Schur 形式。T的对角线中的条目是A的特征值。 - 具有实系数的方阵
A的实数舒尔分解。 如果矩阵的所有特征值均为实值，则可以将矩阵实现为正交矩阵的乘积 A = V●S●V TV和一个上三角矩阵S以及V的转置。S中的块的大小为 1×1 或 2×2。如果块为 1×1，则该值是A的实际特征值之一。 任何 2×2 块代表A的一对复共轭特征值。 我们称S为A的真正 Schur 形式。 - 两个平方矩阵
A和B的复数广义 Schur 分解。 我们可以同时将它们分解为以下形式： A = Q●S●Z H 和 B = Q●T●Z H 具有相同的单一矩阵Q和Z。 矩阵S和T均为上三角，它们对角线元素的比率恰好是A和B的广义特征值。 - 两个实值平方矩阵
A和B的实数广义 Schur 分解。 两者的同时分解可以以 A = Q●S●Z T 和 B = Q●T●Z 的形式实现 T 用于相同的正交矩阵Q和Z。 矩阵S和T是块上三角的，块的大小分别为 1×1 和 2×2。借助这些块，我们可以找到A和B的广义特征值。
模块scipy.linalg 中有四个功能，可为我们提供用于计算以下任何分解的工具：
建设者
|
描述
|
| --- | --- |
| scipy.linalg.schur(A) | 矩阵的舒尔分解 |
| scipy.linalg.rsf2csf(T, Z) | 从真实的 Schur 形式转换为复杂的 Schur 形式 |
| scipy.linalg.qz(A, B) | 两个矩阵的广义 Schur 分解 |
| scipy.linalg.hessenberg(A) | 矩阵的 Hessenberg 形式 |
函数hessenberg 为我们提供了任何 Schur 分解的第一步。 这是形式为 A = Q●U●Q H 的任何方阵A的因式分解，其中Q是一元的，U是上 Hessenberg 矩阵（所有条目在对角线以下均为零）。 该算法基于LAPACK例程GEHRD和GEBAL（用于计算U）和BLAS例程GER和GEMM（用于计算Q）的组合。
函数schur和qz是LAPACK例程GEES和GGES的包装器，用于分别计算平方矩阵的标准和广义 Schur 分解。 我们根据可选参数输出（我们将其设置为'real'或'complex'）选择报告复杂分解还是实分解。 我们还可以对矩阵表示中的特征值进行排序。 我们使用可选参数sort进行操作，并具有以下可能性：
None：如果我们不需要任何排序。 这是默认值。'lhp'：在左侧平面中。'rhp'：在右侧平面'iuc'：单位圆内'ouc'：单位圆以外func：任何称为func的可调用函数均可用于为用户提供自己的排序
摘要
在本章中，我们探讨了数值线性代数的基本原理，这是科学计算中所有过程的核心。 首先将重点放在矩阵和线性运算符的存储以及基本操作上。 我们详细探讨了所有不同的因式分解，重点放在它们的用法上，以找到矩阵方程或特征值问题的解决方案。 在本章中，我们着重将模块scipy.linalg和scipy.sparse的功能链接到库BLAS，LAPACK，ARPACK和SuperLU中的相应例程。 对于我们的实验，我们从现实问题中选择了有趣的矩阵，这些问题是从佛罗里达大学主办的广泛的稀疏矩阵集合中收集的。
在下一章中，我们将解决插值和最小二乘近似的问题。
二、插值和近似
近似理论指出了如何从某个预定类中的另一个函数中找到给定函数的最佳近似，以及这种近似的效果如何。 在本章中，我们将通过两个设置来探索该字段：插值和最小二乘近似。
动机
考虑一个气象实验，该实验测量海上矩形矩形网格上一组浮标的温度。 我们可以通过指示 16×16 位置的网格上的浮标的经度和纬度，以及它们上的随机温度在 36ºF 至 46ºF 之间来模拟这样的实验：
In [1]: import numpy as np, matplotlib.pyplot as plt, \
...: matplotlib.cm as cm; \
...: from mpl_toolkits.basemap import Basemap
In [2]: map1 = Basemap(projection='ortho', lat_0=20, lon_0=-60, \
...: resolution='l', area_thresh=1000.0); \
...: map2 = Basemap(projection='merc', lat_0=20, lon_0=-60, \
...: resolution='l', area_thresh=1000.0, \
...: llcrnrlat=0, urcrnrlat=45, \
...: llcrnrlon=-75, urcrnrlon=-15)
In [3]: longitudes = np.linspace(-60, -30, 16); \
...: latitudes = np.linspace(15, 30, 16); \
...: lons, lats = np.meshgrid(longitudes, latitudes); \
...: temperatures = 10\. * np.random.randn(16, 16) + 36.
In [4]: x1, y1 = map1(lons, lats); \
...: x2, y2 = map2(lons, lats)
In [5]: plt.rc('text', usetex=True); \
...: plt.figure()
In [6]: plt.subplot(121, aspect='equal'); \
...: map1.drawmeridians(np.arange(0, 360, 30)); \
...: map1.drawparallels(np.arange(-90, 90, 15)); \
...: map1.drawcoastlines(); \
...: map1.fillcontinents(color='coral'); \
...: map1.scatter(x1, y1, 15, temperatures, cmap=cm.gray)
In [7]: plt.subplot(122); \
...: map2.drawmeridians(np.arange(0, 360, 30)); \
...: map2.drawparallels(np.arange(-90, 90, 15)); \
...: map2.drawcoastlines(); \
...: map2.fillcontinents(color='coral'); \
...: C = map2.scatter(x2, y2, 15, temperatures, cmap=cm.gray); \
...: Cb = map2.colorbar(C, "bottom", size="5%", pad="2%"); \
...: Cb.set_label(r'$\mbox{}^{\circ} F$'); \
...: plt.show()
我们获得下图：
可以猜测这些浮标之间的温度（不完全准确，但至少在一定程度上），因为温度是地球表面的光滑函数。 让我们假设，我们需要借助三阶分段 2D 多项式进行近似，并且在块彼此相交处具有最大的平滑度。 当然，一个明显的挑战是浮标不是位于平面上，而是在一个很大的球体的表面上。 对于 SciPy 而言，这不是问题。
In [8]: from scipy.interpolate import RectSphereBivariateSpline \
...: as RSBS
In [9]: soln = RSBS(np.radians(latitudes), \
...: np.pi + np.radians(longitudes), \
...: temperatures)
In [10]: long_t = np.linspace(-60, -30, 180); \
....: lat_t = np.linspace(15, 30, 180); \
....: temperatures = soln(np.radians(lat_t), \
....: np.pi + np.radians(long_t))
In [11]: long_t, lat_t = np.meshgrid(long_t, lat_t); \
....: lo1, la1 = map1(long_t, lat_t); \
....: lo2, la2 = map2(long_t, lat_t)
In [12]: plt.figure()
Out[12]: <matplotlib.figure.Figure at 0x10ec28250>
In [13]: plt.subplot(121, aspect='equal'); \
....: map1.drawmeridians(np.arange(0, 360, 30)); \
....: map1.drawparallels(np.arange(-90, 90, 15)); \
....: map1.drawcoastlines(); \
....: map1.fillcontinents(color='coral'); \
....: map1.contourf(lo1, la1, temperatures, cmap=cm.gray)
Out[13]: <matplotlib.contour.QuadContourSet instance at 0x10f63d7e8>
In [14]: plt.subplot(122); \
....: map2.drawmeridians(np.arange(0, 360, 30)); \
....: map2.drawparallels(np.arange(-90, 90, 15)); \
....: map2.drawcoastlines(); \
....: map2.fillcontinents(color='coral'); \
....: C = map2.contourf(lo2, la2, temperatures, cmap=cm.gray); \
....: Cb = map2.colorbar(C, "bottom", size="5%", pad="2%"); \
....: Cb.set_label(r'$\mbox{}^{\circ} F$'); \
....: plt.show()
我们已经通过简单地通过将温度函数表示为与浮标位置上的温度值相符的分段多项式曲面来解决了这个问题。 这在技术上称为球面上的矩形节点网格上具有双变量样条的插值。
在其他情况下，只要结果函数与实际温度更紧密相关，这些位置的精确值就不是很重要。 在这种情况下，我们不希望执行插值，而是要计算具有相同功能类的元素的近似值。
让我们精确地定义两个设置：
插值问题需要三个要素：
- 有限域上的目标函数
f(x)（为方便起见，我们用x表示）。 - 域中的一组有限点：插值的节点，我们用
xi表示。 我们还将需要在这些节点处评估目标函数（可能还有其某些派生函数）。 在本章中，我们用yi表示这些。 - 插值族：具有与目标函数相同的输入/输出结构的函数。
插值问题的目标是通过匹配节点上目标函数的值，使插值成员对目标函数进行近似。
我们在以下设置中探索插值：
- 最近邻插值（任何维度）
- 通过分段线性函数进行插值（任意维度）
- 通过多项式进行单变量插值（Lagrange 和 Hermite 插值）
- 分段多项式的单变量插值
- 通过样条线进行单变量和双变量插值
- 径向基多元插值
提示
我们假设熟悉样条线的理论和应用。 有很多很好的入门资料，但我们建议您使用一些更实用的口味：
- Carl de Boor，花键实用指南。 施普林格（1978）。
- Paul Dierckx，曲线和样条曲线的曲面拟合。 牛津大学出版社，1993 年。
所有插值均通过模块scipy.interpolate进行。 特别是，与样条相关的是 Paul Dierckx FITPACK库中某些例程的一组包装。
要定义近似问题，我们需要以下四个要素：
- 有限域
x上的目标函数f(x)，该函数以尺寸为n的列向量作为输入并输出尺寸为m的列向量 - 一系列近似值
{g[a](x)}：具有与f(x)相同的输入/输出结构，其功能取决于参数a，该参数编码为尺寸为r的列向量 norm是一种用于测量x, ||f(x) - g(x)||任意两个给定功能之间距离的功能
近似问题的目标是找到相对于参数a最小化表达式||f(x) - g[a](x)||的近似成员。 这等效于相对于a为error功能F(a) = ||f(x) - g[a](x)||找到a（局部或全局）最小值。
我们说，如果近似族是一个基本元素的线性组合并且参数a作为系数，则近似是线性； 否则，我们将近似值称为非线性。
在本章中，我们将介绍以下设置中的近似函数：
- 通用线性最小二乘近似（通过求解线性方程组）
- 单变量和二变量样条的最小二乘近似/平滑
- 球上矩形网格上的样条最小二乘近似/平滑
- 通用非线性最小二乘近似（使用 Levenberg-Marquardt 迭代算法）
函数上下文中的最小二乘近似通过几个模块执行：
- 对于一般的线性最小二乘近似，可以始终将问题简化为线性方程组的解。 在这种情况下，我们在上一章研究的
scipy.linalg和scipy.sparse.linalg模块拥有我们需要的所有算法。 如前所述，所需的函数是 Fortran 库BLAS和LAPACK和C库SuperLU中几个例程的包装。 - 对于通过样条线进行线性最小二乘近似的特殊情况，
scipy.interpolate模块执行许多功能（针对所有不同情况），这些功能又是 Paul Diercks 的 Fortran 库FITPACK中的例程包装。 - 对于非线性最小二乘近似，我们使用来自
scipy.optimize模块的函数。 这些函数是 Fortran 库MINPACK中LMDIF和LMDER例程的包装。
提示
有关这些 Fortran 库的更多信息，可以从 Netlib 存储库中的它们的页面中获得很好的参考：
- FITPACK： netlib.org/dierckx/
- MINPACK： netlib.org/minpack/
- FFTPACK： netlib.org/fftpack/
可以从的创建者处找到SuperLU的最佳参考文献，网址为 http://crd-legacy.lbl.gov/~xiaoye/SuperLU/ 。
插补
我们有三种不同的实现方法来处理插值问题：
- 程序模式，该模式计算代表实际解决方案的一组数据点（以
ndarray的形式，具有所需的维数）。 - 在某些特殊情况下，功能性模式为我们提供了代表解决方案的
numpy功能。 - 面向对象的模式，为插值问题创建类。 不同的类具有不同的方法，具体取决于特定种类的插值器所享受的操作。 这种模式的优势在于，通过这些方法，我们可以从解决方案中获取更多信息：不仅是评估或表示，还包括诸如搜索根，计算导数和反导数，错误检查以及计算系数和结之类的相关操作。
代表我们的插值词的方式的选择取决于我们，这主要取决于我们需要多少精度，以及之后需要的信息/操作。
实施细节
在过程模式下，没有太多要添加到实现细节中的内容了。 对于每个插值问题，我们选择一个例程，向该例程中馈送节点xi，这些节点上目标函数（及其可能的导数）的值yi以及要评估插值的域x 。 在某些情况下，如果插值器需要更多结构，我们将提供更多信息。
功能实现甚至更简单：可用时，它们仅需要节点xi的值和在这些节点上的评估yi。
有几种通用的面向对象的插值类。 我们很少篡改它们，而是使用例程在内部创建和操作更合适的子类。 让我们简要介绍一下这些对象：
- 对于通用单变量插值，我们具有
_Interpolator1D类。 可以使用节点集xi以及这些节点上的目标函数的值yi进行初始化。 如有必要，我们也可以使用._set_dtype类方法强制使用yi的数据类型。 如果需要处理插值器的导数，可以使用_Interpolator1DWithDerivatives子类，并使用额外的类方法.derivatives来计算微分的评估。 - 对于样条小于或等于 5 的单变量插值，我们具有
InterpolatedUnivariateSpline类，而该类又是UnivariateSpline类的子类。 它们都是非常丰富的类，具有许多方法不仅可以评估样条曲线或其任何派生类，而且可以计算其派生类和反派生类的样条表示形式。 我们也有方法计算两个点之间的定积分。 也有一些方法可以返回结的位置，样条系数，残差甚至根。 我们至少使用节点xi和适合这些节点yi的值来初始化UnivariateSpline类中的对象。 我们可以选择使用样条曲线的度数初始化对象。 - 对于具有非结构化节点（节点不一定在矩形网格上）的双变量插值，
interp2d类是一个选项，该类可使用 1、3 或 5 阶双变量样条在二维中实现插值。 及其评估。 - 对于在矩形网格上具有节点的双变量样条插值，我们具有
RectBivariateSpline类（当与s = 0参数一起使用时），它是BivariateSpline类的子类。 反过来，BivariateSpline是基类_BivariateSplineBase的子类。 作为其单变量对应物，这是一个非常丰富的类，提供了许多用于评估，提取节点和系数以及计算体积积分或残差的方法。 - 对于多元插值，存在
NDInterpolatorBase类，具有三个子类：NearestNDInterpolator（用于最近邻插值），LinearNDInterpolator（用于分段线性插值）和CloughTocher2DInterpolator（实现分段三次，C1 平滑） ，二维的最小化曲率插值）。 - 为了在球面上的矩形网格上的一组节点上进行插值，需要使用
SphereBivariateSpline类的RectSphereBivariateSpline子类。 我们使用表示球体上节点位置的角度（theta和phi）以及相应的评估来初始化它。 - 对于具有径向函数的多元插值，我们具有
Rbf类。 它相当干燥，因为它仅允许使用评估方法。 它使用节点和评估进行初始化。
单变量插值
下表总结了用 SciPy 编码的不同单变量插值模式，以及我们可能用来解决它们的过程：
|插补模式
|
面向对象的实现
|
程序执行
|
| --- | --- | --- |
| 最近邻居 | interp1d(,kind='nearest') | |
| 拉格朗日多项式。 | BarycentricInterpolator | barycentric_interpolate |
| 埃尔米特多项式。 | KroghInterpolator | krogh_interpolate |
| 分段多项式。 | PiecewisePolynomial | piecewise_polynomial_interpolate |
| 分段线性 | interp1d(,kind='linear') | |
| 通用样条插值 | InterpolatedUnivariateSpline | splrep |
| 零阶样条 | interp1d(,kind='zero') | |
| 线性样条 | interp1d(,kind='slinear') | |
| 二次样条 | interp1d(,kind='quadratic') | |
| 三次样条 | interp1d(,kind='cubic') | |
| 邮编 | PchipInterpolator | pchip_interpolate |
最近邻插值
在一维函数的上下文中，最近邻插值提供了一种解决方案，该解决方案在由节点集的两个连续中点定义的每个子间隔上，围绕每个节点保持不变。 为了计算所需的插值值，我们使用kind='nearest'选项调用通用的scipy.interpolate.interp1d函数。 它仅使用可用的评估方法生成_Interpolator1D类的实例。
以下示例显示了其简单三角函数f(x) = sin(3*x)在 0 到 1 区间上的结果：
In [1]: import numpy as np, matplotlib.pyplot as plt; \
...: from scipy.interpolate import interp1d
In [2]: nodes = np.linspace(0, 1, 5); \
...: print nodes
[ 0\. 0.25 0.5 0.75 1\. ]
In [3]: def f(t): return np.sin(3 * t)
In [4]: x = np.linspace(0, 1, 100) # the domain
In [5]: interpolant = interp1d(nodes, f(nodes), kind='nearest')
In [6]: plt.rc('text', usetex=True)
...: plt.figure(); \
...: plt.axes().set_aspect('equal'); \
...: plt.plot(nodes, f(nodes), 'ro', label='Nodes'); \
...: plt.plot(x, f(x), 'r-', label=r'f(x)=\sin(3x)'); \
...: plt.plot(x, interpolant(x), 'b--', label='Interpolation'); \
...: plt.title("Nearest-neighbor approximation"); \
...: plt.ylim(-0.05, 1.05); \
...: plt.xlim(-0.5, 1.05); \
...: plt.show()
这将产生以下图形：
拉格朗日插值
在拉格朗日插值中，我们寻求与节点集处的目标函数相符的多项式。 在scipy.interpolate模块中，我们有三种方法可以解决此问题：
_Interpolator1D的BarycentricInterpolator子类基于有理函数近似实现了一种非常稳定的算法。 此类有一种评估方法，以及两种添加/更新运行中节点的方法：.add_xi和.set_yi。- 一种程序方案
barycentric_interpolate是上一类的语法糖，评估方法适用于规定的领域。 - 数值不稳定的函数方案
lagrange计算插值多项式的numpy.poly1d实例。 如果节点很少且明智地选择，则此方法可以使我们可靠地处理与目标函数相关的派生，积分和根求解问题。
让我们在臭名昭著的Runge示例中尝试这种插值模式：为函数 f（x）= 1 /（1 + x 2 [ *）*在从-5 到 5 的间隔中，具有两组相等分布的节点：
In [7]: from scipy.interpolate import BarycentricInterpolator, \
...: barycentric_interpolate, lagrange
In [8]: nodes = np.linspace(-5, 5, 11); \
...: x = np.linspace(-5,5,1000); \
...: print nodes
[-5\. -4\. -3\. -2\. -1\. 0\. 1\. 2\. 3\. 4\. 5.]
In [9]: def f(t): return 1\. / (1\. + t**2)
In [10]: interpolant = BarycentricInterpolator(nodes, f(nodes))
In [11]: plt.figure(); \
....: plt.subplot(121, aspect='auto'); \
....: plt.plot(x, interpolant(x), 'b--', \
....: label="Lagrange Interpolation"); \
....: plt.plot(nodes, f(nodes), 'ro', label='nodes'); \
....: plt.plot(x, f(x), 'r-', label="original"); \
....: plt.legend(loc=9); \
....: plt.title("11 equally distributed nodes")
Out[11]: <matplotlib.text.Text at 0x10a5fbe50>
BarycentricInterpolator类允许以最佳方式添加额外的节点并更新插值，而无需从头开始进行重新计算：
In [12]: newnodes = np.linspace(-4.5, 4.5, 10); \
....: print newnodes
[-4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5]
In [13]: interpolant.add_xi(newnodes, f(newnodes))
In [14]: plt.subplot(122, aspect='auto'); \
....: plt.plot(x, interpolant(x), 'b--', \
....: label="Lagrange Interpolation"); \
....: plt.plot(nodes, f(nodes), 'ro', label='nodes'); \
....: plt.plot(x, f(x), 'r-', label="original"); \
....: plt.legend(loc=8); \
....: plt.title("21 equally spaced nodes"); \
....: plt.show()
我们获得以下结果：
Runge示例显示了非常简单的插值的缺点之一。 尽管插值器可以在区间内部精确地近似函数，但它在端点处显示出非常大的偏差。
可以调用初始化插值的相同方法来请求有关它们的信息。 以下简短会议说明了这一点：
In [15]: print interpolant.xi
[-5\. -4\. -3\. -2\. -1\. 0\. 1\. 2\. 3\. 4\. 5\. -4.5 -3.5
–2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5]
In [16]: print interpolant.yi.squeeze()
[ 0.04 0.06 0.1 0.2 0.5 1\. 0.5 0.2 0.1 0.06 0.04
0.05 0.08 0.14 0.31 0.8 0.8 0.31 0.14 0.08 0.05]
程序方案具有更简单的语法，但是缺乏动态更新节点的灵活性：
In [17]: y = barycentric_interpolate(nodes, f(nodes), x)
该功能方案还具有简单的语法：
In [18]: g = lagrange(nodes, f(nodes)); \
....: print g
10 9 8 7 6
3.858e-05 x + 6.268e-19 x - 0.002149 x + 3.207e-17 x + 0.04109 x
5 4 3 2
+ 5.117e-17 x - 0.3302 x - 2.88e-16 x + 1.291 x - 1.804e-16 x
埃尔米特插值
Hermite插值的目的是计算与目标函数及其在一组有限节点中的某些导数一致的多项式。 我们通过两种方案以数字方式完成此任务：
_Interpolator1DWithDerivatives的子类KroghInterpolator，具有.derivative方法来计算插值的任何导数的表示形式，并具有.derivatives方法来对其进行评估。krogh_interpolate函数，是上一类的语法糖，并且在规定的域上应用了评估方法。
让我们以伯恩斯坦的示例展示这些例程：在从-1到1的区间中，使用十个平均分布的节点，将 Hermite 插值计算为绝对值函数，并在每个节点上提供一个导数。
提示
节点需要以递增顺序进行馈送。 对于每个我们要派生导数的节点，我们都会根据需要重复该节点多次。 对于xi中节点的每次出现，我们将yi上的函数及其派生子的评估放在相同的入口级别。
In [19]: from scipy.interpolate import KroghInterpolator
In [20]: nodes = np.linspace(-1, 1, 10); \
....: x = np.linspace(-1, 1, 1000)
In [21]: np.set_printoptions(precision=3, suppress=True)
In [22]: xi = np.repeat(nodes, 2); \
....: print xi; \
....: yi = np.ravel(np.dstack((np.abs(nodes), np.sign(nodes)))); \
....: print yi
[-1\. -1\. -0.778 -0.778 -0.556 -0.556 -0.333 -0.333 -0.111
-0.111 0.111 0.111 0.333 0.333 0.556 0.556 0.778 0.778
1\. 1\. ]
[ 1\. -1\. 0.778 -1\. 0.556 -1\. 0.333 -1\. 0.111
-1\. 0.111 1\. 0.333 1\. 0.556 1\. 0.778 1\.
1\. 1\. ]
In [23]: interpolant = KroghInterpolator(xi, yi)
In [24]: plt.figure(); \
....: plt.axes().set_aspect('equal'); \
....: plt.plot(x, interpolant(x), 'b--', \
....: label='Hermite Interpolation'); \
....: plt.plot(nodes, np.abs(nodes), 'ro'); \
....: plt.plot(x, np.abs(x), 'r-', label='original'); \
....: plt.legend(loc=9); \
....: plt.title('Bernstein example'); \
....: plt.show()
这给出了下图：
分段多项式插值
通过规定几个多项式的阶数和节点的有限集，我们可以构造一个插值器，该插值器在两个连续节点之间的每个子间隔上具有所需顺序的多项式弧。 我们可以通过以下步骤构造具有此特征的插值值：
PiecewisePolynomial的子类PiecewisePolynomial，具有评估插值值及其派生值或附加新节点的方法- 对于分段线性插值的特殊情况，
interp1d实用程序仅使用评估方法创建_Interpolator1D类的实例 piecewise_polynomial_interpolate函数，它是PiecewisePolynomial类的语法糖，其评估方法适用于规定的域
让我们回顾一下本节中的第一个示例。 首先，我们尝试使用interp1d进行分段线性插值。 其次，我们使用PiecewisePolynomial在每个节点上应用具有正确导数的分段二次插值（所有分段的阶数均为 2）。
In [25]: from scipy.interpolate import PiecewisePolynomial
In [26]: nodes = np.linspace(0, 1, 5); \
....: x = np.linspace(0, 1, 100)
In [27]: def f(t): return np.sin(3 * t)
In [28]: interpolant = interp1d(nodes, f(nodes), kind='linear')
In [29]: plt.figure(); \
....: plt.subplot(121, aspect='equal'); \
....: plt.plot(x, interpolant(x), 'b--', label="interpolation"); \
....: plt.plot(nodes, f(nodes), 'ro'); \
....: plt.plot(x, f(x), 'r-', label="original"); \
....: plt.legend(loc=8); \
....: plt.title("Piecewise Linear Interpolation")
Out[29]: <matplotlib.text.Text at 0x107be0390>
In [30]: yi = np.zeros((len(nodes), 2)); \
....: yi[:,0] = f(nodes); \
....: yi[:,1] = 3 * np.cos(3 * nodes); \
....: print yi
[[ 0\. 3\. ]
[ 0.682 2.195]
[ 0.997 0.212]
[ 0.778 -1.885]
[ 0.141 -2.97 ]]
In [31]: interpolant = PiecewisePolynomial(nodes, yi, orders=2)
In [32]: plt.subplot(122, aspect='equal'); \
....: plt.plot(x, interpolant(x), 'b--', label="interpolation"); \
....: plt.plot(nodes, f(nodes), 'ro'); \
....: plt.plot(x, f(x), 'r-', label="original"); \
....: plt.legend(loc=8); \
....: plt.title("Piecewise Quadratic interpolation"); \
....: plt.show()
这给出了下图：
在此图像中，分段二次插值和原始函数实际上是无法区分的。 我们需要去计算差值的绝对值（函数以及它的一阶和二阶导数）才能真正实现计算误差。 以下是近似这些误差的粗略计算，并说明了.derivatives方法的使用：
In [33]: np.abs(f(x) - interpolant(x)).max()
Out[33]: 0.0093371930045896279
In [34]: f1prime = lambda t: 3 * np.cos(3 * t); \
....: np.abs(f1prime(x) - interpolant.derivatives(x)).max()
Out[34]: 10.589218385920123
In [35]: f2prime = lambda t: -9 * np.sin(3 * x); \
....: np.abs(f2prime(x) - interpolant.derivatives(x,der=2)).max()
Out[35]: 9.9980773091170505
分段多项式近似的一个很大优点是在不同子区间上使用不同次数的多项式的灵活性。 例如，对于相同的节点集，我们可以在第一个和最后一个子间隔和三次中使用其他行：
In [36]: interpolant = PiecewisePolynomial(nodes, yi, \
....: orders=[1,3,3,1])
实现这种插值方案的另一个巨大优势是，我们可以轻松添加新节点，而无需从头开始进行重新计算。 例如，要在最后一个节点之后添加新节点，我们发出：
In [37]: interpolant.append(1.25, np.array([f(1.25)]))
样条插值
单变量样条曲线是分段多项式的一种特殊情况。 它们在多项式连接的地方具有高度的平滑度。 可以将这些函数编写为基本样条的线性组合，并且对于给定的度数，平滑度和节点集，只需很少的支持即可。
通过interp1d功能和适当的kind选项，可以使用最多 5 个等级的splines使用单变量样条插值。 此函数使用相应的类方法创建_Interpolator1DWithDerivatives的实例。 通过对 Fortran 库FITPACK中的例程的调用来执行计算。 以下示例显示了不同的可能性：
In [38]: splines = ['zero', 'slinear', 'quadratic', 'cubic', 4, 5]; \
....: g = KroghInterpolator([0,0,0,1,1,1,2,2,2,3,3,3], \
....: [10,0,0,1,0,0,0.25,0,0,0.125,0,0]); \
....: f = lambda t: np.log1p(g(t)); \
....: x = np.linspace(0,3,100); \
....: nodes = np.linspace(0,3,11)
In [39]: plt.figure()
In [40]: for k in xrange(6):
....: interpolant = interp1d(nodes, f(nodes), \
....: kind = splines[k])
....: plt.subplot(2,3,k+1, aspect='equal')
....: plt.plot(nodes, f(nodes), 'ro')
....: plt.plot(x, f(x), 'r-', label='original')
....: plt.plot(x, interpolant(x), 'b--', \
....: label='interpolation')
....: plt.title('{0} spline'.format(splines[k]))
....: plt.legend()
In [41]: plt.show()
这给出了下图：
注意
零样条非常类似于最近邻近似，尽管在这种情况下，插值在两个连续节点的每个选择之间都是恒定的。 线性样条曲线与分段线性插值完全相同。 但是，通过样条曲线执行插值的算法较慢。
对于任何给定的问题设置，都有许多不同的样条插值，它们具有相同的度数，节点和评估。 例如，输出还取决于结的位置和数量。 不幸的是，interp1d功能仅允许控制节点和值。 该算法在结计算方面使用了最简单的设置。
例如，请注意，上一个示例中的三次样条插值不能保留目标函数的单调性。 在这种情况下，可以通过谨慎地限制导数或结的位置来强制插值的单调性。 我们有一个特殊功能，可以通过使用 Fritsch-Carlson 算法实现的分段单调三次 Hermite 插值（ PCHIP ）来完成此任务。 这个简单的算法可以通过_Interpolator1DWithDerivatives的PchipInterpolator子类来实现，也可以通过其等效的过程函数pchip_interpolate来实现。
In [42]: from scipy.interpolate import PchipInterpolator
In [43]: interpolant = PchipInterpolator(nodes, f(nodes))
In [44]: plt.figure(); \
....: plt.axes().set_aspect('equal'); \
....: plt.plot(nodes, f(nodes), 'ro'); \
....: plt.plot(x, f(x), 'r-', label='original'); \
....: plt.plot(x, interpolant(x), 'b--', label='interpolation'); \
....: plt.title('PCHIP interpolation'); \
....: plt.legend(); \
....: plt.show()
这给出了下图：
通用样条插值由InterpolatedUnivariateSpline类处理，在这里我们可以实际控制所有影响样条质量的参数。 在这种情况下，所有计算都由 Fortran 库FITPACK中的例程包装程序执行。 通过scipy.interpolate模块中的一组功能，可以以程序方式访问这些包装器。 下表显示了类方法，相应的过程函数以及它们调用的FITPACK例程之间的匹配：
操作方式
|
面向对象的实现
|
程序
|
FITPACK
|
| --- | --- | --- | --- |
| 插值的实例化 | InterpolatedUnivariateSpline | splrep | CURFIT |
| 报告花键的结 | object.get_knots() | splrep | |
| 报告样条系数 | object.get_coeffs() | splrep | CURFIT |
| 花键的评估 | object() | splev | SPLEV |
| 衍生物 | object.derivative() | splder | |
| 衍生品评估 | object.derivatives() | splev, spalde | SPLDER, SPALDE |
| 反导数 | object.antiderivative() | splantider | |
| 定积分 | object.integral() | splint | SPLINT |
| 根（用于三次样条） | object.roots() | sproot | SPROOT |
提示
所获得的值。 get_coeffs方法是作为 B 样条的线性组合的样条系数。
让我们展示如何通过计算 5 阶相应插值样条曲线的积分来近似估算目标函数图下的面积。
In [45]: from scipy.interpolate import InterpolatedUnivariateSpline \
....: as IUS
In [46]: interpolant = IUS(nodes, f(nodes), k=5)
In [47]: area = interpolant.integral(0,3); \
....: print area
2.14931665485
多元插值
样条曲线的双变量插值可以通过scipy.interpolate模块中的interp2d执行。 这是一个非常简单的类，仅允许使用求值方法，并具有三种基本的样条插值模式编码：线性，三次和五次。 它无法控制节或重。 为了创建双变量样条的表示，interp2d函数从库FITPACK中调用 Fortran 例程SURFIT（可悲的是，实际上并不意味着执行插值！）。 为了用数字方式评估样条，模块调用例程BISPEV。
让我们通过示例显示interp2d的用法。 我们首先构造一个有趣的双变量函数，以在其域上插值 100 个节点的随机选择，并提供可视化效果：
In [1]: import numpy as np, matplotlib.pyplot as plt; \
...: from mpl_toolkits.mplot3d.axes3d import Axes3D
In [2]: def f(x, y): return np.sin(x) + np.sin(y)
In [3]: t = np.linspace(-3, 3, 100); \
...: domain = np.meshgrid(t, t); \
...: X, Y = domain; \
...: Z = f(*domain)
In [4]: fig = plt.figure(); \
...: ax1 = plt.subplot2grid((2,2), (0,0), aspect='equal'); \
...: p = ax1.pcolor(X, Y, Z); \
...: fig.colorbar(p); \
...: CP = ax1.contour(X, Y, Z, colors='k'); \
...: ax1.clabel(CP); \
...: ax1.set_title('Contour plot')
In [5]: nodes = 6 * np.random.rand(100, 2) - 3; \
...: xi = nodes[:, 0]; \
...: yi = nodes[:, 1]; \
...: zi = f(xi, yi)
In [6]: ax2 = plt.subplot2grid((2,2), (0,1), aspect='equal'); \
...: p2 = ax2.pcolor(X, Y, Z); \
...: ax2.scatter(xi, yi, 25, zi) ; \
...: ax2.set_xlim(-3, 3); \
...: ax2.set_ylim(-3, 3); \
...: ax2.set_title('Node selection')
In [7]: ax3 = plt.subplot2grid((2,2), (1,0), projection='3d', \
...: colspan=2, rowspan=2); \
...: ax3.plot_surface(X, Y, Z, alpha=0.25); \
...: ax3.scatter(xi, yi, zi, s=25); \
...: cset = ax3.contour(X, Y, Z, zdir='z', offset=-4); \
...: cset = ax3.contour(X, Y, Z, zdir='x', offset=-5); \
...: ax3.set_xlim3d(-5, 3); \
...: ax3.set_ylim3d(-3, 5); \
...: ax3.set_zlim3d(-4, 2); \
...: ax3.set_title('Surface plot')
In [8]: fig.tight_layout(); \
...: plt.show()
我们获得下图：
然后可以使用以下节点执行分段线性插值：
In [9]: from scipy.interpolate import interp2d
In [10]: interpolant = interp2d(xi, yi, zi, kind='linear')
In [11]: plt.figure(); \
....: plt.axes().set_aspect('equal'); \
....: plt.pcolor(X, Y, interpolant(t, t)); \
....: plt.scatter(xi, yi, 25, zi); \
....: CP = plt.contour(X, Y, interpolant(t, t), colors='k'); \
....: plt.clabel(CP); \
....: plt.xlim(-3, 3); \
....: plt.ylim(-3, 3); \
....: plt.title('Piecewise linear interpolation'); \
....: plt.show()
尽管其名称为interp2d，但此过程不是实际的插值，而是试图拟合数据的粗略近似。 通常，每次运行此代码，都将获得不同质量的结果。 幸运的是，其余的插补例程并非如此！
注意
如果节点的位置不是最优的，我们很可能会收到警告：
Warning: No more knots can be added because the number of B-spline coefficients
already exceeds the number of data points m. Probably causes: either
s or m too small. (fp>s)
kx,ky=1,1 nx,ny=11,14 m=100 fp=0.002836 s=0.000000
提示
注意，在前面的示例中，对插值值的评估是通过调用两个一维数组来执行的。 通常，要评估在矩形网格上用interp2d计算的插值g，该插值可以实现为两个一维数组（大小m的tx和大小n的ty ]），我们发布g(tx, ty)； 这给出了一个二维数组，大小为 m x n 。
当然，结果的质量与节点的密度和结构密切相关。 增加其数量或将其位置置于矩形网格上会改善问题。 在节点形成矩形网格的情况下，借助于RectBivariateSpline类可以使用更快，更准确的方法进行实际插值。 此函数是FITPACK库中 Fortran 例程REGRID的包装。
现在让我们在矩形网格上选择 100 个节点并重新计算，如下所示：
In [12]: ti = np.linspace(-3, 3, 10); \
....: xi, yi = np.meshgrid(ti, ti); \
....: zi = f(xi, yi)
In [13]: from scipy.interpolate import RectBivariateSpline
In [14]: interpolant = RectBivariateSpline(ti, ti, zi, kx=1, ky=1)
In [15]: plt.figure(); \
....: plt.axes().set_aspect('equal'); \
....: plt.pcolor(X, Y, interpolant(t, t)); \
....: CP = plt.contour(X, Y, interpolant(t, t), colors='k'); \
....: plt.clabel(CP); \
....: plt.scatter(xi, yi, 25, zi); \
....: plt.xlim(-3, 3); \
....: plt.ylim(-3, 3); \
....: plt.title('Piecewise linear interpolation, \
....: rectangular grid'); \
....: plt.show()
提示
与interp2d的情况一样，要评估在矩形网格上用RectBivariateSpline计算的插值g，该插值可以实现为两个一维数组（m大小为m和[ n大小的ty），我们发布g(tx, ty)； 这给出了一个二维数组，大小为 m x n 。
现在给出一个实际的插值：
图下的体积积分非常准确（给定域中目标函数的实际积分为零）：
In [16]: interpolant.integral(-3, 3, -3, 3)
Out[16]: 2.636779683484747e-16
在这种情况下，让我们检查一下从此类中获得的一些不同的信息：
-
插值度：
In [17]: interpolant.degrees Out[17]: (1, 1) -
返回的样条曲线近似值的平方残差之和：
In [18]: interpolant.fp Out[18]: 0.0 In [19]: interpolant.get_residual() Out[19]: 0.0 -
插值系数：
In [20]: np.set_printoptions(precision=5, suppress=True) In [21]: print interpolant.get_coeffs() [-0.28224 -0.86421 -1.13653 -0.98259 -0.46831 0.18607 0.70035 0.85429 0.58197 0\. -0.86421 -1.44617 -1.71849 -1.56456 -1.05028 -0.39589 0.11839 0.27232 0\. -0.58197 -1.13653 -1.71849 -1.99082 -1.83688 -1.3226 -0.66821 -0.15394 0\. -0.27232 -0.85429 -0.98259 -1.56456 -1.83688 -1.68294 -1.16867 -0.51428 0\. 0.15394 -0.11839 -0.70035 -0.46831 -1.05028 -1.3226 -1.16867 -0.65439 -0\. 0.51428 0.66821 0.39589 -0.18607 0.18607 -0.39589 -0.66821 -0.51428 -0\. 0.65439 1.16867 1.3226 1.05028 0.46831 0.70035 0.11839 -0.15394 0\. 0.51428 1.16867 1.68294 1.83688 1.56456 0.98259 0.85429 0.27232 0\. 0.15394 0.66821 1.3226 1.83688 1.99082 1.71849 1.13653 0.58197 0\. -0.27232 -0.11839 0.39589 1.05028 1.56456 1.71849 1.44617 0.86421 0\. -0.58197 -0.85429 -0.70035 -0.18607 0.46831 0.98259 1.13653 0.86421 0.28224] -
结的位置：
In [22]: interpolant.get_knots() (array([-3\. , -3\. , -2.33333, -1.66667, -1\. , -0.33333, 0.33333, 1\. , 1.66667, 2.33333, 3\. , 3\. ]), array([-3\. , -3\. , -2.33333, -1.66667, -1\. , -0.33333, 0.33333, 1\. , 1.66667, 2.33333, 3\. , 3\. ]))
使用分段三次样条曲线可以获得更平滑的结果。 在前面的示例中，我们可以通过设置kx = 3和ky = 3来完成此任务：
In [23]: interpolant = RectBivariateSpline(ti, ti, zi, kx=3, ky=3)
In [24]: fig = plt.figure(); \
....: ax = fig.add_subplot(121, projection='3d',aspect='equal'); \
....: ax.plot_surface(X, Y, interpolant(t, t), alpha=0.25, \
....: rstride=5, cstride=5); \
....: ax.scatter(xi, yi, zi, s=25); \
....: C = ax.contour(X, Y, interpolant(t, t), zdir='z', \
....: offset=-4); \
....: C = ax.contour(X, Y, interpolant(t, t), zdir='x',\
....: offset=-5); \
....: ax.set_xlim3d(-5, 3); \
....: ax.set_ylim3d(-3, 5); \
....: ax.set_zlim3d(-4, 2); \
....: ax.set_title('Cubic interpolation, RectBivariateSpline')
在所有可能的三次样条插补中，有一种特殊情况可以最大程度地减小曲率。 通过收敛到解决方案的巧妙的迭代算法，我们可以实现这种特殊情况。 它依赖于以下三个关键概念：
- 使用节点作为顶点的域的 Delaunay 三角剖分
- 使用 Cough-Tocher 方案在每个三角形上支持的 Bezier 三次多项式
- 估计和施加梯度以最小化曲率
可以通过scipy.interpolate模块中的CloughTocher2dInterpolator功能或通过带有method='cubic'选项的同一模块中的黑匣子功能griddata来实现。 让我们比较输出：
In [25]: from scipy.interpolate import CloughTocher2DInterpolator
In [26]: nodes = np.dstack((np.ravel(xi), np.ravel(yi))).squeeze(); \
....: zi = f(nodes[:, 0], nodes[:, 1])
In [27]: interpolant = CloughTocher2DInterpolator(nodes, zi)
In [28]: ax = fig.add_subplot(122, projection='3d', aspect='equal'); \
....: ax.plot_surface(X, Y, interpolant(X, Y), alpha=0.25, \
....: rstride=5, cstride=5); \
....: ax.scatter(xi, yi, zi, s=25); \
....: C = ax.contour(X, Y, interpolant(X, Y), zdir='z', \
....: offset=-4); \
....: C = ax.contour(X, Y, interpolant(X, Y), zdir='x', \
....: offset=-5); \
....: ax.set_xlim3d(-5, 3); \
....: ax.set_ylim3d(-3, 5); \
....: ax.set_zlim3d(-4, 2); \
....: ax.set_title('Cubic interpolation, \
....: CloughTocher2DInterpolator'); \
....: plt.show()
提示
与interp2d和RectBivariateSpline的情况不同，要评估用CloughTocher2DInterpolator在矩形网格X, Y = domain上计算的插值g，我们发布g(X, Y)或g(*domain)。
这给出了下图：
黑盒程序函数griddata还允许我们访问多维的分段线性插值，以及多维最近邻插值。
In [29]: from scipy.interpolate import griddata
In [30]: Z = griddata(nodes, zi, (X, Y), method='nearest')
In [31]: plt.figure(); \
....: plt.axes().set_aspect('equal'); \
....: plt.pcolor(X, Y, Z); \
....: plt.colorbar(); \
....: plt.title('Nearest-neighbors'); \
....: plt.show()
这为我们提供了以下不太令人印象深刻的图表：
还需要考虑另一种插值模式：径向基函数插值。 此处的目的是针对同一函数g插入以(xk, yk)为中心的形式fk(x,y) = g(sqrt((x-xk)**2 + (y-yk)**2))的径向函数的线性组合。 我们可以在七个标准功能g (listed below)中进行选择，甚至可以选择我们自己的：
'multiquadric': g(r) = sqrt((r/self.epsilon)**2 + 1)'inverse': g(r) = 1.0/sqrt((r/self.epsilon)**2 + 1)'gaussian': g(r) = exp(-(r/self.epsilon)**2)'linear': g(r) = r'cubic': g(r) = r**3'quintic': g(r) = r**5'thin_plate': g(r) = r**2 * log(r)
通过Rbf类执行该实现。 可以照常使用节点及其评估值对其进行初始化。 我们还需要包括对径向函数的选择，并且如果需要，还可以包括影响bumps大小的epsilon参数的值。
让我们运行几个插值：首先，通过具有标准偏差epsilon = 2.0的径向高斯曲线，然后使用基于sinc的径向函数。 让我们还回到随机节点：
In [32]: from scipy.interpolate import Rbf
In [33]: nodes = 6 * np.random.rand(100, 2) - 3; \
....: xi = nodes[:, 0]; \
....: yi = nodes[:, 1]; \
....: zi = f(xi, yi)
In [34]: interpolant = Rbf(xi, yi, zi, function='gaussian', \
....: epsilon=2.0)
In [35]: plt.figure(); \
....: plt.subplot(121, aspect='equal'); \
....: plt.pcolor(X, Y, interpolant(X, Y)); \
....: plt.scatter(xi, yi, 25, zi); \
....: plt.xlim(-3, 3); \
....: plt.ylim(-3, 3)
Out[35]: (-3, 3)
In [36]: interpolant = Rbf(xi, yi, zi, function = np.sinc)
In [37]: plt.subplot(122, aspect='equal'); \
....: plt.pcolor(X, Y, interpolant(X, Y)); \
....: plt.scatter(xi, yi, 25, zi); \
....: plt.xlim(-3, 3); \
....: plt.ylim(-3, 3); \
....: plt.show()
尽管节点是非结构化的，但这仍然提供了两个非常准确的插值：
需要考虑的最后一种情况是球体上矩形网格上的二元样条插值。 我们使用RectSphereBivariateSpline函数获得此插值，该函数通过调用FITPACK例程SPGRID（用于计算样条的表示）和BISPEV（用于评估）来实例化SphereBivariateSpline的子类。
实施和评估让人联想到 Fortran 编码方法：
- 要计算样条曲线表示，我们发出
RectSphereBivariateSpline(u, v, data)命令，其中u和v都是严格增加的正值的一维数组，这些正值分别代表了纬度和经度的角度（以弧度为单位）。 节点的位置。 - 在评估时，如果需要在尺寸为 m x n 的细化网格上以二维形式表示插值，我们将发布
object(theta, phi)，其中theta和phi是一维且严格增加的，并且必须包含在以上u和v定义的域中。 输出（尽管您的文档说了什么）是 m x n 数组。
最小二乘近似
在数值上，陈述最小二乘范数的近似问题相对简单。 这是本节的主题。
线性最小二乘近似
在线性最小二乘近似的情况下，始终可以通过求解线性方程组来减少问题，如以下示例所示：
考虑从 0 到 1 的区间中的正弦函数 f（x）= sin（x）。我们选择二阶多项式作为近似值： {a 0 +一个 1 x +一个 2 x 2 } 。 为了计算值 [a 0 ， 1 一个使这个问题最小化的 2 ] ，我们首先形成一个 3×3 矩阵，其中包含成对的点积（两个乘积的积分 基本功能 {1，x，x 2 } 的基本功能）。 由于此问题的性质，我们获得了阶数为 3 的希尔伯特矩阵：
[ < 1, 1 > < 1, x > < 1, x^2 > ] [ 1 1/2 1/3 ]
[ < x, 1 > < x, x > < x, x^2 > ] = [ 1/2 1/3 1/4 ]
[ < x^2, 1 > < x^2, x > < x^2, x^2 > ] [ 1/3 1/4 1/5 ]
系统的右侧是给定间隔中正弦函数与每个基本函数的点积的列向量：
[ < sin(x), 1 > ] [ 1 - cos(1) ]
[ < sin(x), x > ] = [ sin(1) - cos(1) ]
[ < sin(x), x^2 > ] [ 2*sin(1) + cos(1) - 2 ]
我们按以下方式计算系数和相应的近似多项式：
In [1]: import numpy as np, scipy.linalg as spla, \
...: matplotlib.pyplot as plt
In [2]: A = spla.hilbert(3); \
...: b = np.array([1-np.cos(1), np.sin(1)-np.cos(1), \
...: 2*np.sin(1)+ np.cos(1)-2])
In [3]: spla.solve(A, b)
Out[3]: array([-0.00746493, 1.09129978, -0.2354618 ])
In [4]: poly1 = np.poly1d(spla.solve(A, b)[::-1]); \
...: print poly1
2
-0.2355 x + 1.091 x - 0.007465
通常，要解决基于r元素的线性最小二乘近似问题，我们需要解决具有r方程和r不定式的线性方程组的基本系统。 尽管其表面上很简单，但该方法远非完美。 以下是两个主要原因：
- 该系统可能状况不佳，就像前面的示例一样。
- 系数是非永久性的。 系数的值在很大程度上取决于
r。 增加问题的范围会导致产生一组新系数，与先前的系数不同。
注意
有多种方法可以修复系统的故障。 一种标准程序是使用 Gram-Schmidth 和改进的 Gram-Schmidt 正交化方法从原始对象构建正交基础。 这个主题超出了本专论的范围，但是可以在 Walter Gautschi，Birkhäuser，1997 年的数值分析一书的第 1 章中阅读这些方法的良好参考。
始终提供简单线性系统的基础是 B 样条曲线。 所有涉及的系统都是三对角线的，因此无需复杂的操作即可轻松解决。 scipy.interpolate模块中编码的面向对象系统使我们能够在内部执行所有这些计算。 这是所涉及的类和子类的简要枚举：
UnivariateSpline用于一维的样条曲线或任何维的曲线样条。 我们很少直接使用此类，而是使用LSQUnivariateSpline子类。BivariateSpline用于样条线的曲线，表示放置在矩形上的节点上的曲面。 作为其单变量对应项，不得直接使用此类。 相反，我们利用了LSQBivariateSpline子类。- 样条曲线的
SphereBivariateSpline表示放置在球体上的节点上的曲面。 计算必须改为通过LSQSphereBivariateSpline子类进行。
提示
在这三种情况下，基类及其方法在插值问题中都是相对应的。 有关更多信息，请参见插值部分。
让我们通过一些示例来说明这些面向对象的技术：
在最小二乘意义上，用三次样条（k = 3）在相同的域上近似具有相同的sine函数。 首先，请注意，我们必须提供一个边界框，该域上的一组结点以及最小二乘近似的权重w。 我们还可以提供可选的平滑度参数s。 如果s = 0，则获得插值，而对于s的值较大，则可以实现所得样条曲线的不同程度的平滑度。 为了获得可靠的（加权的）最小二乘近似，较好的选择是s = len(w)（默认情况下，例程会执行此操作）。 还请注意，计算出的误差有多小：
In [5]: f = np.sin; \
...: x = np.linspace(0,1,100); \
...: knots = np.linspace(0,1,7)[1:-1]; \
...: weights = np.ones_like(x)
In [6]: from scipy.interpolate import LSQUnivariateSpline
In [7]: approximant = LSQUnivariateSpline(x, f(x), knots, k=3, \
...: w = weights, bbox = [0, 1])
In [8]: spla.norm(f(x) - approximant(x))
Out[8]: 3.370175009262551e-06
提示
计算此近似误差的更方便方法是使用。 get_residual方法如下：
In [9]: approximant.get_residual()**(.5)
Out[9]: 3.37017500928446e-06
在[-3,3] x [-3,3]域上近似二维函数sin(x)+sin(y)。 我们首先选择域的表示形式，网格上的一组 100 个合适的节点以及权重。 由于LSQBivariateSpline函数的所有输入都必须是一维数组，因此我们在调用近似函数之前执行相应的转换：
In [10]: def f(x, y): return np.sin(x) + np.sin(y); \
....: t = np.linspace(-3, 3, 100); \
....: domain = np.meshgrid(t, t); \
....: X, Y = domain; \
....: Z = f(*domain)
In [11]: X = X.ravel(); \
....: Y = Y.ravel(); \
....: Z = Z.ravel()
In [12]: kx = np.linspace(-3,3,12)[1:-1]; \
....: ky = kx.copy(); \
....: weights = np.ones_like(Z);
In [13]: from scipy.interpolate import LSQBivariateSpline
In [14]: approximant = LSQBivariateSpline(X, Y, Z, kx, kx, \
....: w = weights)
In [15]: approximant.get_residual()
Out[15]: 0.0
提示
也可以使用RectBivariateSpline功能执行此计算。 为了实现最小二乘插值，我们提供了足够大的节点（而不是结，因为它们将自动计算），权重w和平滑度参数s。 s = len(w)是一个不错的选择。
非线性最小二乘近似
在非线性最小二乘近似的情况下，我们通常没有简单矩阵表示的奢侈。 取而代之的是，我们使用迭代过程的两个变体，即scipy.optimize模块中托管的 Levenberg-Marquardt 算法。 可以通过leastsq包装器调用这两个版本，它们对应于 Fortran 库MINPACK中的LMDER和LMDIF例程。
下表列出了此功能的所有选项：
|选项
|
描述
|
| --- | --- |
| func | 错误功能F(a) |
| x0 | 最小化的起始估算值，大小为r |
| args | 作为元组的func的额外参数 |
| Dfun | 表示func的雅可比矩阵的函数 |
| full_output | 布尔型 |
| col_deriv | 布尔型 |
| ftol | 期望的相对误差平方和 |
| xtol | 近似解中所需的相对误差 |
| gtol | func和Dfun的列之间需要正交 |
| maxfev | 最大通话次数。 如果为零，则通话次数为100*(r+1) |
| epsfcn | 如果Dfun=None，我们可以指定一个浮点值作为雅可比矩阵的前向差分近似中的步长 |
| factor | 浮点值介于 0.1 到 100 之间，指示初始步长边界 |
| diag | 每个变量的比例因子 |
当我们有一个可靠的误差函数雅可比矩阵时，将使用算法的第一个变体。 如果未提供，则使用该算法的第二种变体，该变体通过前向差近似雅可比行列式。 我们通过几个示例来说明这两种变体。
让我们开始使用该方法重新审视先前的示例，以了解用法和准确性的差异。 我们将把计算重点放在从 0 到 1 的间隔的分区上，该分区具有 100 个均匀间隔的点：
In [16]: from scipy.optimize import leastsq
In [17]: def error_function(a):
....: return a[0] + a[1] * x + a[2] * x**2 - np.sin(x)
In [18]: def jacobian(a): return np.array([np.ones(100), x, x**2])
In [19]: coeffs, success = leastsq(error_function, np.zeros((3,)))
In [20]: poly2 = np.poly1d(coeffs[::-1]); print poly2
2
-0.2354 x + 1.091 x - 0.007232
In [21]: coeffs, success = leastsq(error_function, np.zeros((3,)), \
....: Dfun = jacobian, col_deriv=True)
In [22]: poly3 = np.poly1d(coeffs[::-1]); \
...: print poly3
2
-0.2354 x + 1.091 x - 0.007232
In [23]: map(lambda f: spla.norm(np.sin(x) - f(x)), \
....: [poly1, poly2, poly3])
Out[23]:
[0.028112146265269783, 0.02808377541388009, 0.02808377541388009]
提示
scipy.optimize模块中还有另一个函数可以执行非线性最小二乘近似：curve_fit。 它使用相同的算法，但不是误差函数，而是将其泛型近似值g[a](x)以及自变量x的适当域以及目标函数f的输出传递给该域 。 我们也确实需要输入初始估计。 输出与所需系数一起是所述系数的协方差的估计。
In [23]: from scipy.optimize import curve_fit
In [24]: def approximant(t, a, b, c):
....: return a + b*t + c*t**2
In [25]: curve_fit(approximant, x, np.sin(x), \
....: np.ones((3,)))
(array([-0.007232 , 1.09078356, -0.23537796]),
array([[ 7.03274163e-07, -2.79884256e-06,
2.32064835e-06],
[ -2.79884256e-06, 1.50223585e-05,
-1.40659702e-05],
[ 2.32064835e-06, -1.40659702e-05,
1.40659703e-05]]))
在本节中，我们仅专注于leastsq功能。 这两个功能的目标和编码是相同的，但是leastsq可以按需提供更多信息，并可以更好地控制 Levenberg-Marquardt 算法的不同参数。
现在让我们尝试一些实际的非线性问题：
在第一个示例中，我们将使用有理函数（在每个多项式最多具有 1 的阶数）下从 0 到 1 的区间近似tan(2*x)函数。
In [26]: def error_function(a):
....: return (a[0] + a[1]*x)/(a[2] + a[3]*x) – np.tan(2*x)
In [27]: def jacobian(a):
....: numerator = a[0] + a[1]*x
....: denominator = a[2] + a[3]*x
....: return np.array( [ 1./denominator, x/denominator, \
....: -1.0*numerator/denominator**2, \
....: -1.0*x*numerator/denominator**2 ])
为了显示初始估计的依赖性，我们将尝试三种不同的选择：一种没有意义（所有系数为零），另一种为盲目标准选择（所有条目等于一个），另一种选择 承认tan(2*x)函数具有垂直渐近线这一事实。 我们将假装我们不知道确切的位置，并将其近似为0.78。 然后，我们的第三个初始估计表示一个简单的有理函数，在0.78处有一个渐近线。
显然，错误的初始估算并不能提供任何有用的信息：
In [28]: x1 = np.zeros((4,)); \
....: x2 = np.ones((4,)); \
....: x3 = np.array([1,0,0.78,-1])
In [29]: coeffs, success = leastsq(error_function, x1); \
....: numerator = np.poly1d(coeffs[1::-1]); \
....: denominator = np.poly1d(coeffs[:1:-1]); \
....: print numerator, denominator
0
0
In [30]: coeffs, success = leastsq(error_function, x1, \
....: Dfun=jacobian, col_deriv=True); \
....: numerator = np.poly1d(coeffs[1::-1]); \
....: denominator = np.poly1d(coeffs[:1:-1]); \
....: print numerator, denominator
0
0
使用x2作为初始猜测的这两个近似值都不令人满意：相应的误差很大，并且两个解都不具有从 0 到 1 的区间中的渐近线。
In [31]: coeffs, success = leastsq(error_function, x2); \
....: numerator = np.poly1d(coeffs[1::-1]); \
....: denominator = np.poly1d(coeffs[:1:-1]); \
....: print numerator, denominator; \
....: spla.norm(np.tan(2*x) - numerator(x) / denominator(x))
-9.729 x + 4.28
-1.293 x + 1.986
Out[31]: 220.59056436054016
In [32]: coeffs, success = leastsq(error_function, x2, \
....: Dfun=jacobian, col_deriv=True); \
....: numerator = np.poly1d(coeffs[1::-1]); \
....: denominator = np.poly1d(coeffs[:1:-1]); \
....: print numerator, denominator; \
....: spla.norm(np.tan(2*x) - numerator(x) / denominator(x))
-655.9 x + 288.5
-87.05 x + 133.8
Out[32]: 220.590564978504
使用 x3 作为初始猜测的近似值更接近目标函数，并且两者都具有可接受的渐近线。
In [33]: coeffs, success = leastsq(error_function, x3); \
....: numerator = np.poly1d(coeffs[1::-1]); \
....: denominator = np.poly1d(coeffs[:1:-1]); \
....: print numerator, denominator; \
....: spla.norm(np.tan(2*x) - numerator(x) / denominator(x))
0.01553 x + 0.02421
-0.07285 x + 0.05721
Out[33]: 2.185984698129936
In [34]: coeffs, success = leastsq(error_function, x3, \
....: Dfun=jacobian, col_deriv=True); \
....: numerator = np.poly1d(coeffs[1::-1]); \
....: denominator = np.poly1d(coeffs[:1:-1]); \
....: print numerator, denominator; \
....: spla.norm(np.tan(2*x) - numerator(x) / denominator(x))
17.17 x + 26.76
-80.52 x + 63.24
Out[34]: 2.1859846981334954
当然，我们可以做得更好，但是这些简单的示例现在就足够了。
如果我们希望输出更多信息来监视近似质量，可以将full_output选项设置为True来执行此操作：
In [35]: approximation_info = leastsq(error_function, x3, \
....: full_output=True)
In [36]: coeffs = approximation_info[0]; \
....: print coeffs
[ 0.02420694 0.01553346 0.0572128 -0.07284579]
In [37]: message = approximation_info[-2]; \
....: print message
Both actual and predicted relative reductions in the sum of squares
are at most 0.000000
In [38]: infodict = approximation_info[2]; \
....: print 'The algorithm performed \
....: {0:2d} iterations'.format(infodict['nfev'])
The algorithm performed 97 iterations
尽管从技术上讲，leastsq算法主要处理单变量函数的近似，但可以借助索引，散布，解包（使用特殊的*运算符）和稳定的和来处理多元函数。
提示
使用numpy实例方法sum（或使用numpy函数sum）的ndarray浮点数的总和远非稳定。 我们强烈建议不要将其用于较大数量的数字。 以下示例显示了一种不希望出现的情况，在这种情况下，我们尝试添加 4000 个值：
>>> arr=np.array([1,1e20,1,-1e20]*1000,dtype=np.float64)
>>> arr.sum() # The answer should be, of course, 2000
0.0
为了解决这种情况，我们使用稳定的总和。 在math模块中，为此目的提供了 Shewchuk 算法的实现：
>>> from math import fsum
>>> fsum(arr)
2000.0
有关 Shewchuk 算法的更多信息，以及使用浮点算术进行科学计算时应避免的其他常见陷阱，我们建议使用出色的指南*，David Goldberg 着每位计算机科学家应该了解的浮点算术*。 。 ACM 计算调查，1991 年。 23，第 5-48 页。
最好用一个例子来说明这个过程。 我们首先生成目标函数：32×32 尺寸的图像，除了三个具有不同位置，方差和高度的球形高斯函数外，还包含白噪声。 我们将所有这些值收集在一个 3×4 的数组中，这些数组称为值。 第一和第二列包含中心坐标的x和y值。 第三列包含高度，第四列包含差异。
In [39]: def sphericalGaussian(x0, y0, h, v):
....: return lambda x,y: h*np.exp(-0.5*((x-x0)**2+(y-y0)**2)/v)
....:
In [40]: domain = np.indices((32, 32)); \
....: values = np.random.randn(3,4); \
....: values[:,:2] += np.random.randint(1, 32, size=(3, 2)); \
....: values[:,2] += np.random.randint(1, 64, size=3); \
....: values[:,3] += np.random.randint(1, 16, size=3); \
....: print values
[[ 17.43247918 17.15301326 34.86691265 7.84836966]
[ 5.5450271 20.68753512 34.41364835 4.78337552]
[ 24.44909459 27.28360852 0.62186068 9.15251106]]
In [41]: img = np.random.randn(32,32)
In [42]: for k in xrange(3):
....: img += sphericalGaussian(*values[k])(*domain)
让我们假设我们不知道中心，高度和方差，并希望从目标图像img进行估计。 然后，我们创建一个误差函数来计算 3×4 数组a中压缩的 12 个系数。 请注意numpy函数 ravel 的作用和实例方法的重塑，以确保正确处理数据：
In [43]: from math import fsum
In [44]: def error_function(a):
....: a = a.reshape(3,4)
....: cx = a[:,0] # x-coords
....: cy = a[:,1] # y-coords
....: H = a[:,2] # heights
....: V = a[:,3] # variances
....: guess = np.zeros_like(img)
....: for i in xrange(guess.shape[0]):
....: for j in xrange(guess.shape[1]):
....: arr = H*np.exp(-0.5*((i-cx)**2+(j-cy)**2)/V)
....: guess[i,j] = fsum(arr)
....: return np.ravel(guess-img)
在保证成功的情况下，在这种情况下开始最小二乘法的过程需要一个近似的初步猜测。 对于此特定示例，我们将从名为values的数组中产生初始猜测：
In [45]: x0 = np.vectorize(int)(values); \
....: print x0
[[17 17 34 7]
[ 5 20 34 4]
[24 27 0 9]]
In [46]: leastsq(error_function, x0)
Out[46]:
(array([ 17.43346907, 17.14219682, 34.82077187, 7.85849653,
5.52511918, 20.68319748, 34.28559808, 4.8010449 ,
25.19824918, 24.02286107, 3.87170006, 0.5289382 ]),
1)
现在让我们通过以下过程直观地比较目标图像img及其最小化：
In [47]: coeffs, success = _; \
....: coeffs = coeffs.reshape(3,4)
In [48]: output = np.zeros_like(img)
In [49]: for k in xrange(3):
....: output += sphericalGaussian(*coeffs[k])(*domain)
In [50]: plt.figure(); \
....: plt.subplot(121); \
....: plt.imshow(img); \
....: plt.title('target image'); \
....: plt.subplot(122); \
....: plt.imshow(output); \
....: plt.title('computed approximation'); \
....: plt.show()
这给出了下图：
摘要
在本章中，我们探讨了近似理论领域中的两个基本问题：最小二乘意义上的插值和近似。 我们了解到有三种不同的方式可以解决 SciPy 中的这些问题：
- 一种程序模式，以
ndarrays的形式提供快速的数值解决方案。 - 提供
numpy功能的功能模式作为输出。 - 面向对象的模式，通过不同的类及其方法具有极大的灵活性。 当我们从解决方案中需要其他信息（例如，有关根，系数，结和错误的信息）或相关对象（例如，导数或反导数的表示形式）时，便使用此模式。
我们详细探讨了scipy.interpolate模块中编码的插值的所有不同实现，并且特别了解到与样条相关的实现是 Fortran 库FITPACK中多个例程的包装。
在最小二乘意义上的线性近似情况下，我们了解到，我们可以通过线性方程组（通过上一章中的技术）来实现解决方案，或者在样条近似的情况下，可以通过包装到 Fortran 例程中的方法来实现。 FITPACK库。 所有这些功能都在scipy.interpolate模块中编码。
对于最小二乘意义上的非线性近似，我们发现了scipy.optimize模块中编码的 Levenberg-Marquardt 迭代算法的两个变体。 这些依次从库MINPACK中调用 Fortran 例程LMDER和LMDIF。
在下一章中，我们将掌握差异化和集成的技术和应用。
