### 38.4 Regularized regression

Ordinary weighted least squares models seek a solution vector c which minimizes the residual

    \chi^2 = || y - Xc ||_W^2

where y is the n-by-1 observation vector, X is the n-by-p design matrix, c is the p-by-1 solution vector, W = diag(w_1,...,w_n) is the data weighting matrix, and ||r||_W^2 = r^T W r. In cases where the least squares matrix X is ill-conditioned, small perturbations (ie: noise) in the observation vector could lead to widely different solution vectors 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 `gsl_multifit_linear_tsvd` and `gsl_multifit_wlinear_tsvd`. Another way to help solve ill-posed problems is to include a regularization term in the least squares minimization

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

for a suitably chosen regularization parameter \lambda and matrix L. This type of regularization is known as Tikhonov, or ridge, regression. In some applications, L is chosen as the identity matrix, giving preference to solution vectors c with smaller norms. Including this regularization term leads to the explicit “normal equations” solution

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

which reduces to the ordinary least squares solution when L = 0. In practice, it is often advantageous to transform a regularized least squares system into the form

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

This is known as the Tikhonov “standard form” and has the normal equations solution \tilde{c} = \left( \tilde{X}^T \tilde{X} + \lambda^2 I \right)^{-1} \tilde{X}^T \tilde{y}. For an m-by-p matrix L which is full rank and has m >= p (ie: L is square or has more rows than columns), we can calculate the “thin” QR decomposition of L, and note that ||L c|| = ||R c|| since the Q factor will not change the norm. Since R is p-by-p, we can then use the transformation

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

to achieve the standard form. For a rectangular matrix L with 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 \tilde{X}. The matrix L is often chosen as the identity matrix, or as a first or second finite difference operator, to ensure a smoothly varying coefficient vector c, or as a diagonal matrix to selectively damp each model parameter differently. If L \ne I, the user must first convert the least squares problem to standard form using `gsl_multifit_linear_stdform1` or `gsl_multifit_linear_stdform2`, solve the system, and then backtransform the solution vector to recover the solution of the original problem (see `gsl_multifit_linear_genform1` and `gsl_multifit_linear_genform2`).

In many regularization problems, care must be taken when choosing the regularization parameter \lambda. Since both the residual norm ||y - X c|| and solution norm ||L c|| are being minimized, the parameter \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 ||y - X c|| on the horizontal axis and the solution norm ||L c|| on the vertical axis. This curve nearly always as an 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 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 \lambda which minimizes the GCV function

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

where X_{\lambda}^I is the matrix which relates the solution c_{\lambda} to the right hand side y, ie: 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:

1. Construct the least squares system (X, y, W, L)
2. Transform the system to standard form (\tilde{X},\tilde{y}). This step can be skipped if L = I_p and W = I_n.
3. Calculate the SVD of \tilde{X}.
4. Determine an appropriate regularization parameter \lambda (using for example L-curve or GCV analysis).
5. Solve the standard form system using the chosen \lambda and the SVD of \tilde{X}.
6. Backtransform the standard form solution \tilde{c} to recover the original solution vector 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)
Function: 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 L = diag(l_0,l_1,...,l_{p-1}). The diagonal matrix element l_i is provided by the ith element of the input vector L. The n-by-p least squares matrix X and vector y of length n are then converted to standard form as described above and the parameters (\tilde{X},\tilde{y}) are stored in Xs and ys on output. Xs and ys have the same dimensions as X and y. Optional data weights may be supplied in the vector w of length n. In order to apply this transformation, L^{-1} must exist and so none of the l_i may be zero. After the standard form system has been solved, use `gsl_multifit_linear_genform1` to recover the original solution vector. It is allowed to have X = Xs and y = ys for an in-place transform. In order to perform a weighted regularized fit with L = I, the user may call `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 m-by-p regularization matrix L into a form needed for the later transformation to standard form. L may have any number of rows m. If m \ge p the QR decomposition of L is computed and stored in L on output. If m < p, the QR decomposition of L^T is computed and stored in L on output. On output, the Householder scalars are stored in the vector tau of size MIN(m,p). These outputs will be used by `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)
Function: 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 (X,y,W,L) to standard form (\tilde{X},\tilde{y}) which are stored in Xs and ys respectively. The m-by-p regularization matrix L is specified by the inputs LQR and Ltau, which are outputs from `gsl_multifit_linear_L_decomp`. The dimensions of the standard form parameters (\tilde{X},\tilde{y}) depend on whether m is larger or less than p. For m \ge p, Xs is n-by-p, ys is n-by-1, and M is not used. For m < p, Xs is (n - p + m)-by-m, ys is (n - p + m)-by-1, and M is additional n-by-p workspace, which is required to recover the original solution vector after the system has been solved (see `gsl_multifit_linear_genform2`). Optional data weights may be supplied in the vector w of length n, where 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 \tilde{c} which minimize the cost function \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 \tilde{y} is provided in ys and the matrix of predictor variables \tilde{X} in Xs. The solution vector \tilde{c} is returned in cs, which has length min(m,p). The SVD of Xs must be computed prior to calling this function, using `gsl_multifit_linear_svd`. The regularization parameter \lambda is provided in lambda. The residual norm || \tilde{y} - \tilde{X} \tilde{c} || = ||y - X c||_W is returned in rnorm. The solution norm || \tilde{c} || = ||L c|| is returned in 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 L = diag(\l_0,\l_1,...,\l_{p-1}), this function backtransforms the standard form solution vector cs to recover the solution vector of the original problem c. The diagonal matrix elements l_i are provided in the vector L. It is allowed to have c = 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)
Function: 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 L, specified by (LQR,Ltau), this function backtransforms the standard form solution cs to recover the solution vector of the original problem, which is stored in c, of length p. The original least squares matrix and observation vector are provided in X and y respectively. M is the matrix computed by `gsl_multifit_linear_stdform2`. For weighted fits, the weight vector 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 L = I, this function may be used to convert the system to standard form by applying the weight matrix W = diag(w) to the least squares matrix X and observation vector y. On output, WX is equal to W^{1/2} X and Wy is equal to W^{1/2} y. It is allowed for WX = X and Wy = 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 y and the SVD decomposition of the least squares matrix X, which must be provided to `gsl_multifit_linear_svd` prior to calling this function. The output vectors reg_param, rho, and eta must all be the same size, and will contain the regularization parameters \lambda_i, residual norms ||y - X c_i||, and solution norms || L c_i || which compose the L-curve, where c_i is the regularized solution vector corresponding to \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 \lambda_i are estimated from the singular values of 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 (||y - X c||, ||L c||) defined by the rho and eta input arrays respectively. The corner is defined as the point of maximum curvature of the L-curve in log-log scale. The rho and eta arrays can be outputs of `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 \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 idx. If the input arrays contain colinear points, this function could fail and return `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 (\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 \lambda and solution norms ||L c|| are provided in the reg_param and eta input arrays respectively. The corner is defined as the point of maximum curvature of this alternate L-curve in linear scale. The reg_param and eta arrays can be outputs of `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 \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 idx. If the input arrays contain colinear points, this function could fail and return `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 y. On output, reg_param is set to a vector of regularization parameters in decreasing order and may be of any size. The vector UTy of size p is set to U^T y. The parameter 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 G(\lambda) and stores it in G on output, which must be the same size as reg_param. The inputs reg_param, UTy and delta0 are computed in `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 G(\lambda) and stores it in lambda. The input G is calculated by `gsl_multifit_linear_gcv_curve` and the inputs reg_param, UTy and delta0 are computed by `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 G(\lambda) corresponding to the input 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 `gcv_init`, `gcv_curve`, and `gcv_min` defined above into a single function. The input y is the right hand side vector. On output, reg_param and G, which must be the same size, are set to vectors of \lambda and G(\lambda) values respectively. The output lambda is set to the optimal value of \lambda which minimizes the GCV curve. The minimum value of the GCV curve is returned in 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 L_k of order k on a regular grid of p points and stores it in L. The dimensions of L are (p-k)-by-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 L corresponding to the weighted Sobolov norm ||L c||^2 = \sum_k \alpha_k^2 ||L_k c||^2 where L_k approximates the derivative operator of order k. This regularization norm can be useful in applications where it is necessary to smooth several derivatives of the solution. p is the number of model parameters, kmax is the highest derivative to include in the summation above, and alpha is the vector of weights of size kmax + 1, where alpha[k] = \alpha_k is the weight assigned to the derivative of order k. The output matrix L is size p-by-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 X, defined as the ratio of the smallest and largest singular values, rcond = \sigma_{min}/\sigma_{max}. The routine `gsl_multifit_linear_svd` must first be called to compute the SVD of X.