.. index::

   single: fitting

   single: least squares fit

   single: regression, least squares

   single: weighted linear fits

   single: unweighted linear fits

****************************

Linear Least-Squares Fitting

****************************

This chapter describes routines for performing least squares fits to

experimental data using linear combinations of functions.  The data

may be weighted or unweighted, i.e. with known or unknown errors.  For

weighted data the functions compute the best fit parameters and their

associated covariance matrix.  For unweighted data the covariance

matrix is estimated from the scatter of the points, giving a

variance-covariance matrix.

The functions are divided into separate versions for simple one- or

two-parameter regression and multiple-parameter fits.

.. _sec_lls-overview:

Overview

========

Least-squares fits are found by minimizing :math:`\chi^2`

(chi-squared), the weighted sum of squared residuals over :math:`n`

experimental datapoints :math:`(x_i, y_i)` for the model :math:`Y(c,x)`,

.. math:: \chi^2 = \sum_i w_i (y_i - Y(c, x_i))^2

The :math:`p` parameters of the model are :math:`c = \{c_0, c_1, \dots\}`.

The weight factors :math:`w_i` are given by :math:`w_i = 1/\sigma_i^2`

where :math:`\sigma_i` is the experimental error on the data-point

:math:`y_i`.  The errors are assumed to be

Gaussian and uncorrelated. 

For unweighted data the chi-squared sum is computed without any weight factors. 

.. index::

   single: covariance matrix, linear fits

The fitting routines return the best-fit parameters :math:`c` and their

:math:`p \times p` covariance matrix.  The covariance matrix measures the

statistical errors on the best-fit parameters resulting from the 

errors on the data, :math:`\sigma_i`, and is defined

as

.. only:: not texinfo

   .. math:: C_{ab} = \langle \delta c_a \delta c_b \rangle

.. only:: texinfo

   ::

      C_{ab} = <\delta c_a \delta c_b>

where :math:`\langle \, \rangle`

denotes an average over the Gaussian error distributions of the underlying datapoints.

The covariance matrix is calculated by error propagation from the data

errors :math:`\sigma_i`.  The change in a fitted parameter :math:`\delta c_a`

caused by a small change in the data :math:`\delta y_i` is given

by

.. only:: not texinfo

   .. math:: \delta c_a = \sum_i {\partial c_a \over \partial y_i} \delta y_i

.. only:: texinfo

   ::

      \delta c_a = \sum_i (dc_a/dy_i) \delta y_i

allowing the covariance matrix to be written in terms of the errors on the data,

.. only:: not texinfo

   .. math::

      C_{ab} =  \sum_{i,j} {\partial c_a \over \partial y_i}

                           {\partial c_b \over \partial y_j} 

                           \langle \delta y_i \delta y_j \rangle

.. only:: texinfo

   ::

      C_{ab} = \sum_{i,j} (dc_a/dy_i) (dc_b/dy_j) <\delta y_i \delta y_j>

For uncorrelated data the fluctuations of the underlying datapoints satisfy

.. only:: not texinfo

   .. math:: \langle \delta y_i \delta y_j \rangle = \sigma_i^2 \delta_{ij}

.. only:: texinfo

   ::

      <\delta y_i \delta y_j> = \sigma_i^2 \delta_{ij}

giving a corresponding parameter covariance matrix of

.. only:: not texinfo

   .. math:: C_{ab} = \sum_{i} {1 \over w_i} {\partial c_a \over \partial y_i} {\partial c_b \over \partial y_i} 

.. only:: texinfo

   ::

      C_{ab} = \sum_i (1/w_i) (dc_a/dy_i) (dc_b/dy_i) 

.. index::

   single: variance-covariance matrix, linear fits

When computing the covariance matrix for unweighted data, i.e. data with unknown errors, 

the weight factors :math:`w_i` in this sum are replaced by the single estimate

:math:`w = 1/\sigma^2`, where :math:`\sigma^2` is the computed variance of the

residuals about the best-fit model, :math:`\sigma^2 = \sum (y_i - Y(c,x_i))^2 / (n-p)`.  

This is referred to as the *variance-covariance matrix*.

The standard deviations of the best-fit parameters are given by the

square root of the corresponding diagonal elements of

the covariance matrix, :math:`\sigma_{c_a} = \sqrt{C_{aa}}`.

The correlation coefficient of the fit parameters :math:`c_a` and :math:`c_b`

is given by :math:`\rho_{ab} = C_{ab} / \sqrt{C_{aa} C_{bb}}`.

.. index:: linear regression

Linear regression

=================

The functions in this section are used to fit simple one or two

parameter linear regression models. The functions are declared in

the header file :file:`gsl_fit.h`.

Linear regression with a constant term

--------------------------------------

The functions described in this section can be used to perform

least-squares fits to a straight line model, :math:`Y(c,x) = c_0 + c_1 x`.

.. index::

   single: covariance matrix, from linear regression

.. function:: int gsl_fit_linear (const double * x, const size_t xstride, const double * y, const size_t ystride, size_t n, double * c0, double * c1, double * cov00, double * cov01, double * cov11, double * sumsq)

   This function computes the best-fit linear regression coefficients

   (:data:`c0`, :data:`c1`) of the model :math:`Y = c_0 + c_1 X` for the dataset

   (:data:`x`, :data:`y`), two vectors of length :data:`n` with strides

   :data:`xstride` and :data:`ystride`.  The errors on :data:`y` are assumed unknown so 

   the variance-covariance matrix for the

   parameters (:data:`c0`, :data:`c1`) is estimated from the scatter of the

   points around the best-fit line and returned via the parameters

   (:data:`cov00`, :data:`cov01`, :data:`cov11`).   

   The sum of squares of the residuals from the best-fit line is returned

   in :data:`sumsq`.  Note: the correlation coefficient of the data can be computed using

   :func:`gsl_stats_correlation`, it does not depend on the fit.

.. function:: int gsl_fit_wlinear (const double * x, const size_t xstride, const double * w, const size_t wstride, const double * y, const size_t ystride, size_t n, double * c0, double * c1, double * cov00, double * cov01, double * cov11, double * chisq)

   This function computes the best-fit linear regression coefficients

   (:data:`c0`, :data:`c1`) of the model :math:`Y = c_0 + c_1 X` for the weighted

   dataset (:data:`x`, :data:`y`), two vectors of length :data:`n` with strides

   :data:`xstride` and :data:`ystride`.  The vector :data:`w`, of length :data:`n`

   and stride :data:`wstride`, specifies the weight of each datapoint. The

   weight is the reciprocal of the variance for each datapoint in :data:`y`.

   The covariance matrix for the parameters (:data:`c0`, :data:`c1`) is

   computed using the weights and returned via the parameters

   (:data:`cov00`, :data:`cov01`, :data:`cov11`).  The weighted sum of squares

   of the residuals from the best-fit line, :math:`\chi^2`, is returned in

   :data:`chisq`.

.. function:: int gsl_fit_linear_est (double x, double c0, double c1, double cov00, double cov01, double cov11, double * y, double * y_err)

   This function uses the best-fit linear regression coefficients

   :data:`c0`, :data:`c1` and their covariance

   :data:`cov00`, :data:`cov01`, :data:`cov11` to compute the fitted function

   :data:`y` and its standard deviation :data:`y_err` for the model :math:`Y = c_0 + c_1 X`

   at the point :data:`x`.

Linear regression without a constant term

-----------------------------------------

The functions described in this section can be used to perform

least-squares fits to a straight line model without a constant term,

:math:`Y = c_1 X`.

.. function:: int gsl_fit_mul (const double * x, const size_t xstride, const double * y, const size_t ystride, size_t n, double * c1, double * cov11, double * sumsq)

   This function computes the best-fit linear regression coefficient

   :data:`c1` of the model :math:`Y = c_1 X` for the datasets (:data:`x`,

   :data:`y`), two vectors of length :data:`n` with strides :data:`xstride` and

   :data:`ystride`.  The errors on :data:`y` are assumed unknown so the 

   variance of the parameter :data:`c1` is estimated from

   the scatter of the points around the best-fit line and returned via the

   parameter :data:`cov11`.  The sum of squares of the residuals from the

   best-fit line is returned in :data:`sumsq`.

.. function:: int gsl_fit_wmul (const double * x, const size_t xstride, const double * w, const size_t wstride, const double * y, const size_t ystride, size_t n, double * c1, double * cov11, double * sumsq)

   This function computes the best-fit linear regression coefficient

   :data:`c1` of the model :math:`Y = c_1 X` for the weighted datasets

   (:data:`x`, :data:`y`), two vectors of length :data:`n` with strides

   :data:`xstride` and :data:`ystride`.  The vector :data:`w`, of length :data:`n`

   and stride :data:`wstride`, specifies the weight of each datapoint. The

   weight is the reciprocal of the variance for each datapoint in :data:`y`.

   The variance of the parameter :data:`c1` is computed using the weights

   and returned via the parameter :data:`cov11`.  The weighted sum of

   squares of the residuals from the best-fit line, :math:`\chi^2`, is

   returned in :data:`chisq`.

.. function:: int gsl_fit_mul_est (double x, double c1, double cov11, double * y, double * y_err)

   This function uses the best-fit linear regression coefficient :data:`c1`

   and its covariance :data:`cov11` to compute the fitted function

   :data:`y` and its standard deviation :data:`y_err` for the model :math:`Y = c_1 X`

   at the point :data:`x`.

.. index::

   single: multi-parameter regression

   single: fits, multi-parameter linear

Multi-parameter regression

==========================

This section describes routines which perform least squares fits

to a linear model by minimizing the cost function

.. math:: \chi^2 = \sum_i w_i (y_i - \sum_j X_{ij} c_j)^2 = || y - Xc ||_W^2

where :math:`y` is a vector of :math:`n` observations, :math:`X` is an

:math:`n`-by-:math:`p` matrix of predictor variables, :math:`c`

is a vector of the :math:`p` unknown best-fit parameters to be estimated,

and :math:`||r||_W^2 = r^T W r`.

The matrix :math:`W = \diag(w_1,w_2,...,w_n)`

defines the weights or uncertainties of the observation vector.

This formulation can be used for fits to any number of functions and/or

variables by preparing the :math:`n`-by-:math:`p` matrix :math:`X`

appropriately.  For example, to fit to a :math:`p`-th order polynomial in

:data:`x`, use the following matrix,

.. math:: X_{ij} = x_i^j

where the index :math:`i` runs over the observations and the index

:math:`j` runs from 0 to :math:`p-1`.

To fit to a set of :math:`p` sinusoidal functions with fixed frequencies

:math:`\omega_1`, :math:`\omega_2`, :math:`\ldots`, :math:`\omega_p`, use,

.. math:: X_{ij} = \sin(\omega_j x_i)

To fit to :math:`p` independent variables :math:`x_1`, :math:`x_2`, :math:`\ldots`,

:math:`x_p`, use,

.. math:: X_{ij} = x_j(i)

where :math:`x_j(i)` is the :math:`i`-th value of the predictor variable

:math:`x_j`.

The solution of the general linear least-squares system requires an

additional working space for intermediate results, such as the singular

value decomposition of the matrix :math:`X`.

These functions are declared in the header file :file:`gsl_multifit.h`.

.. type:: gsl_multifit_linear_workspace

   This workspace contains internal variables for fitting multi-parameter models.

.. function:: gsl_multifit_linear_workspace * gsl_multifit_linear_alloc (const size_t n, const size_t p)

   This function allocates a workspace for fitting a model to a maximum of :data:`n`

   observations using a maximum of :data:`p` parameters. The user may later supply

   a smaller least squares system if desired. The size of the workspace is

   :math:`O(np + p^2)`.

.. function:: void gsl_multifit_linear_free (gsl_multifit_linear_workspace * work)

   This function frees the memory associated with the workspace :data:`w`.

.. function:: int gsl_multifit_linear_svd (const gsl_matrix * X, gsl_multifit_linear_workspace * work)

   This function performs a singular value decomposition of the

   matrix :data:`X` and stores the SVD factors internally in :data:`work`.

.. function:: int gsl_multifit_linear_bsvd (const gsl_matrix * X, gsl_multifit_linear_workspace * work)

   This function performs a singular value decomposition of the

   matrix :data:`X` and stores the SVD factors internally in :data:`work`.

   The matrix :data:`X` is first balanced by applying column scaling

   factors to improve the accuracy of the singular values.

.. function:: int gsl_multifit_linear (const gsl_matrix * X, const gsl_vector * y, gsl_vector * c, gsl_matrix * cov, double * chisq, gsl_multifit_linear_workspace * work)

   This function computes the best-fit parameters :data:`c` of the model

   :math:`y = X c` for the observations :data:`y` and the matrix of

   predictor variables :data:`X`, using the preallocated workspace provided

   in :data:`work`.  The :math:`p`-by-:math:`p` variance-covariance matrix of the model parameters

   :data:`cov` is set to :math:`\sigma^2 (X^T X)^{-1}`, where :math:`\sigma` is

   the standard deviation of the fit residuals.

   The sum of squares of the residuals from the best-fit,

   :math:`\chi^2`, is returned in :data:`chisq`. If the coefficient of

   determination is desired, it can be computed from the expression

   :math:`R^2 = 1 - \chi^2 / TSS`, where the total sum of squares (TSS) of

   the observations :data:`y` may be computed from :func:`gsl_stats_tss`.

   The best-fit is found by singular value decomposition of the matrix

   :data:`X` using the modified Golub-Reinsch SVD algorithm, with column

   scaling to improve the accuracy of the singular values. Any components

   which have zero singular value (to machine precision) are discarded

   from the fit.

.. function:: int gsl_multifit_linear_tsvd (const gsl_matrix * X, const gsl_vector * y, const double tol, gsl_vector * c, gsl_matrix * cov, double * chisq, size_t * rank, gsl_multifit_linear_workspace * work)

   This function computes the best-fit parameters :data:`c` of the model

   :math:`y = X c` for the observations :data:`y` and the matrix of

   predictor variables :data:`X`, using a truncated SVD expansion.

   Singular values which satisfy :math:`s_i \le tol \times s_0`

   are discarded from the fit, where :math:`s_0` is the largest singular value.

   The :math:`p`-by-:math:`p` variance-covariance matrix of the model parameters

   :data:`cov` is set to :math:`\sigma^2 (X^T X)^{-1}`, where :math:`\sigma` is

   the standard deviation of the fit residuals.

   The sum of squares of the residuals from the best-fit,

   :math:`\chi^2`, is returned in :data:`chisq`. The effective rank

   (number of singular values used in solution) is returned in :data:`rank`.

   If the coefficient of

   determination is desired, it can be computed from the expression

   :math:`R^2 = 1 - \chi^2 / TSS`, where the total sum of squares (TSS) of

   the observations :data:`y` may be computed from :func:`gsl_stats_tss`.

.. function:: int gsl_multifit_wlinear (const gsl_matrix * X, const gsl_vector * w, const gsl_vector * y, gsl_vector * c, gsl_matrix * cov, double * chisq, gsl_multifit_linear_workspace * work)

   This function computes the best-fit parameters :data:`c` of the weighted

   model :math:`y = X c` for the observations :data:`y` with weights :data:`w`

   and the matrix of predictor variables :data:`X`, using the preallocated

   workspace provided in :data:`work`.  The :math:`p`-by-:math:`p` covariance matrix of the model

   parameters :data:`cov` is computed as :math:`(X^T W X)^{-1}`. The weighted

   sum of squares of the residuals from the best-fit, :math:`\chi^2`, is

   returned in :data:`chisq`. If the coefficient of determination is

   desired, it can be computed from the expression :math:`R^2 = 1 - \chi^2 / WTSS`,

   where the weighted total sum of squares (WTSS) of the

   observations :data:`y` may be computed from :func:`gsl_stats_wtss`.

.. function:: int gsl_multifit_wlinear_tsvd (const gsl_matrix * X, const gsl_vector * w, const gsl_vector * y, const double tol, gsl_vector * c, gsl_matrix * cov, double * chisq, size_t * rank, gsl_multifit_linear_workspace * work)

   This function computes the best-fit parameters :data:`c` of the weighted

   model :math:`y = X c` for the observations :data:`y` with weights :data:`w`

   and the matrix of predictor variables :data:`X`, using a truncated SVD expansion.

   Singular values which satisfy :math:`s_i \le tol \times s_0`

   are discarded from the fit, where :math:`s_0` is the largest singular value.

   The :math:`p`-by-:math:`p` covariance matrix of the model

   parameters :data:`cov` is computed as :math:`(X^T W X)^{-1}`. The weighted

   sum of squares of the residuals from the best-fit, :math:`\chi^2`, is

   returned in :data:`chisq`. The effective rank of the system (number of

   singular values used in the solution) is returned in :data:`rank`.

   If the coefficient of determination is

   desired, it can be computed from the expression :math:`R^2 = 1 - \chi^2 / WTSS`,

   where the weighted total sum of squares (WTSS) of the

   observations :data:`y` may be computed from :func:`gsl_stats_wtss`.

.. function:: int gsl_multifit_linear_est (const gsl_vector * x, const gsl_vector * c, const gsl_matrix * cov, double * y, double * y_err)

   This function uses the best-fit multilinear regression coefficients

   :data:`c` and their covariance matrix

   :data:`cov` to compute the fitted function value

   :data:`y` and its standard deviation :data:`y_err` for the model :math:`y = x.c` 

   at the point :data:`x`.

.. function:: int gsl_multifit_linear_residuals (const gsl_matrix * X, const gsl_vector * y, const gsl_vector * c, gsl_vector * r)

   This function computes the vector of residuals :math:`r = y - X c` for

   the observations :data:`y`, coefficients :data:`c` and matrix of predictor

   variables :data:`X`.

.. function:: size_t gsl_multifit_linear_rank (const double tol, const gsl_multifit_linear_workspace * work)

   This function returns the rank of the matrix :math:`X` which must first have its

   singular value decomposition computed. The rank is computed by counting the number

   of singular values :math:`\sigma_j` which satisfy :math:`\sigma_j > tol \times \sigma_0`,

   where :math:`\sigma_0` is the largest singular value.

.. index::

   single: ridge regression

   single: Tikhonov regression

   single: regression, ridge

   single: regression, Tikhonov

   single: least squares, regularized

.. _sec_regularized-regression:

Regularized regression

======================

Ordinary weighted least squares models seek a solution vector :math:`c`

which minimizes the residual

.. math:: \chi^2 = || y - Xc ||_W^2

where :math:`y` is the :math:`n`-by-:math:`1` observation vector,

:math:`X` is the :math:`n`-by-:math:`p` design matrix, :math:`c` is

the :math:`p`-by-:math:`1` solution vector,

:math:`W = \diag(w_1,...,w_n)` is the data weighting matrix,

and :math:`||r||_W^2 = r^T W r`.

In cases where the least squares matrix :math:`X` is ill-conditioned,

small perturbations (ie: noise) in the observation vector could lead to

widely different solution vectors :math:`c`.

One way of dealing with ill-conditioned matrices is to use a "truncated SVD"

in which small singular values, below some given tolerance, are discarded

from the solution. The truncated SVD method is available using the functions

:func:`gsl_multifit_linear_tsvd` and :func:`gsl_multifit_wlinear_tsvd`. Another way

to help solve ill-posed problems is to include a regularization term in the least squares

minimization

.. math:: \chi^2 = || y - Xc ||_W^2 + \lambda^2 || L c ||^2

for a suitably chosen regularization parameter :math:`\lambda` and

matrix :math:`L`. This type of regularization is known as Tikhonov, or ridge,

regression. In some applications, :math:`L` is chosen as the identity matrix, giving

preference to solution vectors :math:`c` with smaller norms.

Including this regularization term leads to the explicit "normal equations" solution

.. only:: not texinfo

   .. math:: c = \left( X^T W X + \lambda^2 L^T L \right)^{-1} X^T W y

.. only:: texinfo

   ::

      c = ( X^T W X + \lambda^2 L^T L )^-1 X^T W y

which reduces to the ordinary least squares solution when :math:`L = 0`.

In practice, it is often advantageous to transform a regularized least

squares system into the form

.. only:: not texinfo

   .. math:: \chi^2 = || \tilde{y} - \tilde{X} \tilde{c} ||^2 + \lambda^2 || \tilde{c} ||^2

.. only:: texinfo

   ::

      \chi^2 = || y~ - X~ c~ ||^2 + \lambda^2 || c~ ||^2

This is known as the Tikhonov "standard form" and has the normal equations solution

.. only:: not texinfo

   .. math:: \tilde{c} = \left( \tilde{X}^T \tilde{X} + \lambda^2 I \right)^{-1} \tilde{X}^T \tilde{y}

.. only:: texinfo

   ::

      \tilde{c} = ( \tilde{X}^T \tilde{X} + \lambda^2 I )^{-1} \tilde{X}^T \tilde{y}

For an :math:`m`-by-:math:`p` matrix :math:`L` which is full rank and has :math:`m >= p` (ie: :math:`L` is

square or has more rows than columns), we can calculate the "thin" QR decomposition of :math:`L`, and

note that :math:`||L c|| = ||R c||` since the :math:`Q` factor will not change the norm. Since

:math:`R` is :math:`p`-by-:math:`p`, we can then use the transformation

.. only:: not texinfo

   .. math::

      \tilde{X} &= W^{1 \over 2} X R^{-1} \\

      \tilde{y} &= W^{1 \over 2} y \\

      \tilde{c} &= R c

.. only:: texinfo

   ::

      X~ = sqrt(W) X R^-1

      y~ = sqrt(W) y

      c~ = R c

to achieve the standard form. For a rectangular matrix :math:`L` with :math:`m < p`,

a more sophisticated approach is needed (see Hansen 1998, chapter 2.3).

In practice, the normal equations solution above is not desirable due to

numerical instabilities, and so the system is solved using the

singular value decomposition of the matrix :math:`\tilde{X}`.

The matrix :math:`L` is often chosen as the identity matrix, or as a first

or second finite difference operator, to ensure a smoothly varying

coefficient vector :math:`c`, or as a diagonal matrix to selectively damp

each model parameter differently. If :math:`L \ne I`, the user must first

convert the least squares problem to standard form using

:func:`gsl_multifit_linear_stdform1` or :func:`gsl_multifit_linear_stdform2`,

solve the system, and then backtransform the solution vector to recover

the solution of the original problem (see

:func:`gsl_multifit_linear_genform1` and :func:`gsl_multifit_linear_genform2`).

In many regularization problems, care must be taken when choosing

the regularization parameter :math:`\lambda`. Since both the

residual norm :math:`||y - X c||` and solution norm :math:`||L c||`

are being minimized, the parameter :math:`\lambda` represents

a tradeoff between minimizing either the residuals or the

solution vector. A common tool for visualizing the comprimise between

the minimization of these two quantities is known as the L-curve.

The L-curve is a log-log plot of the residual norm :math:`||y - X c||`

on the horizontal axis and the solution norm :math:`||L c||` on the

vertical axis. This curve nearly always as an :math:`L` shaped

appearance, with a distinct corner separating the horizontal

and vertical sections of the curve. The regularization parameter

corresponding to this corner is often chosen as the optimal

value. GSL provides routines to calculate the L-curve for all

relevant regularization parameters as well as locating the corner.

Another method of choosing the regularization parameter is known

as Generalized Cross Validation (GCV). This method is based on

the idea that if an arbitrary element :math:`y_i` is left out of the

right hand side, the resulting regularized solution should predict this element

accurately. This leads to choosing the parameter :math:`\lambda`

which minimizes the GCV function

.. only:: not texinfo

   .. math:: G(\lambda) = {||y - X c_{\lambda}||^2 \over \textrm{Tr}(I_n - X X_{\lambda}^I)^2}

.. only:: texinfo

   ::

      G(\lambda) = (||y - X c_{\lambda}||^2) / Tr(I_n - X X^I)^2

where :math:`X_{\lambda}^I` is the matrix which relates the solution :math:`c_{\lambda}`

to the right hand side :math:`y`, ie: :math:`c_{\lambda} = X_{\lambda}^I y`. GSL

provides routines to compute the GCV curve and its minimum.

For most applications, the steps required to solve a regularized least

squares problem are as follows:

#. Construct the least squares system (:math:`X`, :math:`y`, :math:`W`, :math:`L`)

#. Transform the system to standard form (:math:`\tilde{X}`, :math:`\tilde{y}`). This

   step can be skipped if :math:`L = I_p` and :math:`W = I_n`.

#. Calculate the SVD of :math:`\tilde{X}`.

#. Determine an appropriate regularization parameter :math:`\lambda` (using for example

   L-curve or GCV analysis).

#. Solve the standard form system using the chosen :math:`\lambda` and the SVD of :math:`\tilde{X}`.

#. Backtransform the standard form solution :math:`\tilde{c}` to recover the

   original solution vector :math:`c`.

.. function:: int gsl_multifit_linear_stdform1 (const gsl_vector * L, const gsl_matrix * X, const gsl_vector * y, gsl_matrix * Xs, gsl_vector * ys, gsl_multifit_linear_workspace * work)

              int gsl_multifit_linear_wstdform1 (const gsl_vector * L, const gsl_matrix * X, const gsl_vector * w, const gsl_vector * y, gsl_matrix * Xs, gsl_vector * ys, gsl_multifit_linear_workspace * work)

   These functions define a regularization matrix

   :math:`L = \diag(l_0,l_1,...,l_{p-1})`.

   The diagonal matrix element :math:`l_i` is provided by the

   :math:`i`-th element of the input vector :data:`L`.

   The :math:`n`-by-:math:`p` least squares matrix :data:`X` and

   vector :data:`y` of length :math:`n` are then

   converted to standard form as described above and the parameters

   (:math:`\tilde{X}`, :math:`\tilde{y}`) are stored in :data:`Xs` and :data:`ys`

   on output.  :data:`Xs` and :data:`ys` have the same dimensions as

   :data:`X` and :data:`y`. Optional data weights may be supplied in the

   vector :data:`w` of length :math:`n`. In order to apply this transformation,

   :math:`L^{-1}` must exist and so none of the :math:`l_i`

   may be zero. After the standard form system has been solved,

   use :func:`gsl_multifit_linear_genform1` to recover the original solution vector.

   It is allowed to have :data:`X` = :data:`Xs` and :data:`y` = :data:`ys` for an in-place transform.

   In order to perform a weighted regularized fit with :math:`L = I`, the user may

   call :func:`gsl_multifit_linear_applyW` to convert to standard form.

.. function:: int gsl_multifit_linear_L_decomp (gsl_matrix * L, gsl_vector * tau)

   This function factors the :math:`m`-by-:math:`p` regularization matrix

   :data:`L` into a form needed for the later transformation to standard form. :data:`L`

   may have any number of rows :math:`m`. If :math:`m \ge p` the QR decomposition of

   :data:`L` is computed and stored in :data:`L` on output. If :math:`m < p`, the QR decomposition

   of :math:`L^T` is computed and stored in :data:`L` on output. On output,

   the Householder scalars are stored in the vector :data:`tau` of size :math:`MIN(m,p)`.

   These outputs will be used by :func:`gsl_multifit_linear_wstdform2` to complete the

   transformation to standard form.

.. function:: int gsl_multifit_linear_stdform2 (const gsl_matrix * LQR, const gsl_vector * Ltau, const gsl_matrix * X, const gsl_vector * y, gsl_matrix * Xs, gsl_vector * ys, gsl_matrix * M, gsl_multifit_linear_workspace * work)

              int gsl_multifit_linear_wstdform2 (const gsl_matrix * LQR, const gsl_vector * Ltau, const gsl_matrix * X, const gsl_vector * w, const gsl_vector * y, gsl_matrix * Xs, gsl_vector * ys, gsl_matrix * M, gsl_multifit_linear_workspace * work)

   These functions convert the least squares system (:data:`X`, :data:`y`, :data:`W`, :math:`L`) to standard

   form (:math:`\tilde{X}`, :math:`\tilde{y}`) which are stored in :data:`Xs` and :data:`ys`

   respectively. The :math:`m`-by-:math:`p` regularization matrix :data:`L` is specified by the inputs

   :data:`LQR` and :data:`Ltau`, which are outputs from :func:`gsl_multifit_linear_L_decomp`.

   The dimensions of the standard form parameters (:math:`\tilde{X}`, :math:`\tilde{y}`)

   depend on whether :math:`m` is larger or less than :math:`p`. For :math:`m \ge p`,

   :data:`Xs` is :math:`n`-by-:math:`p`, :data:`ys` is :math:`n`-by-1, and :data:`M` is

   not used. For :math:`m < p`, :data:`Xs` is :math:`(n - p + m)`-by-:math:`m`,

   :data:`ys` is :math:`(n - p + m)`-by-1, and :data:`M` is additional :math:`n`-by-:math:`p` workspace,

   which is required to recover the original solution vector after the system has been

   solved (see :func:`gsl_multifit_linear_genform2`). Optional data weights may be supplied in the

   vector :data:`w` of length :math:`n`, where :math:`W = \diag(w)`.

.. function:: int gsl_multifit_linear_solve (const double lambda, const gsl_matrix * Xs, const gsl_vector * ys, gsl_vector * cs, double * rnorm, double * snorm, gsl_multifit_linear_workspace * work)

   This function computes the regularized best-fit parameters :math:`\tilde{c}`

   which minimize the cost function

   :math:`\chi^2 = || \tilde{y} - \tilde{X} \tilde{c} ||^2 + \lambda^2 || \tilde{c} ||^2`

   which is in standard form. The least squares system must therefore be converted

   to standard form prior to calling this function.

   The observation vector :math:`\tilde{y}` is provided in :data:`ys` and the matrix of

   predictor variables :math:`\tilde{X}` in :data:`Xs`. The solution vector :math:`\tilde{c}` is

   returned in :data:`cs`, which has length min(:math:`m,p`). The SVD of :data:`Xs` must be computed prior

   to calling this function, using :func:`gsl_multifit_linear_svd`.

   The regularization parameter :math:`\lambda` is provided in :data:`lambda`.

   The residual norm :math:`|| \tilde{y} - \tilde{X} \tilde{c} || = ||y - X c||_W`

   is returned in :data:`rnorm`.

   The solution norm :math:`|| \tilde{c} || = ||L c||` is returned in

   :data:`snorm`.

.. function:: int gsl_multifit_linear_genform1 (const gsl_vector * L, const gsl_vector * cs, gsl_vector * c, gsl_multifit_linear_workspace * work)

   After a regularized system has been solved with

   :math:`L = \diag(\l_0,\l_1,...,\l_{p-1})`,

   this function backtransforms the standard form solution vector :data:`cs`

   to recover the solution vector of the original problem :data:`c`. The

   diagonal matrix elements :math:`l_i` are provided in

   the vector :data:`L`. It is allowed to have :data:`c` = :data:`cs` for an

   in-place transform.

.. function:: int gsl_multifit_linear_genform2 (const gsl_matrix * LQR, const gsl_vector * Ltau, const gsl_matrix * X, const gsl_vector * y, const gsl_vector * cs, const gsl_matrix * M, gsl_vector * c, gsl_multifit_linear_workspace * work)

              int gsl_multifit_linear_wgenform2 (const gsl_matrix * LQR, const gsl_vector * Ltau, const gsl_matrix * X, const gsl_vector * w, const gsl_vector * y, const gsl_vector * cs, const gsl_matrix * M, gsl_vector * c, gsl_multifit_linear_workspace * work)

   After a regularized system has been solved with a general rectangular matrix :math:`L`,

   specified by (:data:`LQR`, :data:`Ltau`), this function backtransforms the standard form solution :data:`cs`

   to recover the solution vector of the original problem, which is stored in :data:`c`,

   of length :math:`p`. The original least squares matrix and observation vector are provided in

   :data:`X` and :data:`y` respectively. :data:`M` is the matrix computed by

   :func:`gsl_multifit_linear_stdform2`. For weighted fits, the weight vector

   :data:`w` must also be supplied.

.. function:: int gsl_multifit_linear_applyW (const gsl_matrix * X, const gsl_vector * w, const gsl_vector * y, gsl_matrix * WX, gsl_vector * Wy)

   For weighted least squares systems with :math:`L = I`, this function may be used to

   convert the system to standard form by applying the weight matrix :math:`W = \diag(w)`

   to the least squares matrix :data:`X` and observation vector :data:`y`. On output, :data:`WX`

   is equal to :math:`W^{1/2} X` and :data:`Wy` is equal to :math:`W^{1/2} y`. It is allowed

   for :data:`WX` = :data:`X` and :data:`Wy` = :data:`y` for an in-place transform.

.. function:: int gsl_multifit_linear_lcurve (const gsl_vector * y, gsl_vector * reg_param, gsl_vector * rho, gsl_vector * eta, gsl_multifit_linear_workspace * work)

   This function computes the L-curve for a least squares system

   using the right hand side vector :data:`y` and the SVD decomposition

   of the least squares matrix :data:`X`, which must be provided

   to :func:`gsl_multifit_linear_svd` prior to

   calling this function. The output vectors :data:`reg_param`,

   :data:`rho`, and :data:`eta` must all be the same size, and will

   contain the regularization parameters :math:`\lambda_i`, residual norms

   :math:`||y - X c_i||`, and solution norms :math:`|| L c_i ||`

   which compose the L-curve, where :math:`c_i` is the regularized

   solution vector corresponding to :math:`\lambda_i`.

   The user may determine the number of points on the L-curve by

   adjusting the size of these input arrays. The regularization

   parameters :math:`\lambda_i` are estimated from the singular values

   of :data:`X`, and chosen to represent the most relevant portion of

   the L-curve.

.. function:: int gsl_multifit_linear_lcorner (const gsl_vector * rho, const gsl_vector * eta, size_t * idx)

   This function attempts to locate the corner of the L-curve

   :math:`(||y - X c||, ||L c||)` defined by the :data:`rho` and :data:`eta`

   input arrays respectively. The corner is defined as the point of maximum

   curvature of the L-curve in log-log scale. The :data:`rho` and :data:`eta`

   arrays can be outputs of :func:`gsl_multifit_linear_lcurve`. The

   algorithm used simply fits a circle to 3 consecutive points on the L-curve

   and uses the circle's radius to determine the curvature at

   the middle point. Therefore, the input array sizes must be

   :math:`\ge 3`. With more points provided for the L-curve, a better

   estimate of the curvature can be obtained. The array index

   corresponding to maximum curvature (ie: the corner) is returned

   in :data:`idx`. If the input arrays contain colinear points,

   this function could fail and return :macro:`GSL_EINVAL`.

.. function:: int gsl_multifit_linear_lcorner2 (const gsl_vector * reg_param, const gsl_vector * eta, size_t * idx)

   This function attempts to locate the corner of an alternate L-curve

   :math:`(\lambda^2, ||L c||^2)` studied by Rezghi and Hosseini, 2009.

   This alternate L-curve can provide better estimates of the

   regularization parameter for smooth solution vectors. The regularization

   parameters :math:`\lambda` and solution norms :math:`||L c||` are provided

   in the :data:`reg_param` and :data:`eta` input arrays respectively. The

   corner is defined as the point of maximum curvature of this

   alternate L-curve in linear scale. The :data:`reg_param` and :data:`eta`

   arrays can be outputs of :func:`gsl_multifit_linear_lcurve`. The

   algorithm used simply fits a circle to 3 consecutive points on the L-curve

   and uses the circle's radius to determine the curvature at

   the middle point. Therefore, the input array sizes must be

   :math:`\ge 3`. With more points provided for the L-curve, a better

   estimate of the curvature can be obtained. The array index

   corresponding to maximum curvature (ie: the corner) is returned

   in :data:`idx`. If the input arrays contain colinear points,

   this function could fail and return :macro:`GSL_EINVAL`.

.. function:: int gsl_multifit_linear_gcv_init(const gsl_vector * y, gsl_vector * reg_param, gsl_vector * UTy, double * delta0, gsl_multifit_linear_workspace * work)

   This function performs some initialization in preparation for computing

   the GCV curve and its minimum. The right hand side vector is provided

   in :data:`y`. On output, :data:`reg_param` is set to a vector of regularization

   parameters in decreasing order and may be of any size. The vector

   :data:`UTy` of size :math:`p` is set to :math:`U^T y`. The parameter

   :data:`delta0` is needed for subsequent steps of the GCV calculation.

.. function:: int gsl_multifit_linear_gcv_curve(const gsl_vector * reg_param, const gsl_vector * UTy, const double delta0, gsl_vector * G, gsl_multifit_linear_workspace * work)

   This funtion calculates the GCV curve :math:`G(\lambda)` and stores it in

   :data:`G` on output, which must be the same size as :data:`reg_param`. The

   inputs :data:`reg_param`, :data:`UTy` and :data:`delta0` are computed in

   :func:`gsl_multifit_linear_gcv_init`.

.. function:: int gsl_multifit_linear_gcv_min(const gsl_vector * reg_param, const gsl_vector * UTy, const gsl_vector * G, const double delta0, double * lambda, gsl_multifit_linear_workspace * work)

   This function computes the value of the regularization parameter

   which minimizes the GCV curve :math:`G(\lambda)` and stores it in

   :data:`lambda`. The input :data:`G` is calculated by

   :func:`gsl_multifit_linear_gcv_curve` and the inputs

   :data:`reg_param`, :data:`UTy` and :data:`delta0` are computed by

   :func:`gsl_multifit_linear_gcv_init`.

.. function:: double gsl_multifit_linear_gcv_calc(const double lambda, const gsl_vector * UTy, const double delta0, gsl_multifit_linear_workspace * work)

   This function returns the value of the GCV curve :math:`G(\lambda)` corresponding

   to the input :data:`lambda`.

.. function:: int gsl_multifit_linear_gcv(const gsl_vector * y, gsl_vector * reg_param, gsl_vector * G, double * lambda, double * G_lambda, gsl_multifit_linear_workspace * work)

   This function combines the steps :code:`gcv_init`, :code:`gcv_curve`,

   and :code:`gcv_min` defined above into a single function. The input

   :data:`y` is the right hand side vector. On output, :data:`reg_param` and

   :data:`G`, which must be the same size, are set to vectors of

   :math:`\lambda` and :math:`G(\lambda)` values respectively. The

   output :data:`lambda` is set to the optimal value of :math:`\lambda`

   which minimizes the GCV curve. The minimum value of the GCV curve is

   returned in :data:`G_lambda`.

.. function:: int gsl_multifit_linear_Lk (const size_t p, const size_t k, gsl_matrix * L)

   This function computes the discrete approximation to the derivative operator :math:`L_k` of

   order :data:`k` on a regular grid of :data:`p` points and stores it in :data:`L`. The dimensions of :data:`L` are

   :math:`(p-k)`-by-:math:`p`.

.. function:: int gsl_multifit_linear_Lsobolev (const size_t p, const size_t kmax, const gsl_vector * alpha, gsl_matrix * L, gsl_multifit_linear_workspace * work)

   This function computes the regularization matrix :data:`L` corresponding to the

   weighted Sobolov norm

   :math:`||L c||^2 = \sum_k \alpha_k^2 ||L_k c||^2` where :math:`L_k` approximates

   the derivative operator of order :math:`k`. This regularization norm can be useful

   in applications where it is necessary to smooth several derivatives of the solution.

   :data:`p` is the number of model parameters, :data:`kmax` is the highest derivative

   to include in the summation above, and :data:`alpha` is the vector of weights of

   size :data:`kmax` + 1, where :code:`alpha[k]` = :math:`\alpha_k` is the weight

   assigned to the derivative of order :math:`k`.  The output matrix :data:`L` is size

   :data:`p`-by-:data:`p` and upper triangular.

.. function:: double gsl_multifit_linear_rcond (const gsl_multifit_linear_workspace * work)

   This function returns the reciprocal condition number of the least squares matrix :math:`X`,

   defined as the ratio of the smallest and largest singular values,

   rcond = :math:`\sigma_{min}/\sigma_{max}`.

   The routine :func:`gsl_multifit_linear_svd` must first be called to compute the SVD

   of :math:`X`.

.. index::

   single: robust regression

   single: regression, robust

   single: least squares, robust

Robust linear regression

========================

Ordinary least squares (OLS) models are often heavily influenced by the presence of outliers.

Outliers are data points which do not follow the general trend of the other observations,

although there is strictly no precise definition of an outlier. Robust linear regression

refers to regression algorithms which are robust to outliers. The most common type of

robust regression is M-estimation. The general M-estimator minimizes the objective function

.. math:: \sum_i \rho(e_i) = \sum_i \rho (y_i - Y(c, x_i))

where :math:`e_i = y_i - Y(c, x_i)` is the residual of the ith data point, and

:math:`\rho(e_i)` is a function which should have the following properties:

* :math:`\rho(e) \ge 0`

* :math:`\rho(0) = 0`

* :math:`\rho(-e) = \rho(e)`

* :math:`\rho(e_1) > \rho(e_2)` for :math:`|e_1| > |e_2|`

The special case of ordinary least squares is given by :math:`\rho(e_i) = e_i^2`.

Letting :math:`\psi = \rho'` be the derivative of :math:`\rho`, differentiating

the objective function with respect to the coefficients :math:`c`

and setting the partial derivatives to zero produces the system of equations

.. math:: \sum_i \psi(e_i) X_i = 0

where :math:`X_i` is a vector containing row :math:`i` of the design matrix :math:`X`.

Next, we define a weight function :math:`w(e) = \psi(e)/e`, and let

:math:`w_i = w(e_i)`:

.. math:: \sum_i w_i e_i X_i = 0

This system of equations is equivalent to solving a weighted ordinary least squares

problem, minimizing :math:`\chi^2 = \sum_i w_i e_i^2`. The weights however, depend

on the residuals :math:`e_i`, which depend on the coefficients :math:`c`, which depend

on the weights. Therefore, an iterative solution is used, called Iteratively Reweighted

Least Squares (IRLS).

#. Compute initial estimates of the coefficients :math:`c^{(0)}` using ordinary least squares

#. For iteration :math:`k`, form the residuals :math:`e_i^{(k)} = (y_i - X_i c^{(k-1)})/(t \sigma^{(k)} \sqrt{1 - h_i})`,

   where :math:`t` is a tuning constant depending on the choice of :math:`\psi`, and :math:`h_i` are the

   statistical leverages (diagonal elements of the matrix :math:`X (X^T X)^{-1} X^T`). Including :math:`t`

   and :math:`h_i` in the residual calculation has been shown to improve the convergence of the method.

   The residual standard deviation is approximated as :math:`\sigma^{(k)} = MAD / 0.6745`, where MAD is the

   Median-Absolute-Deviation of the :math:`n-p` largest residuals from the previous iteration.

#. Compute new weights :math:`w_i^{(k)} = \psi(e_i^{(k)})/e_i^{(k)}`.

#. Compute new coefficients :math:`c^{(k)}` by solving the weighted least squares problem with

   weights :math:`w_i^{(k)}`.

#. Steps 2 through 4 are iterated until the coefficients converge or until some maximum iteration

   limit is reached. Coefficients are tested for convergence using the critera:

  .. only:: not texinfo

     .. math:: |c_i^{(k)} - c_i^{(k-1)}| \le \epsilon \times \hbox{max}(|c_i^{(k)}|, |c_i^{(k-1)}|)

  .. only:: texinfo

     ::

        |c_i^(k) - c_i^(k-1)| <= \epsilon * max(|c_i^(k)|, |c_i^(k-1)|)

  for all :math:`0 \le i < p` where :math:`\epsilon` is a small tolerance factor.

The key to this method lies in selecting the function :math:`\psi(e_i)` to assign

smaller weights to large residuals, and larger weights to smaller residuals. As

the iteration proceeds, outliers are assigned smaller and smaller weights, eventually

having very little or no effect on the fitted model.

.. type:: gsl_multifit_robust_workspace

   This workspace is used for robust least squares fitting.

.. function:: gsl_multifit_robust_workspace * gsl_multifit_robust_alloc (const gsl_multifit_robust_type * T, const size_t n, const size_t p)

   This function allocates a workspace for fitting a model to :data:`n`

   observations using :data:`p` parameters. The size of the workspace

   is :math:`O(np + p^2)`. The type :data:`T` specifies the

   function :math:`\psi` and can be selected from the following choices.

   .. type:: gsl_multifit_robust_type

      .. var:: gsl_multifit_robust_default

         This specifies the :data:`gsl_multifit_robust_bisquare` type (see below) and is a good

         general purpose choice for robust regression.

      .. var:: gsl_multifit_robust_bisquare

         This is Tukey's biweight (bisquare) function and is a good general purpose choice for

         robust regression. The weight function is given by

         .. only:: not texinfo

            .. math::

               w(e) =

               \left\{

                 \begin{array}{cc}

                   (1 - e^2)^2, & |e| \le 1 \\