目录
- 一、参考
- 二、非线性最小二乘
- 三、GSL库中非线性最小二乘拟合部分
- 1. gsl_multifit_nlinear_parameters结构体
- 2. gsl_multifit_nlinear_type结构体
- 3. gsl_multifit_nlinear_alloc函数
- 4. gsl_multifit_nlinear_default_parameters函数
- 5. gsl_multifit_nlinear_init函数
- 6. gsl_multifit_nlinear_free函数
- 7. gsl_multifit_nlinear_name函数
- 8. gsl_multifit_nlinear_trs_name函数
- 9. gsl_multifit_nlinear_fdf结构体
- 10. gsl_multifit_nlinear_driver函数
- 11. gsl_multifit_nlinear_jac函数
- 12. gsl_multifit_nlinear_covar函数
- 四 、VS2019下的GSL库安装
- 五、非线性最小二乘的代码实现
一、参考
《数值最优化算法与理论》
数值最优化—概述
数值最优化—无约束问题的下降算法与线性搜索
数值最优化—无约束问题最速下降法和Newton法
数值最优化—无约束问题共轭梯度法
数值最优化—非线性方程组与最小二乘问题
二、非线性最小二乘
如下特殊形式的无约束问题：
m
i
n
f
(
x
)
=
1
2
∑
i
=
1
m
F
i
(
x
)
2
min f(x) = \frac 1 2 \sum ^m _{i=1} F_i(x)^2
minf(x)=21i=1∑mFi(x)2
其中，
F
i
:
R
n
→
R
(
i
=
1
,
2
,
⋅
⋅
⋅
,
m
)
F_i:R^n \to R(i=1,2,···,m)
Fi:Rn→R(i=1,2,⋅⋅⋅,m)连续可微。称该问题为最小二乘问题。非线性最小二乘问题包含非线性方程组作为其特殊情形，即
m
=
n
m=n
m=n。且该问题的最优解处的目标函数值为0。
当 F i ( i = 1 , 2 , ⋅ ⋅ ⋅ , m ) F_i(i=1,2,···,m) Fi(i=1,2,⋅⋅⋅,m)都是线性函数时，该问题称为线性最小二乘问题。
当问题中的至少有一个 F i F_i Fi是非线性函数，问题称为非线性最小二乘问题。
三、GSL库中非线性最小二乘拟合部分
本章描述多维非线性最小二乘拟合的函数。求解非线性最小二乘问题一般有两类算法，即行搜索法和信赖域法。GSL目前只实现信赖域法，并为用户提供对迭代中间步骤的完全访问。用户还能够调优一些参数，这些参数会影响算法的底层参数，有助于加快当前特定问题的收敛速度。GSL为非线性最小二乘拟合提供了两个独立的接口。第一个是为小到中等规模的问题设计的，第二个是为非常大的问题设计的，这些问题可能有也可能没有显著的稀疏结构。
头文件gsl_multifit_nlinear.h包含多维非线性拟合函数的原型，以及与小型到中型系统相关的声明。
头文件gsl_multilarge_nlinear.h包含多维非线性拟合函数的原型，以及与大型系统相关的声明。
1. gsl_multifit_nlinear_parameters结构体
用户可以在分配内存时，对迭代的几乎所有方面进行调优。定义在gsl_multifit_nlinear.h，用户可以修改gsl_multifit_nlinear_parameters结构体中参数，其定义如下：
typedef struct
{
const gsl_multifit_nlinear_trs *trs; /* trust region subproblem method */
const gsl_multifit_nlinear_scale *scale; /* scaling method */
const gsl_multifit_nlinear_solver *solver; /* solver method */
gsl_multifit_nlinear_fdtype fdtype; /* finite difference method */
double factor_up; /* factor for increasing trust radius */
double factor_down; /* factor for decreasing trust radius */
double avmax; /* max allowed |a|/|v| */
double h_df; /* step size for finite difference Jacobian */
double h_fvv; /* step size for finite difference fvv */
} gsl_multifit_nlinear_parameters;
① const gsl_multifit_nlinear_trs *trs;：决定了用于解决信赖域子问题的方法，可以从以下选项中选择：
/* trust region subproblem methods */
GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_lm; //Levenberg-Marquardt算法
GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_lmaccel; //带有测地加速度的Levenberg-Marquardt算法
GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_dogleg; //狗腿算法
GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_ddogleg; //双狗腿算法
GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_subspace2D; //二维子空间算法
② const gsl_multifit_nlinear_scale *scale;：决定了对角缩放矩阵D，可以从以下选项中选择：
/* scaling matrix strategies */
GSL_VAR const gsl_multifit_nlinear_scale * gsl_multifit_nlinear_scale_levenberg; //由Levenberg提出的阻尼策略。
GSL_VAR const gsl_multifit_nlinear_scale * gsl_multifit_nlinear_scale_marquardt; //由Marquartdt提出的阻尼策略。
GSL_VAR const gsl_multifit_nlinear_scale * gsl_multifit_nlinear_scale_more; //由More提出的阻尼策略，本库默认选择。
③ const gsl_multifit_nlinear_solver *solver;：决定了如何求解信赖域子问题中线性最小二乘系统，可以从以下选项中选择：
/* linear solvers */
GSL_VAR const gsl_multifit_nlinear_solver * gsl_multifit_nlinear_solver_cholesky; //使用Cholesky分解来求解系统。
GSL_VAR const gsl_multifit_nlinear_solver * gsl_multifit_nlinear_solver_qr; //使用雅可比矩阵的QR分解来求解系统。
GSL_VAR const gsl_multifit_nlinear_solver * gsl_multifit_nlinear_solver_svd; //使用雅可比矩阵的奇异值分解来求解系统。
④ gsl_multifit_nlinear_fdtype fdtype;：指定了近似雅克比矩阵时使用前向差分还是中心差分。这只在解析雅克比矩阵时用，不提供给求解器使用。可以从以下选项中选择：
GSL_MULTIFIT_NLINEAR_FWDIFF, //前向有限差分
GSL_MULTIFIT_NLINEAR_CTRDIFF //中心有限差分
⑤ double factor_up; ：当接受一步时，信赖域半径将增加这个因子。默认值3。
⑥ double factor_down; ：当拒绝一步时，信赖域半径将减少这个因子。默认值2。
⑦ double avmax; ：用测地加速度求解非线性最小二乘时，设置加速度项与速度项的比值的阈值。默认0.75。
⑧ double h_df; ：用有限差分逼近雅克比矩阵的步长。默认
ϵ
\sqrt {\epsilon}
ϵ。
⑨ double h_fvv; ：用测地加速度时，并使用有限差分时估计第二个方向导数时的步长。默认0.02。
2. gsl_multifit_nlinear_type结构体
指定了用于解决非线性最小二乘问题的算法类型。目前GSL库只实现了这一种算法：
/* top-level algorithms */
GSL_VAR const gsl_multifit_nlinear_type * gsl_multifit_nlinear_trust; //信赖域算法
3. gsl_multifit_nlinear_alloc函数
返回一个指针，指向工作空间。
gsl_multifit_nlinear_type 和 gsl_multifit_nlinear_parameters结构体如上。
n：拟合的数据点个数。
p：待拟合的参数个数。
GSL_FUN gsl_multifit_nlinear_workspace *
gsl_multifit_nlinear_alloc (const gsl_multifit_nlinear_type * T,
const gsl_multifit_nlinear_parameters * params,
size_t n, size_t p);
4. gsl_multifit_nlinear_default_parameters函数
返回一组用于解决非线性最小二乘问题推荐的默认参数。用户可以在这之后对每个参数进行优化，以提高针对特定问题的性能。
GSL_FUN gsl_multifit_nlinear_parameters gsl_multifit_nlinear_default_parameters(void);
5. gsl_multifit_nlinear_init函数
使用fdf 和 初始参数x，初始化现有的工作空间w。
GSL_FUN int
gsl_multifit_nlinear_init (const gsl_vector * x,
gsl_multifit_nlinear_fdf * fdf,
gsl_multifit_nlinear_workspace * w);
6. gsl_multifit_nlinear_free函数
释放工作空间w的内存。
GSL_FUN void gsl_multifit_nlinear_free (gsl_multifit_nlinear_workspace * w);
7. gsl_multifit_nlinear_name函数
返回工作空间w的求解器名称。调试用。
GSL_FUN const char *
gsl_multifit_nlinear_name (const gsl_multifit_nlinear_workspace * w);
8. gsl_multifit_nlinear_trs_name函数
返回工作空间w的求解信赖域子问题的算法名称。调试用。
GSL_FUN const char *
gsl_multifit_nlinear_trs_name (const gsl_multifit_nlinear_workspace * w);
9. gsl_multifit_nlinear_fdf结构体
typedef struct
{
int (* f) (const gsl_vector * x, void * params, gsl_vector * f); //待拟合的函数
int (* df) (const gsl_vector * x, void * params, gsl_matrix * df); //设置雅克比矩阵。如果设为NULL，雅克比矩阵将使用函数f的有限差分近似在内部计算。
int (* fvv) (const gsl_vector * x, const gsl_vector * v, void * params,
gsl_vector * fvv); //测地加速度相关
size_t n; /* number of functions （数据点个数）*/
size_t p; /* number of independent variables （待拟合的函数参数个数）*/
void * params; /* user parameters */
size_t nevalf; /* number of function evaluations （无需用户设置） */
size_t nevaldf; /* number of Jacobian evaluations（无需用户设置） */
size_t nevalfvv; /* number of fvv evaluations（无需用户设置） */
} gsl_multifit_nlinear_fdf;
10. gsl_multifit_nlinear_driver函数
使用工作空间w下的求解器进行最多maxiter次的迭代。每次迭代后，以xtol,gtol,ftol来测试系统的收敛性。
用户可以提供一个回调函数callback，在每次迭代后，用户可以保存或打印相关的值。当callback设置为NULL时，代表不回调。
GSL_FUN int
gsl_multifit_nlinear_driver (const size_t maxiter,
const double xtol,
const double gtol,
const double ftol,
void (*callback)(const size_t iter, void *params,
const gsl_multifit_nlinear_workspace *w),
void *callback_params,
int *info,
gsl_multifit_nlinear_workspace * w);
11. gsl_multifit_nlinear_jac函数
返回一个指针，指向工作空间w的雅克比矩阵。
GSL_FUN gsl_matrix *
gsl_multifit_nlinear_jac (const gsl_multifit_nlinear_workspace * w);
12. gsl_multifit_nlinear_covar函数
利用雅克比矩阵J计算最佳拟合参数的协方差矩阵，并将其存储在covar中。
GSL_FUN int
gsl_multifit_nlinear_covar (const gsl_matrix * J, const double epsrel,
gsl_matrix * covar);
四 、VS2019下的GSL库安装
在VS2019打开Nuget包管理器。搜索gsl。安装gsl-msvc-x64包。如下：
VS2019同样选择x64平台开始调试。在工程中直接include gsl的头文件即可开始使用。
如果在其他工程中使用gsl库，可以在电脑中找到VS2019的Nuget的GSL库的下载位置，将lib和头文件拷贝到需要的工程目录下。
五、非线性最小二乘的代码实现
main.cpp：
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlinear.h>
#define N 100 /* number of data points to fit */
#define TMAX (3.0) /* time variable in [0,TMAX] */
struct data {
size_t n;
double * t;
double * y;
};
int
expb_f (const gsl_vector * x, void *data,
gsl_vector * f)
{
size_t n = ((struct data *)data)->n;
double *t = ((struct data *)data)->t;
double *y = ((struct data *)data)->y;
double A = gsl_vector_get (x, 0);
double lambda = gsl_vector_get (x, 1);
double b = gsl_vector_get (x, 2);
size_t i;
for (i = 0; i < n; i++)
{
/* Model Yi = A * exp(-lambda * t_i) + b */
double Yi = A * exp (-lambda * t[i]) + b;
gsl_vector_set (f, i, Yi - y[i]);
}
return GSL_SUCCESS;
}
int
expb_df (const gsl_vector * x, void *data,
gsl_matrix * J)
{
size_t n = ((struct data *)data)->n;
double *t = ((struct data *)data)->t;
double A = gsl_vector_get (x, 0);
double lambda = gsl_vector_get (x, 1);
size_t i;
for (i = 0; i < n; i++)
{
/* Jacobian matrix J(i,j) = dfi / dxj, */
/* where fi = (Yi - yi)/sigma[i], */
/* Yi = A * exp(-lambda * t_i) + b */
/* and the xj are the parameters (A,lambda,b) */
double e = exp(-lambda * t[i]);
gsl_matrix_set (J, i, 0, e);
gsl_matrix_set (J, i, 1, -t[i] * A * e);
gsl_matrix_set (J, i, 2, 1.0);
}
return GSL_SUCCESS;
}
void
callback(const size_t iter, void *params,
const gsl_multifit_nlinear_workspace *w)
{
gsl_vector *f = gsl_multifit_nlinear_residual(w);
gsl_vector *x = gsl_multifit_nlinear_position(w);
double rcond;
/* compute reciprocal condition number of J(x) */
gsl_multifit_nlinear_rcond(&rcond, w);
fprintf(stderr, "iter %2zu: A = %.4f, lambda = %.4f, b = %.4f, cond(J) = %8.4f, |f(x)| = %.4f\n",
iter,
gsl_vector_get(x, 0),
gsl_vector_get(x, 1),
gsl_vector_get(x, 2),
1.0 / rcond,
gsl_blas_dnrm2(f));
}
int
main (void)
{
const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
gsl_multifit_nlinear_workspace *w;
gsl_multifit_nlinear_fdf fdf;
gsl_multifit_nlinear_parameters fdf_params =
gsl_multifit_nlinear_default_parameters();
const size_t n = N;
const size_t p = 3;
gsl_vector *f;
gsl_matrix *J;
gsl_matrix *covar = gsl_matrix_alloc (p, p);
double t[N], y[N], weights[N];
struct data d = { n, t, y };
double x_init[3] = { 1.0, 1.0, 0.0 }; /* starting values */
gsl_vector_view x = gsl_vector_view_array (x_init, p);
gsl_vector_view wts = gsl_vector_view_array(weights, n);
gsl_rng * r;
double chisq, chisq0;
int status, info;
size_t i;
const double xtol = 1e-8;
const double gtol = 1e-8;
const double ftol = 0.0;
gsl_rng_env_setup();
r = gsl_rng_alloc(gsl_rng_default);
/* define the function to be minimized */
fdf.f = expb_f;
fdf.df = expb_df; /* set to NULL for finite-difference Jacobian */
fdf.fvv = NULL; /* not using geodesic acceleration */
fdf.n = n;
fdf.p = p;
fdf.params = &d;
/* this is the data to be fitted */
for (i = 0; i < n; i++)
{
double ti = i * TMAX / (n - 1.0);
double yi = 1.0 + 5 * exp (-1.5 * ti);
double si = 0.1 * yi;
double dy = gsl_ran_gaussian(r, si);
t[i] = ti;
y[i] = yi + dy;
weights[i] = 1.0 / (si * si);
printf ("data: %g %g %g\n", ti, y[i], si);
};
/* allocate workspace with default parameters */
w = gsl_multifit_nlinear_alloc (T, &fdf_params, n, p);
/* initialize solver with starting point and weights */
gsl_multifit_nlinear_winit (&x.vector, &wts.vector, &fdf, w);
/* compute initial cost function */
f = gsl_multifit_nlinear_residual(w);
gsl_blas_ddot(f, f, &chisq0);
/* solve the system with a maximum of 100 iterations */
status = gsl_multifit_nlinear_driver(100, xtol, gtol, ftol,
callback, NULL, &info, w);
/* compute covariance of best fit parameters */
J = gsl_multifit_nlinear_jac(w);
gsl_multifit_nlinear_covar (J, 0.0, covar);
/* compute final cost */
gsl_blas_ddot(f, f, &chisq);
#define FIT(i) gsl_vector_get(w->x, i)
#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
fprintf(stderr, "summary from method '%s/%s'\n",
gsl_multifit_nlinear_name(w),
gsl_multifit_nlinear_trs_name(w));
fprintf(stderr, "number of iterations: %zu\n",
gsl_multifit_nlinear_niter(w));
fprintf(stderr, "function evaluations: %zu\n", fdf.nevalf);
fprintf(stderr, "Jacobian evaluations: %zu\n", fdf.nevaldf);
fprintf(stderr, "reason for stopping: %s\n",
(info == 1) ? "small step size" : "small gradient");
fprintf(stderr, "initial |f(x)| = %f\n", sqrt(chisq0));
fprintf(stderr, "final |f(x)| = %f\n", sqrt(chisq));
{
double dof = n - p;
double c = GSL_MAX_DBL(1, sqrt(chisq / dof));
fprintf(stderr, "chisq/dof = %g\n", chisq / dof);
fprintf (stderr, "A = %.5f +/- %.5f\n", FIT(0), c*ERR(0));
fprintf (stderr, "lambda = %.5f +/- %.5f\n", FIT(1), c*ERR(1));
fprintf (stderr, "b = %.5f +/- %.5f\n", FIT(2), c*ERR(2));
}
fprintf (stderr, "status = %s\n", gsl_strerror (status));
gsl_multifit_nlinear_free (w);
gsl_matrix_free (covar);
gsl_rng_free (r);
return 0;
}
运行结果如图：
程序中：
① 待拟合的曲线函数为 Y = A e − λ t + b Y = A e^{-\lambda t} + b Y=Ae−λt+b。
② 输入了100个加了高斯误差的数据点 t ，并且每个数据点具有权值。
③ 设置了误差为 1 0 − 8 10^{-8} 10−8。最大迭代次数100次。
④ 设置了初始参数为A = 1.0。 λ \lambda λ=1.0。b=1.0。
⑤ 最后求解出的曲线函数的参数为A = 4.79653。 λ \lambda λ=1.43937。b=1.00368。
拟合效果绘制曲线如下图：
84
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
