Data and Model Fitting
Mathwrist C++ API Series

(January 20, 2023)

Abstract

This document details the C++ programming interface of data and model fitting features in Mathwrist’s Numerical Programming Library (NPL). Please refer to our “Data and Model Fitting” white paper for technical overview.

1 Linear Least Square

Assume a $n$ dimensional linear model $y=\mathbf{x}^{T}\beta:\mathbb{R}^{n}\rightarrow\mathbb{R}$ parameterized by coefficient vector $\beta$ . Given $m$ number of observations $\left(y_{i},\mathbf{x}_{i}\right),i=0,\cdots,m-1$ , we provides two least square based methods to recover parameters $\beta$ in “linear_ls.h” header file.

1.1 Generalized Ridge Regression

Function linear_fit::grr() under namesapce linear_fit solves a generalized ridge regression problem as the following,

\arg\min_{\beta}\left(\mathbf{y}-\mathbf{X}\beta\right)^{T}\left(\mathbf{y}-% \mathbf{X}\beta\right)+\lambda\beta^{T}\mathbf{\Omega}\beta

(1)

1 void linear_fit::grr(

2 const Matrix& X,

3 const Matrix& y,

4 const Matrix& Omega,

5 double lambda,

6 Matrix& beta);

Parameter Description:

1.

X is the $m\times n$ matrix $\mathbf{X}$ in equation (1). The $i$ -th row of $\mathbf{X}$ holds the $i$ -th observation of n-dimensional independent variables.
2.

y is the $m\times 1$ vector $\mathbf{y}$ in equation (1). The $i$ -th element of $\mathbf{y}$ is the observed linear model value.
3.

Omega is a $n\times n$ regularization matrix $\mathbf{\Omega}$ in equation (1). It needs be at least positive definite.
4.

lambda is the penalty factor $\lambda$ in equation (1). If a positive value $\lambda>0$ is supplied, this penalty factor will be used to solve (1). Otherwise, we will use a Generalized Cross Validation (GCV) method to compute an optimal $\lambda^{*}$ based on the noisy level of the data.
5.

beta will holds the n-d linear model parameters $\beta$ on funciton return.

1.2 Constrained Least Square

Function linear_fit::lsq under namespace linear_fit minimizes the sum of square residuals (SSR) subject to general linear constraints and simple bounds.

$\displaystyle\arg\min_{\beta}\left(\mathbf{y}-\mathbf{X}\beta\right)^{T}\left(% \mathbf{y}-\mathbf{X}\beta\right)$	s.t.	(2)
$\displaystyle\mathbf{b}_{l}\leq\mathbf{A}\beta\leq\mathbf{b}_{u}\mbox{ and }$		(3)
$\displaystyle\mathbf{l}\leq\beta\leq\mathbf{u}$		(4)

1 void linear_fit::lsq(

2 const Matrix& X,

3 const Matrix& y,

4 const LinearProg::Constraints& c,

5 Matrix& beta);

Parameter Description:

1.

X is the $m\times n$ matrix $\mathbf{X}$ in equation (2). The $i$ -th row of $\mathbf{X}$ holds the $i$ -th observation of n-dimensional independent variables.
2.

y is the $m\times 1$ vector $\mathbf{y}$ in equation (2). The $i$ -th element of $\mathbf{y}$ is the observed linear model value.
3.

c is the linear constraint specification that includes general linear constraints in equation (LABEL:eqn:linear-lsq-general-constraints) and simple bounds in equation (LABEL:eqn:linear-lsq-simple-bound). Please refer to our Linear Programming (LP) API for details of linear constraint specification.
4.

beta will holds the n-d linear model parameters $\beta$ on funciton return.

2 Nonlinear Least Square

Let $y=h(\mathbf{x};\beta):\mathbb{R}^{n}\rightarrow\mathbb{R}$ be a $n$ -dimensional nonlinear function with model parameter $\beta$ . Nonlinear least square fit estimates parameter $\beta$ given $m$ number of observations $(y_{i},\mathbf{x}_{i})$ , $i=0,\cdots,m-1$ by minimizing the $l_{2}$ norm of residual vector,

\arg\min_{\beta}\psi(\beta)=\frac{1}{2}\|\mathbf{r}(\beta)\|_{2}^{2}

(5)

, where the $i$ -th element $r_{i}(\beta)=h(\mathbf{x}_{i};\beta)-y_{i}$ .

In header file “non_linear_ls.h”, we provide three least square based methods, all derived from the NonlinearLS base class, to calibrate model parameters $\beta$ . Each method is represented by a solver class designed with the following inheritance hierarchy,

1 class NonlinearLS : public IterativeMethod

2 {

3 };

4 class NLS_GaussNewton : public NonlinearLS

5 {

6 };

7 class NLS_Levenberg : public NonlinearLS

8 {

9 };

10 class NLS_General : public NonlinearLS

11 {

12 };

2.1 Termination Control

The exact termination behavior depends on the actual fitting methodology implemented in a derived least square solver class. Please refer to our “Data and Model Fitting” white paper for the technical details of termination conditions. In general, a nonlinear least square solver terminates at a major iteration whenever any of the following conditions is satisfied,

1.

the residual satisfies $|\mathbf{r}(\beta)|_{\infty}<\epsilon_{r}$ , where $\epsilon_{r}$ is the residual tolerance.
2.

the max number of iterations $M$ is reached.
3.

both the stationary and convergence conditions are satisfied with $\epsilon_{c}$ , where $\epsilon_{c}$ is the convergence tolerance.

The max number of iterations $M$ and convergence tolerance $\epsilon_{c}$ are control functions defined in IterativeMethod base class. Please refer to our “1-d and n-d Functions” API documentation for the details of IterativeMethod.

The residual tolerance $\epsilon_{r}$ control functions are public member functions of the NonlinearLS base class.

1 void set_residual_tolerance(double tol);

2 double residual_tolerance() const;

Description:

1.

Set the tolerance $\epsilon_{r}$ to be stored as a class data member of NonlinearLS.
2.

Return the tolerance $\epsilon_{r}$ that is used to test termination condition.

2.2 Max Step Change

1 void set_max_step(double magnitude);

2 double max_step() const;

Description:

Due to the nonlinear nature of function $h(\mathbf{x};\beta)$ , we want to restrict the max amount of step change $\Delta$ at each iteration of the least square solver, details as the following:

•

In modified Gauss-Newton method, $\Delta$ determines the max step change in a line search method.
•

In Levenberg-Marquardt method, $\Delta$ determines the max trust region radius.
•

In regularized and constrained method, $\Delta$ determines the max step change in a line search method.

The following two functions are public member functions of the NonlinearLS base class.

1.

Set the max step change $\Delta$ passed from parameter magnitude.
2.

Return the max step change $\Delta$ used by the solver.

2.3 Solver Operator

1 virtual void operator()(const VtrValueFunctionND& residuals,

2 Matrix& beta) = 0;

Description:

This is the public interface function declared in NonlinearLS base class and implemented by concrete nonlinera least square derived classes.

1.

Parameter residuals is the residual function $\mathbf{r}(\beta)$ .
2.

Parameter beta holds the initial guess of model parameter $\beta$ on call and holds the calibrated parameters $\beta$ on return.

2.4 Modified Gauss-Newton

1 NLS_GaussNewton();

2 void set_bump_size(double dx);

3 double bump_size() const;

Description:

These functions are public function of the derived class NLS_GaussNewton that uses a modified Gauss-Newton method to solve the nonlinear least square problem.

1.

Constructor, create a solver object.
2.

Set the bump size that is used to compute the second order Hessian term $\mathbf{Q}(\beta)$ by finite (center) difference. Details in white paper.
3.

Return the bump size to compute $\mathbf{Q}(\beta)$ by finite difference.

2.5 Levenberg-Marquardt

1 NLS_Levenberg(double max_radius, size_t max_lambda_iters);

2 void set_max_lambda_iterations(size_t n);

3 size_t max_lambda_iterations() const;

Description:

These functions are public function of the derived class NLS_Levenberg that uses the Levenberg-Marquardt method to solve the nonlinear least square problem.

1.

Constructor, create a solver object. Parameter max_radius determines the max trust region radius. Parameter max_lambda_iters determines the max number of iterations in root finding $\lambda$ in trust region method. Details in white paper.
2.

Set the max number of iterations allowed to root find $\lambda$ in trust region method.
3.

Return the max number of iterations allowed to root find $\lambda$ in trust region method.

2.6 Regularization and Constraints

1 NLS_General(const Matrix* _O=0);

2 NLS_General(const LinearProg::Constraints& c, const Matrix* _O=0);

3 void set_reg_penalty(double lambda=0);

Description:

These functions are public function of the derived class NLS_General that adds a regularization term, and optionally additional constraints, to solve the nonlinear least square problem. The regularization matrix $\mathbf{\Omega}$ needs be at least positive semi-definite.

\arg\min_{\beta}\psi(\beta)=\frac{1}{2}\|\mathbf{r}(\beta)\|_{2}^{2}+\lambda% \beta^{T}\mathbf{\Omega}\beta\mbox{ s.t. }

(6)

	$\displaystyle\mathbf{b}_{l}\leq\mathbf{A}\beta\leq\mathbf{b}_{u}\mbox{ and }$		(7)
	$\displaystyle\mathbf{l}\leq\beta\leq\mathbf{u}$		(8)

1.

Constructor, create a solver object. Parameter _O, if supplied, is the regularization matrix $\mathbf{\Omega}$ in equation (6). If a null pointer is supplied, $\mathbf{\Omega}$ is taken as an identity matrix.
2.

Constructor, create a solver object. Parameter c is the linear constraint specification in equation (LABEL:eqn:nlsq-general-constraint). Please refer to our LP API for constraint specification details. Parameter _O is a pointer to the regularization matrix $\mathbf{\Omega}$ as the first constructor.
3.

Store the regularization penalty $\lambda$ in equation (6) as a data member of the solver. If $\lambda\leq 0$ is supplied, the algorithm will use the Generalized Cross Validation (GCV) method to compute an optimal penalty factor $\lambda^{*}$ .

3 Curve Fitting

We represent a 1-d curve function $f(x;\theta)=\phi(x)^{T}\theta$ as the linear combination of a basis vector $\phi(x)^{T}=(\phi_{0}(x),\cdots,\phi_{n-1}(x))$ and coefficient vector $\theta^{T}=(\theta_{0}(x),\cdots,\theta_{n-1}(x))$ .

Mathwrist NPL offers two choices of basis functions $\phi(x)$ . Class SmoothSplineCurve uses B-spline as basis. Class SmoothChebyshevCurve uses Chebyshev polynomial of the first kind as basis. They are both derived from SmoothCurveFitting, which in turn inherits from the top level LinearModelFitting base class.

1 class LinearModelFitting

2 {

3 };

4 class SmoothCurveFitting : public LinearModelFitting

5 {

6 };

7 class SmoothSplineCurve : public SmoothCurveFitting

8 {

9 };

10 class SmoothChebyshevCurve : public SmoothCurveFitting

11 {

12 };

The relevant C++ header files are:

1.

linear_model_fitting.h
2.

smooth_curve_fitting.h
3.

smooth_spline_curve.h
4.

smooth_chebyshev_curve.h

3.1 Smooth Curve Fitting

This section documents the public interface of the SmoothCurveFitting base class. Derived class SmoothSplineCurve and SmoothChebyshevCurve share the same API as the base class.

3.1.1 Curve Domain

1 virtual double x_lo() const = 0;

2 virtual double x_hi() const = 0;

Description:

The curve is defined over a finite domain $[a,b]$ . x_lo() and x_hi() are implemented by concreted derived classes to return the boundary $a$ and $b$ respectively.

3.1.2 Fitted Curve

1 virtual const Function1D& curve() const = 0;

Description:

This function is implemented by concrete derived classes to retrieve the underlying curve object $f(x;\theta)$ after the fitting.

3.1.3 Roughness Measure

1 enum Regularization

2 {

3 INT_SQUARE_D0,

4 INT_SQUARE_D1,

5 INT_SQUARE_D2,

6 INT_SQUARE_D3,

7 SUM_SQUARE_DIFF_1,

8 SUM_SQUARE_DIFF_2,

9 L2_PARAM

10 };

11 void set_reg_form(Regularization reg);

Description:

The above enumeration lists available regularization choices of curve roughness measure $\mathbf{R}(\theta)$ . Users can call the property setting function set_reg_form() to choose the roughness measure to be used in curve fitting.

1.

INT_SQUARE_D0: $\mathbf{R}(\theta)=\int_{a}^{b}f(x;\theta)^{2}dx$
2.

INT_SQUARE_D1: $\mathbf{R}(\theta)=\int_{a}^{b}f^{\prime}(x;\theta)^{2}dx$
3.

INT_SQUARE_D2: $\mathbf{R}(\theta)=\int_{a}^{b}f^{\prime\prime}(x;\theta)^{2}dx$
4.

INT_SQUARE_D3: $\mathbf{R}(\theta)=\int_{a}^{b}f^{\prime\prime\prime}(x;\theta)^{2}dx$
5.

SUM_SQUARE_DIFF_1: $\mathbf{R}(\theta)=\sum_{k}\left(\frac{f^{\prime}(x_{k};\theta)-f^{\prime}(x_{% k+1};\theta)}{x_{k+1}-x_{k}}\right)^{2}$ over predefined curve intervals $(x_{k},x_{k+1})$ .
6.

SUM_SQUARE_DIFF_2: $\mathbf{R}(\theta)=\sum_{k}\left(\frac{f^{\prime\prime}(x_{k};\theta)-f^{% \prime\prime}(x_{k+1};\theta)}{x_{k+1}-x_{k}}\right)^{2}$ over predefined curve intervals $(x_{k},x_{k+1})$ .
7.

L2_PARAM: $\mathbf{R}(\theta)=|\theta|_{2}$ .

3.1.4 Roughness Leverage

It is possible to partition $[a,b]$ into sub intervals $(x_{0},x_{1},\cdots,x_{k+1})$ and apply different roughness weights in different sub interval.

\mathbf{R}(\theta)=\sum_{i=0}^{k}w_{i}\int_{x_{i}}^{x_{i+1}}\left(f^{(d)}(x;% \theta)\right)^{2}dx

The sub interval weights $\{w_{i}\}$ play the role as micro leverage factors and can be set through a piecewise constant function $w(x)$ .

1 void set_leverage(const FunctPiecewiseConst& leverage);

2 bool has_leverage() const;

3 void remove_leverage();

Description:

1.

Set the leverage factors $\{w_{i}\}$ through $w(x)$ . Parameter leverage is a reference to $w(x)$ .
2.

Return if the current curve fitting object has micro leverage factors.
3.

Clear all micro leverage factors stored in the current fitting object.

3.1.5 Shape Constraints

Users can guide the curve fitting by restricting $l_{k}\leq f^{(d)}(x_{k};\theta)\leq u_{k}$ , where $f^{(d)},d=0,1,2$ is the $d$ -th derivative of the curve function $f(x;\theta)$ wrt $x$ . Effectively, it applies a numerical bound to curve value $f(x;\theta)$ , slope $f^{\prime}(x;\theta)$ and curvature $f^{\prime\prime}(x;\theta)$ at a given $x_{k}$ . The constraint specification is defined as a nested struct SmoothCurveFitting::Constraint.

1 struct Constraint

2 {

3 enum Type

4 {

5 VALUE,

6 SLOPE,

7 CURVATURE,

8 };

10 Bound value;

11 double x;

12 Type type;

14 Constraint();

15 Constraint(const Bound& b, double xx, Type t);

16 bool operator< (const Constraint& other) const;

17 };

Description:

1.

Enumeration Type indicates the constraint type $f(x;\theta)$ , $f^{\prime}(x;\theta)$ or $f^{\prime\prime}(x;\theta)$ .
2.

Data member Bound value is the numeric bound of the constraint.
3.

Data member double x is the location of the constrained point $x$ .
4.

Data member Type type specifies the constraint type.
5.

The first constructor creates an empty constraint object.
6.

The second constructor initializes the constraint object with numerical bound parmater b, location parameter xx and type paramter t.
7.

The comparison operator for internal use.

Users can add or remove constraints by calling the following member functions of class SmoothCurveFitting.

1 void add_bc_lo(const Bound& bv, Constraint::Type bc);

2 void add_bc_hi(const Bound& bv, Constraint::Type bc);

3 void add_interior_constr(const Constraint& c);

4 void clear_constr();

Description:

1.

Add a constraint at lower boundary $a$ . Parameter bv is the numberic bound of this constraint. Parameter bc indicates the constraint type.
2.

Add a constraint at upper boundary $b$ . Parameter bv is the numberic bound of this constraint. Parameter bc indicates the constraint type.
3.

Add a constraint at an interior point. Parameter c is the constraint object.
4.

Remove all constraints stored in the curve fitting object.

3.2 Data Fitting Operators

Given $m$ -element vectors $\mathbf{x}$ and $\mathbf{y}$ , where $(x_{i},y_{i}),i=0,\cdots,m-1$ are the observed values of $x$ and function $f(x)$ , we construct a $m\times n$ basis matrix $\mathbf{\Phi}(\mathbf{x})$ where its $(i,j)$ -th element is $\Phi_{i,j}(\mathbf{x})=\phi_{j}(x_{i})$ . Effectively $\mathbf{\Phi}(\mathbf{x})\theta-\mathbf{y}$ is the vector of fitting errors (residuals). Based on the user’s choice of curve roughness measure, we construct a $n\times n$ regularization (roughness) matrix $\mathbf{\Omega}$ .

Two data fitting methodologies are provided to recover parameter $\theta$ from the curve data. They are implemented by two public operators in the LinearModelFitting base class. Derived class SmoothSplineCurve and SmoothChebyshevCurve share the same fitting operators.

3.2.1 Regularized Least Squre

1 void operator()(const Matrix& Xy, double lambda = 0);

Description:

This functions recovers $\theta$ by solving the following regularized linear least square problem.

\arg\min_{\theta}\left(\mathbf{y}-\mathbf{\Phi}(\mathbf{x})\theta\right)^{T}% \left(\mathbf{y}-\mathbf{\Phi}(\mathbf{x})\theta\right)+\lambda\theta^{T}% \mathbf{\Omega}\theta

(9)

, possibly subject to the shape constraints,

l_{k}\leq f^{(d)}(x_{k};\theta)\leq u_{k},\forall k

(10)

Parameters:

1.

Parameter Xy is a $m\times 2$ matrix. Its first and second columns hold the observation vector $\mathbf{x}$ and $\mathbf{y}$ .
2.

Parameter lambda is the regularization penalty $\lambda$ in equation (9). If $\lambda\leq 0$ is supplied, we compute an optimal penalty factor $\lambda^{*}$ using the Generalized Cross Validation (GCV) method.

3.2.2 Error as Constraints

1 void operator()(

2 const Matrix& Xy,

3 const Matrix& eps,

4 size_t max_iter,

5 double elastic_penalty = 0);

Description:

This functions recovers $\theta$ by minimizing the curve roughness subject to fitting errors constraints,

\displaystyle\arg\min_{\theta}\theta^{T}\mathbf{\Omega}\theta\mbox{ s.t.}-% \epsilon<\mathbf{\Phi}(\mathbf{x})\theta-\mathbf{y}<\epsilon

(11)

and possibly additional curve shape constraints in equation (10).

Parameters:

1.

Parameter Xy is a $m\times 2$ matrix. Its first and second columns hold the observation vector $\mathbf{x}$ and $\mathbf{y}$ .
2.

Parameter eps is a $m\times 1$ vector of tolerance $\epsilon$ as fitting error constraints in equation (11).
3.

Parameter max_iter is the max number of iterations $M$ to solve (11).
4.

If parameter elastic_penalty $>0$ is provided, the fitting error constraints are treated as elastic constraints.

3.3 Model Fitting Operators

In the context of model fitting, we design a mathematical model $g(x;\theta)=\phi(x)^{T}\theta$ to be a smooth curve. The dependent variable $y$ is connected to independent variable $x$ by a known nonlinear mapping function through the model, i.e. $y=f(g(x;\theta),x)$ .

Given model parameter $\theta$ , let $\mathbf{r}(\theta)$ be the vector of residual functions where the $i$ -th element $r_{i}(\theta)=f(g(x_{i};\theta),x_{i})-y_{i}$ is the residual error for observation $(x_{i},y_{i})$ . The LinearModelFitting base class provides two model fitting operators. Derived class SmoothSplineCurve and SmoothChebyshevCurve share the same fitting operators.

3.3.1 Regularized Least Squre

1 void operator()(

2 const VtrValueFunctionND& residuals,

3 size_t max_iter,

4 bool exception_on_max_iter,

5 double residual_tolerance,

6 double converge_tolerance,

7 double max_step_move,

8 double lambda = 0);

Description:

This fitting operator calibrates model parameter $\theta$ by solving the following regularized nonlinear least square problem.

\arg\min_{\theta}|\mathbf{r}(\theta)|_{2}^{2}+\lambda\theta^{T}\mathbf{\Omega}\theta

(12)

and possibly additional curve shape constraints in equation (10).

Parameters:

1.

residuals is the vector of residual functions $\mathbf{r}(\theta)$ .
2.

max_iter is the max number of iterations allowed to solve problem 12).
3.

exception_on_max_iter controls whether an exception should be thrown if the max number of iterations are reached.
4.

residual_tolerance is the residual tolerance $\epsilon_{r}$ for termination control. Details in white paper.
5.

converge_tolerance is the convergency tolerance $\epsilon_{c}$ for termination control. Details in white paper.
6.

max_step_move is the max step change $\Delta x$ at each iteration. Details in white paper.
7.

lambda is the penalty factor $\lambda$ in equation (12). If $\lambda\leq 0$ is supplied, GCV method is used.

3.3.2 Error as Constraints

1 void operator()(

2 const VtrValueFunctionND& residuals,

3 const Matrix& eps,

4 size_t max_iter,

5 bool exception_on_max_iter,

6 double converge_tolerance,

7 double max_step_move);

Description:

This fitting operator calibrates model parameter $\theta$ by minimizing the curve roughness subject to residual error constraints,

\arg\min_{\theta}\theta^{T}\mathbf{\Omega}\theta\mbox{ s.t. }-\epsilon<\mathbf% {r}(\theta)<\epsilon

(13)

and possibly additional curve shape constraints in equation (10).

Parameters:

1.

residuals is the vector of residual functions $\mathbf{r}(\theta)$ .
2.

eps is the vector of residual tolerance in equation (13).
3.

max_iter is the max number of iterations allowed to solve problem 13).
4.

exception_on_max_iter controls whether an exception should be thrown if the max number of iterations are reached.
5.

converge_tolerance is the convergency tolerance $\epsilon_{c}$ for termination control. Details in white paper.
6.

max_step_move is the max step change $\Delta x$ at each iteration. Details in white paper.

3.4 Model Parameter Access

The following functions are public interface functions declared in base class LinearModelFitting and implemented in concrete curve subclasses.

1 virtual size_t n_parameters() const = 0;

2 virtual void set_parameters(const double* _beta) = 0;

3 virtual void set_parameters(const Matrix& beta) = 0;

Description:

1.

Return the size of model parameter $\theta$ .
2.

Explicitly set the values of model parameter $\theta$ passed from a point parameter _beta.
3.

Explicitly set the values of model parameter $\theta$ passed from a matrix object reference parameter beta.

3.5 Smooth Spline Curve

Two constructors are provided to create a smooth curve using B-spline as the basis functions $\phi(x)$ .

1 SmoothSplineCurve(

2 const ftl::Buffer<double>& knots,

3 size_t degree,

4 size_t n_constraints,

5 Regularization reg = INT_SQUARE_D2);

7 SmoothSplineCurve(

8 const double* knots,

9 size_t n_knots,

10 size_t degree,

11 size_t n_constraints,

12 Regularization reg = INT_SQUARE_D2);

Parameters:

1.

The buffer parameter knots in the first constructor specifies the knot points $(x_{0},\cdots,x_{m-1})$ of B-spline basis. $x_{0}$ is the lower boundary point $a$ . $x_{m-1}$ is the upper boundary point $b$ .
2.

The pointer parameter knots and size parameter n_knots in the second constructor together specify the knot points $(x_{0},\cdots,x_{m-1})$ of B-spline basis. $x_{0}$ is the lower boundary point $a$ . $x_{m-1}$ is the upper boundary point $b$ .
3.

degree specifies the polynomial degree of B-spline basis.
4.

n_constraints specifies the max number of curve shape constraints to be added for fitting.
5.

reg specifies the choice of curve roughness measure.

User can directly access the fitted B-spline curve object from the member function of SmoothSplineCurve class.

1 const BSplineCurve& spline() const;

3.6 Smooth Chebyshev Curve

The following constructor is provided to create a smooth curve using Chebyshev polynomial of the first kind as the basis functions $\phi(x)$ .

1 SmoothChebyshevCurve(

2 size_t n_degree,

3 size_t n_constraints,

4 double a,

5 double b,

6 Regularization reg = INT_SQUARE_D2);

Parameters:

1.

n_degree is the highest degree of Chebyshev polynomial basis to be used.
2.

n_constraints specifies the max number of curve shape constraints to be added for fitting.
3.

a is the lower boundary of curve domain $[a,b]$ .
4.

b is the upper boundary of curve domain $[a,b]$ .
5.

reg specifies the choice of curve roughness measure.

User can directly access the fitted Chebyshev curve object from the member function of SmoothChebyshevCurve class.

1 const ChebyshevCurve& cheby() const;

4 Surface Fitting

We represent a 2-d surface function $f(x,y;\mathbf{\Theta})$ as the tensor product of two basis vectors $\phi(x)$ and $\psi(y)$ ,

f(x,y;\mathbf{\Theta})=\phi(x)^{T}\mathbf{\Theta}\psi(y)

, where $\mathbf{\Theta}$ is the coefficient matrix of the tensor product. Both $\phi(x)$ and $\psi(y)$ are of the same type of basis functions.

Mathwrist NPL offers three concrete surfaces types based on the choices of basis functions. Class SmoothSplineSurface uses B-spline as basis. Class SmoothChebyshevSurface uses Chebyshev polynomial of the first kind as basis. A special case of B-spline basis is when $\phi(x)$ and $\psi(y)$ are piecewise linear basis. This special case is separately treated in class SmoothBilinearSurface. All these concrete surfaces are derived from SmoothSurfaceFitting class that in turn inherits from top level LinearModelFitting class.

1 class LinearModelFitting

2 {

3 };

4 class SmoothSurfaceFitting : public LinearModelFitting

5 {

6 };

7 class SmoothSplineSurface : public SmoothSurfaceFitting

8 {

9 };

10 class SmoothChebyshevSurface : public SmoothSurfaceFitting

11 {

12 };

13 class SmoothBilinearSurface : public SmoothSurfaceFitting

14 {

15 };

The relevant C++ header files are:

1.

linear_model_fitting.h
2.

smooth_surface_fitting.h
3.

smooth_spline_surface.h
4.

smooth_chebyshev_surface.h
5.

smooth_bilinear_surface.h

4.1 Smooth Surface Fitting

This section documents the public interface of the SmoothSurfaceFitting base class, shared by SmoothSplineSurface, SmoothChebyshevSurface and SmoothBilinearSurface derived classes.

4.1.1 Surface Domain

1 virtual double x_lo() const = 0;

2 virtual double x_hi() const = 0;

3 virtual double y_lo() const = 0;

4 virtual double y_hi() const = 0;

Description:

The surface is defined over a rectangular domain $[a,b]\times[c,d]$ . x_lo() and x_hi() return the boundary $a$ and $b$ respectively. y_lo() and y_hi() return the boundary $c$ and $d$ respectively. These virtual functions are implemented by concrete sub-classes.

4.1.2 Fitted Surface

1 virtual const FunctionND& surf() const = 0;

Description:

This function is implemented by concrete sub-classes to retrieve the underlying surface object $f(x,y;\mathbf{\Theta})$ after fitting.

4.1.3 Roughness Measure

1 enum Regularization

2 {

3 DIRICHLET_ENERGY,

4 THIN_PLATE_ENERGY,

5 FROBENIUS_NORM

6 };

7 void set_reg_form(Regularization reg);

Description:

The above enumeration lists available regularization choices of surface roughness measure $\mathbf{R}(\mathbf{\Theta})$ . Users can call the property setting function set_reg_form() to choose the roughness measure to be used in surface fitting.

1.

DIRICHLET_ENERGY:

$\mathbf{R}(\mathbf{\Theta};a,b,c,d)=\int_{a}^{b}\int_{c}^{d}\|\nabla f(x,y;% \mathbf{\Theta})\|^{2}dydx$ .
2.

THIN_PLATE_ENERGY:

$\mathbf{R}(\mathbf{\Theta};a,b,c,d)=\int_{a}^{b}\int_{c}^{d}\left(\nabla_{xx}f% (x,y;\mathbf{\Theta})^{2}+2\nabla_{xy}f(x,y;\mathbf{\Theta})^{2}+\nabla_{yy}f(% x,y;\mathbf{\Theta})^{2}\right)dydx$
3.

FROBENIUS_NORM: $\mathbf{R}(\mathbf{\Theta})$ is the Frobenius norm of $\mathbf{\Theta}$ .

4.1.4 Roughness Leverage

It is possible to partition the whole surface domain $[a,b]\times[c,d]$ into sub areas $[x_{i},x_{i+1}]\times[y_{j},y_{j+1}]$ and apply different roughness weights $w_{i,j}$ in different sub areas.

\mathbf{R}(\mathbf{\Theta};a,b,c,d)=\sum_{i=0}^{k}\sum_{j=0}^{l}w_{i,j}\mathbf% {R}(\mathbf{\Theta};x_{i},x_{i+1},y_{j},y_{j+1})

The sub area weights $\{w_{i,j}\}$ play the role as micro leverage factors and can be set through a 2-d piecewise constant function $w(x,y)$ .

1 void set_leverage(const FunctPiecewise2DConst& leverage);

2 bool has_leverage() const;

3 void remove_leverage();

Description:

1.

Set the leverage factors $\{w_{i,j}\}$ through $w(x,y)$ . Parameter leverage is a reference to $w(x,y)$ .
2.

Return if the current surface fitting object has micro leverage factors.
3.

Clear all micro leverage factors stored in the current fitting object.

4.1.5 Shape Constraints

Users can guide the surface fitting by imposing constraints

l_{k}\leq\delta^{(r,s)}(f(x_{k},y_{k};\mathbf{\Theta}))\leq u_{k}

, where $\delta^{(r,s)}(f(x,y;\mathbf{\Theta}))$ is the operation to take $r$ -th order partial derivative wrt $x$ and then $s$ -th order partial derivative wrt $y$ . The supported values of $(r,s)$ and effective constraint types are the following:

1.

$(r=0,s=0)$ : function value $f(x_{k},y_{k},\mathbf{\Theta})$ is in range $(l_{k},u_{k})$ ;
2.

$(r=1,s=0)$ : partial delta $\frac{\partial}{\partial x}f(x_{k},y_{k},\mathbf{\Theta})$ is in range $(l_{k},u_{k})$ ;
3.

$(r=0,s=1)$ : partial delta $\frac{\partial}{\partial y}f(x_{k},y_{k},\mathbf{\Theta})$ is in range $(l_{k},u_{k})$ ;
4.

$(r=2,s=0)$ : gamma $\frac{\partial^{2}}{\partial^{2}x}f(x_{k},y_{k},\mathbf{\Theta})$ is in range $(l_{k},u_{k})$ ;
5.

$(r=0,s=2)$ : gamma $\frac{\partial^{2}}{\partial^{2}y}f(x_{k},y_{k},\mathbf{\Theta})$ is in range $(l_{k},u_{k})$ ;
6.

$(r=1,s=1)$ : cross gamma $\frac{\partial^{2}}{\partial x\partial y}f(x_{k},y_{k},\mathbf{\Theta})$ is in range $(l_{k},u_{k})$ ;

The constraint specification is defined as a nested struct SmoothSurfaceFitting::Constraint.

1 struct Constraint

2 {

3 enum Type

4 {

5 VALUE,

6 DELTA_X,

7 DELTA_Y,

8 GAMMA_X,

9 GAMMA_Y,

10 GAMMA_CROSS

11 };

13 Bound value;

14 double x;

15 double y;

16 Type type;

18 Constraint();

19 Constraint(const Bound& b, double xv, double yv, Type t);

20 bool operator< (const Constraint& other) const;

21 };

Description:

1.

Enumeration Type indicates the constraint type, effectively the $\delta^{(r,s)}$ operation.
2.

Data member Bound value is the numerical bound of the constraint.
3.

Data member double x is the $x$ coordinate of the constrained point.
4.

Data member double y is the $y$ coordinate of the constrained point.
5.

Data member Type type specifies the constraint type.
6.

The first constructor creates an empty constraint object.
7.

The second constructor initializes the constraint object with numerical bound parameter b, $x$ coordinate parameter xv, $y$ coordinate parameter yv and type parameter t.
8.

The comparison operator for internal use.

Users can add or remove constraints by the following functions. All constraints will be stored in SmoothSurfaceFitting object and applied to surface fitting.

1 void add_constr(const Constraint& c);

2 void clear_constr();

4.2 Data Fitting Operators

Given $N$ number of data points $(x_{k},y_{k},z_{k})$ , $k=0,\cdots,N-1$ , two data fitting methodologies are provided to recover parameter $\mathbf{\Theta}$ . They are implemented by two public operators in the LinearModelFitting base class.

4.2.1 Regularized Least Squre

1 void operator()(const Matrix& Xy, double lambda = 0);

Description:

This functions recovers $\mathbf{\Theta}$ by solving the following regularized linear least square problem.

\arg\min_{\mathbf{\Theta}}\sum_{k=0}^{N-1}\left(z_{k}-f(\mathbf{\Theta};x_{k},% y_{k})\right)^{2}+\lambda\mathbf{R}(\mathbf{\Theta}))

(14)

, possibly subject to the shape constraints,

l_{k}\leq\delta^{(r_{k},s_{k})}(f(x_{k},y_{k};\mathbf{\Theta}))\leq u_{k},\forall k

(15)

Parameters:

1.

Parameter Xy is a $N\times 3$ matrix. The $k$ -th row of the matrix holds the observation data point $(x_{k},y_{k},z_{k})$ .
2.

Parameter lambda is the regularization penalty $\lambda$ in equation (14). If $\lambda\leq 0$ is supplied, we compute an optimal penalty factor $\lambda^{*}$ using the Generalized Cross Validation (GCV) method.

4.2.2 Error as Constraints

1 void operator()(

2 const Matrix& Xy,

3 const Matrix& eps,

4 size_t max_iter,

5 double elastic_penalty = 0);

Description:

This functions recovers $\mathbf{\Theta}$ by minimizing the surface roughness subject to fitting errors constraints,

		$\displaystyle\arg\min_{\mathbf{\Theta}}\mathbf{R}(\mathbf{\Theta})$	$\displaystyle s.t.$		(16)
		$\displaystyle-\epsilon<f(\mathbf{\Theta};x_{k},y_{k})-z_{k}<\epsilon,\forall k$			(17)

and possibly additional surface shape constraints in equation (15).

Parameters:

1.

Parameter Xy is a $N\times 3$ matrix. The $k$ -th row of the matrix holds the observation data point $(x_{k},y_{k},z_{k})$ .
2.

Parameter eps is a $N\times 1$ vector of tolerance $\epsilon$ as fitting error constraints in equation (LABEL:eqn:surface-optimization).
3.

Parameter max_iter is the max number of iterations $M$ to solve (LABEL:eqn:surface-optimization).
4.

If parameter elastic_penalty $>0$ is provided, the fitting error constraints are treated as elastic constraints.

4.3 Model Fitting Operators

In the context of model fitting, we design a mathematical model $g(x,y;\mathbf{\Theta})$ to be a smooth surface. The dependent variable $z$ is connected to independent variables $(x,y)$ by a known nonlinear mapping function through the model, i.e. $z=f(g(x,y;\mathbf{\Theta}),x,y)$ .

Given $N$ number of data points $(x_{k},y_{k},z_{k})$ , $k=0,\cdots,N-1$ , let $\mathbf{r}(\mathbf{\Theta})$ be the vector of residual functions where the $k$ -th element $r_{k}(\mathbf{\Theta})=f(g(x_{k},y_{k};\mathbf{\Theta}),x_{k},y_{k})-z_{k}$ is the residual error for observation $(x_{k},y_{k},z_{k})$ . The LinearModelFitting base class provides two model fitting operators.

4.3.1 Regularized Least Squre

1 void operator()(

2 const VtrValueFunctionND& residuals,

3 size_t max_iter,

4 bool exception_on_max_iter,

5 double residual_tolerance,

6 double converge_tolerance,

7 double max_step_move,

8 double lambda = 0);

Description:

This fitting operator calibrates model parameter $\mathbf{\Theta}$ by solving the following regularized nonlinear least square problem.

\arg\min_{\mathbf{\Theta}}|\mathbf{r}(\mathbf{\Theta})|_{2}^{2}+\lambda\mathbf% {R}(\mathbf{\Theta})

(18)

and possibly additional surface shape constraints in equation (15).

Parameters:

1.

residuals is the vector of residual functions $\mathbf{r}(\mathbf{\Theta})$ .
2.

max_iter is the max number of iterations allowed to solve problem 18).
3.

exception_on_max_iter controls whether an exception should be thrown if the max number of iterations are reached.
4.

residual_tolerance is the residual tolerance $\epsilon_{r}$ for termination control. Details in white paper.
5.

converge_tolerance is the convergency tolerance $\epsilon_{c}$ for termination control. Details in white paper.
6.

max_step_move is the max step change $\Delta x$ at each iteration. Details in white paper.
7.

lambda is the penalty factor $\lambda$ in equation (18). If $\lambda\leq 0$ is supplied, GCV method is used.

4.3.2 Error as Constraints

1 void operator()(

2 const VtrValueFunctionND& residuals,

3 const Matrix& eps,

4 size_t max_iter,

5 bool exception_on_max_iter,

6 double converge_tolerance,

7 double max_step_move);

Description:

This fitting operator calibrates model parameter $\mathbf{\Theta}$ by minimizing the surface roughness subject to residual error constraints,

\arg\min_{\mathbf{\Theta}}\mathbf{R}(\mathbf{\Theta})\mbox{ s.t. }-\epsilon<% \mathbf{r}(\mathbf{\Theta})<\epsilon

(19)

and possibly additional surface shape constraints in equation (15).

Parameters:

1.

residuals is the vector of residual functions $\mathbf{r}(\mathbf{\Theta})$ .
2.

eps is the vector of residual tolerance in equation (19).
3.

max_iter is the max number of iterations allowed to solve problem 19).
4.

exception_on_max_iter controls whether an exception should be thrown if the max number of iterations are reached.
5.

converge_tolerance is the convergency tolerance $\epsilon_{c}$ for termination control. Details in white paper.
6.

max_step_move is the max step change $\Delta x$ and $\Delta y$ at each iteration. Details in white paper.

4.4 Model Parameter Access

The following functions are public interface functions declared in base class LinearModelFitting and implemented in concrete smooth surface sub-classes.

1 virtual size_t n_parameters() const = 0;

2 virtual void set_parameters(const double* _beta) = 0;

3 virtual void set_parameters(const Matrix& beta) = 0;

Description:

1.

Return the number of model parameter $\mathbf{\Theta}$ .
2.

Explicitly set the values of model parameter $\mathbf{\Theta}$ passed from a point parameter _beta.
3.

Explicitly set the values of model parameter $\mathbf{\Theta}$ passed from a matrix object reference parameter beta.

4.5 Smooth Spline Surface

Two constructors are provided to create a smooth surface using B-spline as the basis functions for $\phi(x)$ and $\psi(y)$ .

1 SmoothSplineSurface(

2 const ftl::Buffer<double>& knots_x,

3 const ftl::Buffer<double>& knots_y,

4 size_t degree_x,

5 size_t degree_y,

6 size_t n_constraints,

7 Regularization reg = DIRICHLET_ENERGY);

9 SmoothSplineSurface(

10 const double* knots_x,

11 const double* knots_y,

12 size_t n_knots_x,

13 size_t n_knots_y,

14 size_t degree_x,

15 size_t degree_y,

16 size_t n_constraints,

17 Regularization reg = DIRICHLET_ENERGY);

Parameters:

1.

The buffer reference parametes knots_x and knots_y in the first constructor specifies the knot points of B-spline basis $\phi(x)$ and $\psi(y)$ .
2.

In the second constructor, the pointer parameter knots_x and size parameter n_knots_x together specify the B-spline basis knot points for $\phi(x)$ . The pointer parameter knots_y and size parameter n_knots_y together specify the B-spline basis knot points for $\psi(y)$ .
3.

degree_x specifies the polynomial degree of B-spline basis $\phi(x)$ .
4.

degree_y specifies the polynomial degree of B-spline basis $\psi(y)$ .
5.

n_constraints specifies the max number of surface shape constraints to be added for fitting.
6.

reg specifies the choice of surface roughness measure.

Users can directly access the fitted B-spline surface object from the member function of SmoothSplineSurface class.

1 const BSplineSurface& spline();

4.6 Smooth Chebyshev Surface

The following constructors are provided to create a smooth surface using Chebyshev polynomial of the first kind as the basis functions for $\phi(x)$ and $\psi(y)$ .

1 SmoothChebyshevSurface(

2 double a,

3 double b,

4 double c,

5 double d,

6 size_t degree_x,

7 size_t degree_y,

8 size_t n_constraints,

9 Regularization reg = DIRICHLET_ENERGY);

Parameters:

1.

a is the lower boundary of surface $x$ coordinate domain $[a,b]$ .
2.

b is the upper boundary of surface $x$ coordinate domain $[a,b]$ .
3.

c is the lower boundary of surface $y$ coordinate domain $[c,d]$ .
4.

d is the upper boundary of surface $y$ coordinate domain $[c,d]$ .
5.

degree_x specifies the highest degree of Chebyshev polynomial basis $\phi(x)$ .
6.

degree_y specifies the highest degree of Chebyshev polynomial basis $\psi(y)$ .
7.

n_constraints specifies the max number of surface shape constraints to be added for fitting.
8.

reg specifies the choice of surface roughness measure.

Users can directly access the fitted Chebyshev surface object from the member function of SmoothChebyshevSurface class.

1 const ChebyshevSurface& cheby() const;

4.7 Smooth Bilinear Surface

Two constructors are provided to create a smooth surface for the special case when $\phi(x)$ and $\psi(y)$ are linear B-spline basis.

1 SmoothBilinearSurface(

2 const ftl::Buffer<double>& knots_x,

3 const ftl::Buffer<double>& knots_y,

4 size_t n_constraints,

5 SmoothSurfaceFitting::Regularization reg = FROBENIUS_NORM);

7 SmoothBilinearSurface(

8 const double* knots_x,

9 const double* knots_y,

10 size_t n_knots_x,

11 size_t n_knots_y,

12 size_t n_constraints,

13 SmoothSurfaceFitting::Regularization reg = FROBENIUS_NORM);

Parameters:

1.

The buffer reference parametes knots_x and knots_y in the first constructor specifies the knot points of basis $\phi(x)$ and $\psi(y)$ .
2.

In the second constructor, the pointer parameter knots_x and size parameter n_knots_x together specify the knot points for $\phi(x)$ . The pointer parameter knots_y and size parameter n_knots_y together specify the knot points for $\psi(y)$ .
3.

n_constraints specifies the max number of surface shape constraints to be added for fitting.
4.

reg specifies the choice of surface roughness measure. Thin-plate energy is not supported for bilinear surface.

Users can directly access the fitted bilinear surface object from the member function of SmoothBilinearSurface class.

1 const BilinearSurface& bilinear() const;

Data and Model Fitting Mathwrist C++ API Series

Abstract

1 Linear Least Square

1.1 Generalized Ridge Regression

1.2 Constrained Least Square

2 Nonlinear Least Square

2.1 Termination Control

2.2 Max Step Change

2.3 Solver Operator

2.4 Modified Gauss-Newton

2.5 Levenberg-Marquardt

2.6 Regularization and Constraints

3 Curve Fitting

3.1 Smooth Curve Fitting

3.1.1 Curve Domain

3.1.2 Fitted Curve

3.1.3 Roughness Measure

3.1.4 Roughness Leverage

3.1.5 Shape Constraints

3.2 Data Fitting Operators

3.2.1 Regularized Least Squre

3.2.2 Error as Constraints

3.3 Model Fitting Operators

3.3.1 Regularized Least Squre

3.3.2 Error as Constraints

3.4 Model Parameter Access

3.5 Smooth Spline Curve

3.6 Smooth Chebyshev Curve

4 Surface Fitting

4.1 Smooth Surface Fitting

4.1.1 Surface Domain

4.1.2 Fitted Surface

4.1.3 Roughness Measure

4.1.4 Roughness Leverage

4.1.5 Shape Constraints

4.2 Data Fitting Operators

4.2.1 Regularized Least Squre

4.2.2 Error as Constraints

4.3 Model Fitting Operators

4.3.1 Regularized Least Squre

4.3.2 Error as Constraints

4.4 Model Parameter Access

4.5 Smooth Spline Surface

4.6 Smooth Chebyshev Surface

4.7 Smooth Bilinear Surface

Data and Model Fitting
Mathwrist C++ API Series