利用matlab实现非线性拟合（三维、高维、参数方程）

2021.06 更新。补充了非线性拟合中，关于最小二乘定义的问题，放在了最后一章。 2021.12更新。应评论区建议，补充了3.3章绘图用的代码。

本文首发于 matlab爱好者 微信公众号，欢迎关注。

惯例声明：本人没有相关的工程应用经验，只是纯粹对相关算法感兴趣才写此博客。所以如果有错误，欢迎在评论区指正，不胜感激。本文主要关注于算法的实现，对于实际应用等问题本人没有任何经验，所以也不再涉及。

本人在学习matlab匿名函数时，作为练习有感而发，写下下面内容，非专业，难免有误。

一般而言，通过已有的数据点去推导其它数据点，常见的方法有插值和拟合。插值适用性较广，尤其是线性插值或样条插值已被广泛的应用。但是通过已知的函数去拟合数据，是连接理论与实验重要的桥梁，这一点是插值无法替代的。

日常学习工作中，经常会遇到下面这种问题：想要用某个具体函数去拟合自己的数据，明明知道这个函数的具体形式，却不知道其中的参数怎么选取。本文就简单介绍一下matlab环境下，如何进行非线性拟合。

多项式拟合就是利用下面形式的方程去拟合数据： y = a + b x + c x 2 + d x 3 + . . . y=a +bx+cx^2+dx^3+... y=a+bx+cx2+dx3+... matlab中可以用polyfit()函数进行多项式拟合。下面举一个小例子：

对于已有的数据点，我们采用4阶多项式拟合。其中已知函数的的表达式为

在此基础上添加了一些噪声点。拟合曲线依然采用4阶进行拟合，结果如下。

可以看到拟合曲线与理论曲线基本一致，说明这种方法能够较好的拟合出原始数据的趋势。源代码如下：

对于多项式拟合，不是阶数越多越好。比如还是这个上面这个例子，阶数越多，曲线越能够穿过每一个点，但是曲线的形状与理论曲线偏离越大。所以选择最适合的才是最好的。

线性拟合就是能够把拟合函数写成下面这种形式的： y = p 1 f 1 ( x ) + p 2 f 2 ( x ) + p 3 f 3 ( x ) + . . . y=p_1 f_1(x) +p_2 f_2(x)+p_3 f_{3}(x)+... y=p1​f1​(x)+p2​f2​(x)+p3​f3​(x)+...

其中f(x)是关于x的函数，其表达式是已知的。p是常数系数，这个是未知的。

对于这种形式的拟合，matlab内部有一个及其强悍的函数，可以自动输出p的解，并且满足最小二乘。这个函数就是\。没错，就是斜杠。这个符号通常用于求解方程AX=B的情况，我们用X=A\B可以求出未知数X。我们利用当A行和列不等时，输出X的最小二乘这个特性，就可以求出相应的最佳拟合。

还是举个例子

虽然看上去很非线性，但是x和y都是已知的，a、b、c是未知的，而且是线性的，所以可以被非常简单的拟合出来。假设a=2.5，b=0.5，c=-1，加入随机扰动。拟合的最终效果为：

最终得到的拟合参数为：a=2.47，b=0.47，c=-0.66。考虑到随机噪声的影响，与原始数据相差不大，源代码如下：

事实上，其实常用的拟合方式中，有很多都是线性拟合，比如多项式拟合，傅里叶拟合等。对于多项式拟合，只需要把拟合的部分替换成x，x.^2，x.^3这种形式。对于傅里叶级数，只需要把拟合的部分替换为sin(x)，sin(2*x)，sin(3*x)这种形式就行。

下面展示一个利用线性拟合，进行不同频率的三角波级数拟合正弦函数的例子：

前面介绍了线性的拟合方法。如果一个非线性方程，可以化为上面线性方程中公式中给出的样子，那么我们也可以套用线性的方法去求解。

比如下面的方程： y = a ∗ exp ⁡ ( − b x ) y=a*\exp(-bx) y=a∗exp(−bx) 经过取对数变换，可以变为下面等效的线性形式： lg ⁡ ( y ) = lg ⁡ ( a ) + b ∗ ( − x ) \lg(y)=\lg(a)+b*(-x) lg(y)=lg(a)+b∗(−x) 这样求出来之后，再反变换回去，就可以得到原方程的系数了。 代码如下：

matlab自带了一个Curve Fitting App，它在matlab集成的App里面。 界面左边输入数据，右侧选择方法。 界面里常用的拟合方式都有，而且直接展示拟合效果，非常方便。非常适合鼠标直接拖拖拽拽点点点的操作方式。

matlab中，fit()函数是一个比较通用的用于函数拟合的函数。它的优点就是非常的全面，可以用于各种种类的拟合。上面的App里，很多拟合种类都是间接调用了fit函数来实现的拟合。

对于非线性拟合，可以使用fit()函数中的Nonlinear Least Squares方法。其大概原理为，首先确定一个初始的点，计算该点的雅可比矩阵J，来近似线性化的判断该点周围的趋势，并将这个点向更小的方向移动。

因此，这个方法的一个缺点在于，对于初始点的选取非常敏感，最终结果只能在初始点附近的局部最小值点上，而不能保证全局最小值。

相比于前面的fit()函数，nlinfit()函数是matlab专门的非线性拟合函数。对于非稳健估计，采用的是Levenberg-Marquardt（LM）方法，也叫阻尼最小二乘法。对于稳健估计，采用的是Iteratively Reweighted Least Squares方法，也就是在Least Squares基础上，对每一个拟合点的权重进行调整的一种方法。这两者方法也都是基于雅克比矩阵的方法。

lsqnonlin()也是matlab中自带的一个非线性拟合函数。它给出了两种计算非线性拟合的方法，一种是比较经典之前介绍过的LM方法，一种是信赖域方法。信赖域法（trust region reflective）是通过Hessian 矩阵，逐步试探邻域内的最小化，来求解问题的。相比较之前的那些雅克比相关的方法，信赖域法会占用更多内存和速度，所以适用于中小规模的矩阵。

lsqcurvefit()函数和lsqnonlin()内容上相似，只是引用格式上有所不同。

这也是一个求解非线性方程的函数，可以求解方程组或者矩阵形式，功能非常强大。默认的算法为trust-region-dogleg，俗称狗腿法，属于信赖域法。这里用到的功能比较基础，所以也不过多介绍了。

说了那么多，发现逐渐从如何非线性拟合，陷入到了最优化的深坑里。而且前面的那么多方法，很多都解决不了陷入到局部最优解的问题。实际上，这种问题如果进入了最优化领域，很多智能算法也可以被考虑进来。所以我也把粒子群PSO算法加入到了里面，尝试将结果收敛到全局最优解。

第一个例子是 y=a.*exp(-b\*(x-c).^2)+d，一个简单的高斯函数形式的非线性方程，其参数给定为：

第二个例子是一个指数增长的正弦函数，在很多线性系统中都可以测量到这种信号。函数的形式为：

y=a*x+b*sin(c*x).*exp(d*x)+e 。其给定的参数为：

这个函数的拟合具有一定难度，拟合过程中会遇到非常多的局部解。

由于初始点比较随机，所以最终结果每次都会不一样。 其中前面的几种方法对于初始值的敏感度比较高，如果初始值选的比较接近原始解，也是可以得到较好的结果。其中nlinfit函数经常会报错，容错率较低。而PSO算法经常能够收敛到最优解（虽然不是每次都可以，偶尔也会陷入局部解）。

代码如下：

下图展示了PSO求解过程中逐渐收敛到全局最优解的过程。

之前的文章中的数据具有一 一对应的特点，所以严格来讲并不是普遍的二维拟合。对于一些复杂的二维函数，比如椭圆，可能原本的拟合就会失效。

对于这种一般性质的方程或方程组，将原本输入方程整理为f(x,y,z,…)=0的形式。衡量拟合程度的优化函数，就直接取函数f(xi,yi,zi,…)的值即可。

下面演示最终的两个例子：

第一个是三维直线，采用两平面式描述。

总共2个方程，维度为3维，第一个方程有3个参数，第二个方程也有3个参数。离散点已知的条件下，三维直线的两平面表达式不唯一。

最终拟合效果如下：

第二个是二维椭圆，方程为：

总共1个方程，维度为2维。方程共有5个参数。

最终拟合效果如下：

高维参数方程的拟合比较困难。因为原本方程中x、y、z的坐标点都是已知的。但是参数方程中，x、y、z的坐标点已知，但是与参数u、v往往未知。所以相当于原本的方程中引入了额外的未知信息。

但是基本思路和普通方程是一样的。离散点距离假想方程的距离，需要用该点距方程上每个点的距离中的最小值，来近似判断。

依然是展示两个例子。

第一个是计算三维李萨如图形。参数方程为：

方程为三维参数方程，只有一个参数u。第一个方程有1个未知量A，第二个方程有1个未知量B，第三方程有1个未知量C。

最终拟合效果如下。由于李萨如图形中，只要频率的比例固定，图案就会固定。所以最终ABC的值不唯一，但是它们的比例肯定唯一。

第二个例子是一个三维旋转曲面。参数方程为：

方程为三维参数方程，有2个参数u、v。第一个方程有2个未知量AB，第二个方程有2个未知量CD，第三方程有0个未知量。

最终拟合效果如下：

这两个例子的演示代码如下。这里参数方程的Dis_P()由于频繁的计算点与点之间的距离，所以运算速度比第一章单纯函数的Dis()较慢。

主要函数就是Dis和Dis_P，准备工作就是把方程和离散点整理成可以输入的形式。优化用到的函数就是PSO（particleswarm），需要更改未知参数数量和范围就可以。

这一部分的拟合优度或者置信区间就没办法拿出来说了，毕竟有试凑的嫌疑。适合单纯的验证。

另外这里和最小二乘求出来的结果会有所不同。这主要是定义的问题。其中最小二乘指的是距离y方向上最小的二乘距离。而这一章节中定义的最小二乘，是指的点到拟合曲线距离的最小二乘（有点类似于主成分分析）。虽然一般情况下两者并不会有太大差别，但是如果有误差的话，还请注意这一点，

如下图，这里是一个椭圆分布的离散点。 如果将方程的形式定义为y=ax+b，要求距离y最小，则最终结果为红色的线。如果将方程定义为ax+by+c=0，要求距离直线的距离最小，则最终拟合结果为黄色的线。这里无所谓对错，只是要求的最小定义不同。

单独在后面补充了这么一个章节，希望在应用中注意。

上面的图，绘图代码如下：

信赖域法 https://zhuanlan.zhihu.com/p/205555114