Data and Model Fitting
Mathwrist White Paper Series

Based on the linear regression model assumption, we can write $𝐲=𝐗\beta +ϵ$, where $ϵ$ is the residual (noise) vector. In generalized ridge regression (2), the penalty factor $\lambda$ could be an experimental choice. However, it is often desired to automatically choose the “optimal” penalty ${\lambda }^{*}$ based on the noisy level of $ϵ$.

Generalized Cross Validation (GCV) is a very successful technique to compute ${\lambda }^{*}$ based on the assumption that $ϵ$ followings a Gaussian distribution. In GCV, the optimal penalty factor ${\lambda }^{*}$ is the solution of the following 1-dimensional minimization problem.

, where $𝐇(\lambda )$ is the unique symmetric influence matrix,

If a parameter $\lambda \le 0$ is passed to function linear_fit::grr(), we will use GCV to compute ${\lambda }^{*}$ to solve the generalized ridge regression (2).

The linear model parameter $\beta$ maybe imposed to general linear constraints and simple bounds. Mathwrist NPL provides a function linear_fit::lsq() to solve the following linearly constrained linear least square problem,

, which effectively is a convex quadratic programming (QP) problem.

The classic Gauss-Newton method approximates the true Hessian matrix in equation (11) by $𝐉{(\beta )}^T𝐉(\beta )$. It uses a line search algorithm and iteratively computes a Newton search direction $𝐩$ at each step,

This Hessian approximation has quite a few limitations. First, if the Jacobian matrix $𝐉(\beta )$ doesn’t have full rank, it produces unstable model calibration. Second, if the residual $𝐫(\beta )$ is naturally large or non-negligible at certain point of the model calibration process, ignoring the second order term $𝐐(\beta )$ in equation (11) produces incorrect search direction $𝐩$.

To address this, we implemented a modified Gauss-Newton method that follows a similar idea as [2], chapter 4. Effectively, at the $k$-th iteration in the line search, we recompute the SVD decomposition $𝐉(\beta )={\mathrm{𝐔𝐒𝐕}}^T$. It can be shown that the Newton direction $𝐩$ wrt the true Hessian solves

, where $\overline{𝐫}(\beta )={𝐔}^T𝐫(\beta )$. Let ${𝐒}_d$ be the leading submatrix of $d$ number of dominant singulars in $𝐒$. Accordingly, let ${𝐕}_d$ be the first $d$ columns of $𝐕$, in other words the principle components. We test whether $\sqrt{{|𝐫(\beta )|}_{\mathrm{\infty }}}$ is small enough relative to the smallest singulars in ${𝐒}_d$. If so, we can ignore the second order term $𝐐(\beta )$ from the true Hessian ${\nabla }^2\psi (\beta )$ and just compute a direction $𝐩={𝐕}_d{𝐩}_d$ spanned by the principle components. Let ${\overline{𝐫}}_d(\beta )$ be the first $d$ elements of $\overline{𝐫}(\beta )$. Equation (12) reduces to,

If $𝐐(\beta )$ cannot be ignored, we compute it using finite difference and solve the direction $𝐩$ in the full space of $𝐕$, $𝐩=𝐕\overline{𝐩}$. Equation (12) becomes to,

The second order term $𝐐(\beta )$ could be indefinite. We use modified Cholesky to solve equation (13). This is similar to the modified Newton method documented in our unconstrained optimization white paper.

Our modified Gauss-Newton method terminates whenever any of the following conditions are satisfied in a line search iteration,

1. 1.

   the residual satisfies , where ${ϵ}_r$ is the residual tolerance.

2. 2.

   the max number of iterations $M$ is reached.

3. 3.

   the gradient in equation (9) satisfies , where ${ϵ}_c$ is the convergence tolerance.

Levenberg-Marquardt (LM) method is a special case of the trust region algorithm that uses ${𝐉}^T(\beta )𝐉(\beta )$ to approximate the true Hessian (11). At the $k$-th trust region iteration, we solve a sub problem

${𝐩}^{*}$ is a solution of the trust region subproblem (14) if and only if $\exists \lambda \ge 0$ such that

1. 1.

   $\left({𝐉}^T(\beta )𝐉(\beta )+\lambda 𝐈\right)$ is positive semidefinite and

2. 2.

   $\left({𝐉}^T(\beta )𝐉(\beta )+\lambda 𝐈\right){𝐩}^{*}=-{𝐉}^T(\beta )𝐫(\beta )$ and

3. 3.

   $\lambda \left(\mathrm{\Delta }-\parallel {𝐩}^{*}\parallel \right)=0$

The first condition is automatically satisfied in the case of LM method. Write $𝐩(\lambda )$ as a function of $\lambda$ computed from the second condition,

If , $𝐩(\lambda =0)$ is an exact solution of trust region sub problem (14). Otherwise we can always find a $\lambda \in (0,\mathrm{\infty })$ such that $\parallel 𝐩(\lambda )\parallel ={\mathrm{\Delta }}_k$. The main cost of root finding $\lambda$ is to repeatedly solve equation (15). Let $𝐉(\beta )$ have a QR decomposition $𝐉(\beta )=𝐐\left(\begin{array}{c}\hfill 𝐑\hfill \\ \hfill 0\hfill \end{array}\right)$. Based on the idea in [1], section 10.2, we can economically obtain an upper triangular ${𝐑}_{\lambda }$ from $𝐑$ such that

Also, as suggested in [1] section 4.2, we perform root finding on a more numerically stable formulation $\frac{1}{{\mathrm{\Delta }}_k}-\frac{1}{\parallel 𝐩(\lambda )\parallel }=0$. The LM method terminates whenever any of the following conditions are satisfied after solving a trust region sub problem (14),

1. 1.

   the residual satisfies , where ${ϵ}_r$ is the residual tolerance.

2. 2.

   the max number of iterations $M$ is reached.

3. 3.

   the gradient equation (9) satisfies , where ${ϵ}_c$ is the convergence tolerance.

In practice, it is often desired to regularize the model parameters $\beta$ and perhaps impose additional constraints. Mathwrist NPL offers a general nonlinear least square fit method that recovers model parameters from the following optimization problem,

, where again $r_i(\beta )=h({𝐱}_i;\beta )-y_i$ and the second term is a regularization term. $𝛀$ is at least positive semidefinete and typically choosen to make the model “smooth”, or close to certain desired level. The penalty factor $\lambda >0$ could be a user input or automatically computed by the Generalized Cross Validation (GCV) method.

If the constraint specification in (18) is missing, our general nonlinear least square uses the line search modified Newton method to solve an unconstrained optimiation problem. Otherwise, it solves a linearly constrained optimization problem using a combination of line search modified Newton method and active-set method. The algorithm terminates whenever any of the following conditions are satisfied at a major iteration,

1. 1.

   the residual satisfies , where ${ϵ}_r$ is the residual tolerance.

2. 2.

   the max number of iterations $M$ is reached.

3. 3.

   both the stationary and convergence conditions are satisfied with ${ϵ}_c$, where ${ϵ}_c$ is the convergence tolerance.

In the first approach, we construct a basis matrix $𝚽(x)$ where the $(i,j)$-th element of the matrix is ${\mathrm{\Phi }}_{i,j}(x)={\varphi }_j(x_i)$. The sum of square of residuals (SSR) is

We recover the curve model parameter $\theta$ from data $(𝐱,𝐲)$ with desired “smoothness” property by solving a regularized least square problem,

The second term in formulation (19) is a regularization term that penalizes the roughness of curve $f(x)$. The roughness matrix $𝛀$ needs be at least positive semi-definite. The penalty factor $\lambda >0$ could be a user input or automatically computed by GCV method.

Alternatively, we can formulate the smooth curve fitting problem as minimizing the roughness of the curve and subject to the constraints that fitting errors are under certain tolerance level.

We provide two choices of basis functions $\varphi (x)$, namely B-spline and Chebyshev polynomial of the first kind. B-spline basis are piecewise polynomials of the same degree but different “locality”. Chebyshev polynomials are orthogonal polynomials of different degree but on the same global canonical domain $[-1,1]$. Internally, we map physical domain $[a,b]$ to canonical domain $[-1,1]$. These two different basis choices are represent by C++ class SmoothSplineCurve and SmoothChebyshevCurve respectively.

Because B-spline basis are overlapped piecewise polynomials, the placement of knot points has significant impact to the result of fitting. In general, it is desired to have data points uniformally scattered across all pieces. If the data implies dramatic curvature changes, adding an extra connecting piece could help smooth transition. Based on our experience, B-spline could be vulnerable to data noises. Careful attention to regularization and perhaps with additional shape constraints are often needed.

Chebyshev polynomial has many nice theoretical properties. Most relevant to data fitting, it is very close to the ideal “minimax” approximation polynomial. In other words, Chebyshev polynomials are somewhat naturally “smooth” and numerically more stable. Each Chebyshev basis function is a “global” polynomial over $[-1,1]$. We just need determine the number $n$ of basis to be used in $f(x)$, which is the highest degree of Chebyshev polynomial series. Increasing $n$ of course increases the freedom of fitting. Practically, very large degree $n$, i.e. greater than 7 or 8, is often unnecessary.

The regularization choice in formulation (19) and objective function in formulation (21) are to make the curve function $f(x;\theta )$ “smoother”. If one chooses $𝛀$ being the identity matrix, the roughness measure is to reduce the $l_2$ norm of basis coefficients ${\parallel \theta \parallel }_2^2$, which tends to produce a flat curve closer and closer to $f(x;\theta )=0$ as we increase the penalty factor $\lambda \to \mathrm{\infty }$.

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

The total roughness measure $𝐑(\theta )$ then is a weighted sum of local roughness. The sub interval weights $𝐰=\{w_i\}$ play the role as micro leverage factors. Users can set the leverage factors $𝐰$ through a piecewise constant function $w(x)$.

In both formulations (19) and (21), we may impose additional shape constraints at selected points $x_k,k=0,\mathrm{\cdots },K$,

The shape constraints (22) effectively restrict the function value $f(x_k;\theta )$, slope $f^{\prime }(x_k;\theta )$ or curvature $f^{\prime \prime }(x_k;\theta )$ to be bounded within a certain range. For example, we could impose $f^{\prime \prime }(a;\theta )=0$ and $f^{\prime \prime }(b;\theta )=0$ constraints to fit a curve that satisfies natural boundary conditions.

Here, we design a mathematical model $g(x;\theta )=\varphi {(x)}^T\theta$ that is to be approximated by 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)$.

Let $𝐫(\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 formulation (19) changes to the following regularized nonlinear least square problem,

And formulation (21) now becomes to an nonlinear programming problem,

Again, additional curve shape constraint (22) may be imposed to both formulation (23) and (25). These two model calibration methods internally use our linearly constraint (LC) optimization solver and nonlinear programming (NLP) solver respectively. Please refer to our LC and NLP white paper for further details.

The basis functions $\varphi (x)$ and $\psi (y)$ are of the same type, either B-spline or Chebyshev polynomials of the first kind. These two basis choices are represented by two different C++ surface class types, SmoothSplineSurface and SmoothChebyshevSurface respectively. As a special case of B-spline basis, when we choose $\varphi (x)$ and $\psi (y)$ to be piecewise linear functions, function $f(x,y)$ is a bi-linear surface. This special case deserves separate treatment. It is represented by a separate surface C++ class SmoothBilinearSurface for more efficient calculation.

We offer 3 different surface roughness measures, Frobenius norm, Dirichlet energy and Thin-plate energy. Frobenius norm is equivalent to the $l_2$ norm in curve fitting situation.

Let $[a,b]\times [c,d]$ be the domain of the surface. Dirichlet energy is defined as the square integral of the gradient norm,

The standard “thin-plate energy” is a rotation-invariant measure defined as the following,

Thin-plate energy roughness is not applicable to bilinear surface.

Users can partition the whole surface domain $[a,b]\times [c,d]$ into sub areas $[x_i,x_{i+1}]\times [y_j,y_{j+1}]$, $i=0,\mathrm{\cdots },k$, $j=0,\mathrm{\cdots },l$ and supply a 2-d piecewise constant function $w(x,y)$ to surface fitting.

Let $w_{i,j}$ be the weights computed from $w(x,y)$ for area $[x_i,x_{i+1}]\times [y_j,y_{j+1}]$. We then compute the total roughness measure as the weighted sum of roughness over those sub surface areas.

In both formulation (27) and (29), it is possible to impose constraints directly to certain points $(x_k,y_k)$ on the surface. The supported constraint types are:

1. 1.

   Function value $f(x_k,y_k,𝚯)$ is bounded;

2. 2.

   Partial delta $\frac{\partial }{\partial x}f(x_k,y_k,𝚯)$ or $\frac{\partial }{\partial y}f(x_k,y_k,𝚯)$ is bounded;

3. 3.

   Gamma $\frac{{\partial }^2}{{\partial }^{2}x}f(x_k,y_k,𝚯)$ or $\frac{{\partial }^2}{{\partial }^{2}y}f(x_k,y_k,𝚯)$ is bounded;

4. 4.

   Cross gamma $\frac{{\partial }^2}{\partial x\partial y}f(x_k,y_k,𝚯)$ is bounded;

Let ${\delta }^{(r,s)}(f(x,y;𝚯))$ be the differential operator relevant to the supported constraint types, i.e. ${\delta }^{(0,0)}$ for function value, ${\delta }^{(1,0)}$ for partial delta $\frac{\partial }{\partial x}$. The optimization problem (27) and (29) may be subject to surface shape constraints,

Similar to nonlinear curve model fitting, here we have a known 2-d nonlinear mapping function $z=f(g(x,y;𝚯),x,y)$ that connects independent variables $(x,y)$ to function value $z$ through a model $g(x,y;𝚯)$. It is the model function $g(x,y;𝚯)$ that is approximated by a smooth surface, i.e. $g(x,y;𝚯)=\varphi {(x)}^T𝚯\psi (y)$.

Let $𝐫(𝚯)$ be the vector of residual functions where the $i$-th element $r_i(𝚯)=f(g(x_i,y_i;𝚯),x_i,y_i)-z_i$ is the residual error for observation $(x_i,y_i,z_i)$. Users can choose to calibrate surface model parameter $𝚯$ either by a regularized nonlinear least square fit,

, or by solving an nonlinear programming problem,

Again, in both formulation (31) and (33), additional shape constraints (30) can be imposed.