File: Regularized-regression.html

package info (click to toggle)
gsl-ref-html 2.3-1
links: PTS
area: non-free
in suites: bookworm, bullseye, buster, forky, sid, trixie
size: 6,876 kB
ctags: 4,574
sloc: makefile: 35
file content (430 lines) | stat: -rw-r--r-- 29,600 bytes
<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/html4/loose.dtd">
<html>
<!-- Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016 The GSL Team.

Permission is granted to copy, distribute and/or modify this document
under the terms of the GNU Free Documentation License, Version 1.3 or
any later version published by the Free Software Foundation; with the
Invariant Sections being "GNU General Public License" and "Free Software
Needs Free Documentation", the Front-Cover text being "A GNU Manual",
and with the Back-Cover Text being (a) (see below). A copy of the
license is included in the section entitled "GNU Free Documentation
License".

(a) The Back-Cover Text is: "You have the freedom to copy and modify this
GNU Manual." -->
<!-- Created by GNU Texinfo 5.1, http://www.gnu.org/software/texinfo/ -->
<head>
<title>GNU Scientific Library &ndash; Reference Manual: Regularized regression</title>

<meta name="description" content="GNU Scientific Library &ndash; Reference Manual: Regularized regression">
<meta name="keywords" content="GNU Scientific Library &ndash; Reference Manual: Regularized regression">
<meta name="resource-type" content="document">
<meta name="distribution" content="global">
<meta name="Generator" content="makeinfo">
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<link href="index.html#Top" rel="start" title="Top">
<link href="Function-Index.html#Function-Index" rel="index" title="Function Index">
<link href="Least_002dSquares-Fitting.html#Least_002dSquares-Fitting" rel="up" title="Least-Squares Fitting">
<link href="Robust-linear-regression.html#Robust-linear-regression" rel="next" title="Robust linear regression">
<link href="Multi_002dparameter-regression.html#Multi_002dparameter-regression" rel="previous" title="Multi-parameter regression">
<style type="text/css">
<!--
a.summary-letter {text-decoration: none}
blockquote.smallquotation {font-size: smaller}
div.display {margin-left: 3.2em}
div.example {margin-left: 3.2em}
div.indentedblock {margin-left: 3.2em}
div.lisp {margin-left: 3.2em}
div.smalldisplay {margin-left: 3.2em}
div.smallexample {margin-left: 3.2em}
div.smallindentedblock {margin-left: 3.2em; font-size: smaller}
div.smalllisp {margin-left: 3.2em}
kbd {font-style:oblique}
pre.display {font-family: inherit}
pre.format {font-family: inherit}
pre.menu-comment {font-family: serif}
pre.menu-preformatted {font-family: serif}
pre.smalldisplay {font-family: inherit; font-size: smaller}
pre.smallexample {font-size: smaller}
pre.smallformat {font-family: inherit; font-size: smaller}
pre.smalllisp {font-size: smaller}
span.nocodebreak {white-space:nowrap}
span.nolinebreak {white-space:nowrap}
span.roman {font-family:serif; font-weight:normal}
span.sansserif {font-family:sans-serif; font-weight:normal}
ul.no-bullet {list-style: none}
-->
</style>


</head>

<body lang="en" bgcolor="#FFFFFF" text="#000000" link="#0000FF" vlink="#800080" alink="#FF0000">
<a name="Regularized-regression"></a>
<div class="header">
<p>
Next: <a href="Robust-linear-regression.html#Robust-linear-regression" accesskey="n" rel="next">Robust linear regression</a>, Previous: <a href="Multi_002dparameter-regression.html#Multi_002dparameter-regression" accesskey="p" rel="previous">Multi-parameter regression</a>, Up: <a href="Least_002dSquares-Fitting.html#Least_002dSquares-Fitting" accesskey="u" rel="up">Least-Squares Fitting</a> &nbsp; [<a href="Function-Index.html#Function-Index" title="Index" rel="index">Index</a>]</p>
</div>
<hr>
<a name="Regularized-regression-1"></a>
<h3 class="section">38.4 Regularized regression</h3>
<a name="index-ridge-regression"></a>
<a name="index-Tikhonov-regression"></a>
<a name="index-regression_002c-ridge"></a>
<a name="index-regression_002c-Tikhonov"></a>
<a name="index-least-squares_002c-regularized"></a>

<p>Ordinary weighted least squares models seek a solution vector <em>c</em>
which minimizes the residual
</p><div class="example">
<pre class="example">\chi^2 = || y - Xc ||_W^2
</pre></div>
<p>where <em>y</em> is the <em>n</em>-by-<em>1</em> observation vector,
<em>X</em> is the <em>n</em>-by-<em>p</em> design matrix, <em>c</em> is
the <em>p</em>-by-<em>1</em> solution vector,
<em>W =</em> diag<em>(w_1,...,w_n)</em> is the data weighting matrix,
and <em>||r||_W^2 = r^T W r</em>.
In cases where the least squares matrix <em>X</em> is ill-conditioned,
small perturbations (ie: noise) in the observation vector could lead to
widely different solution vectors <em>c</em>.
One way of dealing with ill-conditioned matrices is to use a &ldquo;truncated SVD&rdquo;
in which small singular values, below some given tolerance, are discarded
from the solution. The truncated SVD method is available using the functions
<code>gsl_multifit_linear_tsvd</code> and <code>gsl_multifit_wlinear_tsvd</code>. Another way
to help solve ill-posed problems is to include a regularization term in the least squares
minimization
</p><div class="example">
<pre class="example">\chi^2 = || y - Xc ||_W^2 + \lambda^2 || L c ||^2
</pre></div>
<p>for a suitably chosen regularization parameter <em>\lambda</em> and
matrix <em>L</em>. This type of regularization is known as Tikhonov, or ridge,
regression. In some applications, <em>L</em> is chosen as the identity matrix, giving
preference to solution vectors <em>c</em> with smaller norms.
Including this regularization term leads to the explicit &ldquo;normal equations&rdquo; solution
</p><div class="example">
<pre class="example">c = ( X^T W X + \lambda^2 L^T L )^-1 X^T W y
</pre></div>
<p>which reduces to the ordinary least squares solution when <em>L = 0</em>.
In practice, it is often advantageous to transform a regularized least
squares system into the form
</p><div class="example">
<pre class="example">\chi^2 = || y~ - X~ c~ ||^2 + \lambda^2 || c~ ||^2
</pre></div>
<p>This is known as the Tikhonov &ldquo;standard form&rdquo; and has the normal equations solution
<em>\tilde{c} = \left( \tilde{X}^T \tilde{X} + \lambda^2 I \right)^{-1} \tilde{X}^T \tilde{y}</em>.
For an <em>m</em>-by-<em>p</em> matrix <em>L</em> which is full rank and has <em>m &gt;= p</em> (ie: <em>L</em> is
square or has more rows than columns), we can calculate the &ldquo;thin&rdquo; QR decomposition of <em>L</em>, and
note that <em>||L c||</em> = <em>||R c||</em> since the <em>Q</em> factor will not change the norm. Since
<em>R</em> is <em>p</em>-by-<em>p</em>, we can then use the transformation
</p><div class="example">
<pre class="example">X~ = sqrt(W) X R^-1
y~ = sqrt(W) y
c~ = R c
</pre></div>
<p>to achieve the standard form. For a rectangular matrix <em>L</em> with <em>m &lt; p</em>,
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 <em>\tilde{X}</em>.
The matrix <em>L</em> is often chosen as the identity matrix, or as a first
or second finite difference operator, to ensure a smoothly varying
coefficient vector <em>c</em>, or as a diagonal matrix to selectively damp
each model parameter differently. If <em>L \ne I</em>, the user must first
convert the least squares problem to standard form using
<code>gsl_multifit_linear_stdform1</code> or <code>gsl_multifit_linear_stdform2</code>,
solve the system, and then backtransform the solution vector to recover
the solution of the original problem (see
<code>gsl_multifit_linear_genform1</code> and <code>gsl_multifit_linear_genform2</code>).
</p>
<p>In many regularization problems, care must be taken when choosing
the regularization parameter <em>\lambda</em>. Since both the
residual norm <em>||y - X c||</em> and solution norm <em>||L c||</em>
are being minimized, the parameter <em>\lambda</em> 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 <em>||y - X c||</em>
on the horizontal axis and the solution norm <em>||L c||</em> on the
vertical axis. This curve nearly always as an <em>L</em> 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.
</p>
<p>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 <em>y_i</em> is left out of the
right hand side, the resulting regularized solution should predict this element
accurately. This leads to choosing the parameter <em>\lambda</em>
which minimizes the GCV function
</p>
<div class="example">
<pre class="example">G(\lambda) = (||y - X c_{\lambda}||^2) / Tr(I_n - X X^I)^2
</pre></div>

<p>where <em>X_{\lambda}^I</em> is the matrix which relates the solution <em>c_{\lambda}</em>
to the right hand side <em>y</em>, ie: <em>c_{\lambda} = X_{\lambda}^I y</em>. GSL
provides routines to compute the GCV curve and its minimum.
</p>
<p>For most applications, the steps required to solve a regularized least
squares problem are as follows:
</p>
<ol>
<li> Construct the least squares system (<em>X</em>, <em>y</em>, <em>W</em>, <em>L</em>)

</li><li> Transform the system to standard form (<em>\tilde{X}</em>,<em>\tilde{y}</em>). This
step can be skipped if <em>L = I_p</em> and <em>W = I_n</em>.

</li><li> Calculate the SVD of <em>\tilde{X}</em>.

</li><li> Determine an appropriate regularization parameter <em>\lambda</em> (using for example
L-curve or GCV analysis).

</li><li> Solve the standard form system using the chosen <em>\lambda</em> and the SVD of <em>\tilde{X}</em>.

</li><li> Backtransform the standard form solution <em>\tilde{c}</em> to recover the
original solution vector <em>c</em>.

</li></ol>

<dl>
<dt><a name="index-gsl_005fmultifit_005flinear_005fstdform1"></a>Function: <em>int</em> <strong>gsl_multifit_linear_stdform1</strong> <em>(const gsl_vector * <var>L</var>, const gsl_matrix * <var>X</var>, const gsl_vector * <var>y</var>, gsl_matrix * <var>Xs</var>, gsl_vector * <var>ys</var>, gsl_multifit_linear_workspace * <var>work</var>)</em></dt>
<dt><a name="index-gsl_005fmultifit_005flinear_005fwstdform1"></a>Function: <em>int</em> <strong>gsl_multifit_linear_wstdform1</strong> <em>(const gsl_vector * <var>L</var>, const gsl_matrix * <var>X</var>, const gsl_vector * <var>w</var>, const gsl_vector * <var>y</var>, gsl_matrix * <var>Xs</var>, gsl_vector * <var>ys</var>, gsl_multifit_linear_workspace * <var>work</var>)</em></dt>
<dd><p>These functions define a regularization matrix
<em>L =</em> diag<em>(l_0,l_1,...,l_{p-1})</em>.
The diagonal matrix element <em>l_i</em> is provided by the
<em>i</em>th element of the input vector <var>L</var>.
The <em>n</em>-by-<em>p</em> least squares matrix <var>X</var> and
vector <var>y</var> of length <em>n</em> are then
converted to standard form as described above and the parameters
(<em>\tilde{X}</em>,<em>\tilde{y}</em>) are stored in <var>Xs</var> and <var>ys</var>
on output.  <var>Xs</var> and <var>ys</var> have the same dimensions as
<var>X</var> and <var>y</var>. Optional data weights may be supplied in the
vector <var>w</var> of length <em>n</em>. In order to apply this transformation,
<em>L^{-1}</em> must exist and so none of the <em>l_i</em>
may be zero. After the standard form system has been solved,
use <code>gsl_multifit_linear_genform1</code> to recover the original solution vector.
It is allowed to have <var>X</var> = <var>Xs</var> and <var>y</var> = <var>ys</var> for an in-place transform.
In order to perform a weighted regularized fit with <em>L = I</em>, the user may
call <code>gsl_multifit_linear_applyW</code> to convert to standard form.
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005fmultifit_005flinear_005fL_005fdecomp"></a>Function: <em>int</em> <strong>gsl_multifit_linear_L_decomp</strong> <em>(gsl_matrix * <var>L</var>, gsl_vector * <var>tau</var>)</em></dt>
<dd><p>This function factors the <em>m</em>-by-<em>p</em> regularization matrix
<var>L</var> into a form needed for the later transformation to standard form. <var>L</var>
may have any number of rows <em>m</em>. If <em>m \ge p</em> the QR decomposition of
<var>L</var> is computed and stored in <var>L</var> on output. If <em>m &lt; p</em>, the QR decomposition
of <em>L^T</em> is computed and stored in <var>L</var> on output. On output,
the Householder scalars are stored in the vector <var>tau</var> of size <em>MIN(m,p)</em>.
These outputs will be used by <code>gsl_multifit_linear_wstdform2</code> to complete the
transformation to standard form.
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005fmultifit_005flinear_005fstdform2"></a>Function: <em>int</em> <strong>gsl_multifit_linear_stdform2</strong> <em>(const gsl_matrix * <var>LQR</var>, const gsl_vector * <var>Ltau</var>, const gsl_matrix * <var>X</var>, const gsl_vector * <var>y</var>, gsl_matrix * <var>Xs</var>, gsl_vector * <var>ys</var>, gsl_matrix * <var>M</var>, gsl_multifit_linear_workspace * <var>work</var>)</em></dt>
<dt><a name="index-gsl_005fmultifit_005flinear_005fwstdform2"></a>Function: <em>int</em> <strong>gsl_multifit_linear_wstdform2</strong> <em>(const gsl_matrix * <var>LQR</var>, const gsl_vector * <var>Ltau</var>, const gsl_matrix * <var>X</var>, const gsl_vector * <var>w</var>, const gsl_vector * <var>y</var>, gsl_matrix * <var>Xs</var>, gsl_vector * <var>ys</var>, gsl_matrix * <var>M</var>, gsl_multifit_linear_workspace * <var>work</var>)</em></dt>
<dd><p>These functions convert the least squares system (<var>X</var>,<var>y</var>,<var>W</var>,<em>L</em>) to standard
form (<em>\tilde{X}</em>,<em>\tilde{y}</em>) which are stored in <var>Xs</var> and <var>ys</var>
respectively. The <em>m</em>-by-<em>p</em> regularization matrix <var>L</var> is specified by the inputs
<var>LQR</var> and <var>Ltau</var>, which are outputs from <code>gsl_multifit_linear_L_decomp</code>.
The dimensions of the standard form parameters (<em>\tilde{X}</em>,<em>\tilde{y}</em>)
depend on whether <em>m</em> is larger or less than <em>p</em>. For <em>m \ge p</em>,
<var>Xs</var> is <em>n</em>-by-<em>p</em>, <var>ys</var> is <em>n</em>-by-1, and <var>M</var> is
not used. For <em>m &lt; p</em>, <var>Xs</var> is <em>(n - p + m)</em>-by-<em>m</em>,
<var>ys</var> is <em>(n - p + m)</em>-by-1, and <var>M</var> is additional <em>n</em>-by-<em>p</em> workspace,
which is required to recover the original solution vector after the system has been
solved (see <code>gsl_multifit_linear_genform2</code>). Optional data weights may be supplied in the
vector <var>w</var> of length <em>n</em>, where <em>W =</em> diag(w).
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005fmultifit_005flinear_005fsolve"></a>Function: <em>int</em> <strong>gsl_multifit_linear_solve</strong> <em>(const double <var>lambda</var>, const gsl_matrix * <var>Xs</var>, const gsl_vector * <var>ys</var>, gsl_vector * <var>cs</var>, double * <var>rnorm</var>, double * <var>snorm</var>, gsl_multifit_linear_workspace * <var>work</var>)</em></dt>
<dd><p>This function computes the regularized best-fit parameters <em>\tilde{c}</em>
which minimize the cost function
<em>\chi^2 = || \tilde{y} - \tilde{X} \tilde{c} ||^2 + \lambda^2 || \tilde{c} ||^2</em> which is
in standard form. The least squares system must therefore be converted
to standard form prior to calling this function.
The observation vector <em>\tilde{y}</em> is provided in <var>ys</var> and the matrix of
predictor variables <em>\tilde{X}</em> in <var>Xs</var>. The solution vector <em>\tilde{c}</em> is
returned in <var>cs</var>, which has length min(<em>m,p</em>). The SVD of <var>Xs</var> must be computed prior
to calling this function, using <code>gsl_multifit_linear_svd</code>.
The regularization parameter <em>\lambda</em> is provided in <var>lambda</var>.
The residual norm <em>|| \tilde{y} - \tilde{X} \tilde{c} || = ||y - X c||_W</em> is returned in <var>rnorm</var>.
The solution norm <em>|| \tilde{c} || = ||L c||</em> is returned in
<var>snorm</var>.
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005fmultifit_005flinear_005fgenform1"></a>Function: <em>int</em> <strong>gsl_multifit_linear_genform1</strong> <em>(const gsl_vector * <var>L</var>, const gsl_vector * <var>cs</var>, gsl_vector * <var>c</var>, gsl_multifit_linear_workspace * <var>work</var>)</em></dt>
<dd><p>After a regularized system has been solved with
<em>L =</em> diag<em>(\l_0,\l_1,...,\l_{p-1})</em>,
this function backtransforms the standard form solution vector <var>cs</var>
to recover the solution vector of the original problem <var>c</var>. The
diagonal matrix elements <em>l_i</em> are provided in
the vector <var>L</var>. It is allowed to have <var>c</var> = <var>cs</var> for an
in-place transform.
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005fmultifit_005flinear_005fgenform2"></a>Function: <em>int</em> <strong>gsl_multifit_linear_genform2</strong> <em>(const gsl_matrix * <var>LQR</var>, const gsl_vector * <var>Ltau</var>, const gsl_matrix * <var>X</var>, const gsl_vector * <var>y</var>, const gsl_vector * <var>cs</var>, const gsl_matrix * <var>M</var>, gsl_vector * <var>c</var>, gsl_multifit_linear_workspace * <var>work</var>)</em></dt>
<dt><a name="index-gsl_005fmultifit_005flinear_005fwgenform2"></a>Function: <em>int</em> <strong>gsl_multifit_linear_wgenform2</strong> <em>(const gsl_matrix * <var>LQR</var>, const gsl_vector * <var>Ltau</var>, const gsl_matrix * <var>X</var>, const gsl_vector * <var>w</var>, const gsl_vector * <var>y</var>, const gsl_vector * <var>cs</var>, const gsl_matrix * <var>M</var>, gsl_vector * <var>c</var>, gsl_multifit_linear_workspace * <var>work</var>)</em></dt>
<dd><p>After a regularized system has been solved with a general rectangular matrix <em>L</em>,
specified by (<var>LQR</var>,<var>Ltau</var>), this function backtransforms the standard form solution <var>cs</var>
to recover the solution vector of the original problem, which is stored in <var>c</var>,
of length <em>p</em>. The original least squares matrix and observation vector are provided in
<var>X</var> and <var>y</var> respectively. <var>M</var> is the matrix computed by
<code>gsl_multifit_linear_stdform2</code>. For weighted fits, the weight vector
<var>w</var> must also be supplied.
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005fmultifit_005flinear_005fapplyW"></a>Function: <em>int</em> <strong>gsl_multifit_linear_applyW</strong> <em>(const gsl_matrix * <var>X</var>, const gsl_vector * <var>w</var>, const gsl_vector * <var>y</var>, gsl_matrix * <var>WX</var>, gsl_vector * <var>Wy</var>)</em></dt>
<dd><p>For weighted least squares systems with <em>L = I</em>, this function may be used to
convert the system to standard form by applying the weight matrix <em>W =</em> diag(<var>w</var>)
to the least squares matrix <var>X</var> and observation vector <var>y</var>. On output, <var>WX</var>
is equal to <em>W^{1/2} X</em> and <var>Wy</var> is equal to <em>W^{1/2} y</em>. It is allowed
for <var>WX</var> = <var>X</var> and <var>Wy</var> = <var>y</var> for an in-place transform.
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005fmultifit_005flinear_005flcurve"></a>Function: <em>int</em> <strong>gsl_multifit_linear_lcurve</strong> <em>(const gsl_vector * <var>y</var>, gsl_vector * <var>reg_param</var>, gsl_vector * <var>rho</var>, gsl_vector * <var>eta</var>, gsl_multifit_linear_workspace * <var>work</var>)</em></dt>
<dd><p>This function computes the L-curve for a least squares system
using the right hand side vector <var>y</var> and the SVD decomposition
of the least squares matrix <var>X</var>, which must be provided
to <code>gsl_multifit_linear_svd</code> prior to
calling this function. The output vectors <var>reg_param</var>,
<var>rho</var>, and <var>eta</var> must all be the same size, and will
contain the regularization parameters <em>\lambda_i</em>, residual norms
<em>||y - X c_i||</em>, and solution norms <em>|| L c_i ||</em>
which compose the L-curve, where <em>c_i</em> is the regularized
solution vector corresponding to <em>\lambda_i</em>.
The user may determine the number of points on the L-curve by
adjusting the size of these input arrays. The regularization
parameters <em>\lambda_i</em> are estimated from the singular values
of <var>X</var>, and chosen to represent the most relevant portion of
the L-curve.
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005fmultifit_005flinear_005flcorner"></a>Function: <em>int</em> <strong>gsl_multifit_linear_lcorner</strong> <em>(const gsl_vector * <var>rho</var>, const gsl_vector * <var>eta</var>, size_t * <var>idx</var>)</em></dt>
<dd><p>This function attempts to locate the corner of the L-curve
<em>(||y - X c||, ||L c||)</em> defined by the <var>rho</var> and <var>eta</var>
input arrays respectively. The corner is defined as the point of maximum
curvature of the L-curve in log-log scale. The <var>rho</var> and <var>eta</var>
arrays can be outputs of <code>gsl_multifit_linear_lcurve</code>. The
algorithm used simply fits a circle to 3 consecutive points on the L-curve
and uses the circle&rsquo;s radius to determine the curvature at
the middle point. Therefore, the input array sizes must be
<em>\ge 3</em>. 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 <var>idx</var>. If the input arrays contain colinear points,
this function could fail and return <code>GSL_EINVAL</code>.
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005fmultifit_005flinear_005flcorner2"></a>Function: <em>int</em> <strong>gsl_multifit_linear_lcorner2</strong> <em>(const gsl_vector * <var>reg_param</var>, const gsl_vector * <var>eta</var>, size_t * <var>idx</var>)</em></dt>
<dd><p>This function attempts to locate the corner of an alternate L-curve
<em>(\lambda^2, ||L c||^2)</em> 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 <em>\lambda</em> and solution norms <em>||L c||</em> are provided
in the <var>reg_param</var> and <var>eta</var> input arrays respectively. The
corner is defined as the point of maximum curvature of this
alternate L-curve in linear scale. The <var>reg_param</var> and <var>eta</var>
arrays can be outputs of <code>gsl_multifit_linear_lcurve</code>. The
algorithm used simply fits a circle to 3 consecutive points on the L-curve
and uses the circle&rsquo;s radius to determine the curvature at
the middle point. Therefore, the input array sizes must be
<em>\ge 3</em>. 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 <var>idx</var>. If the input arrays contain colinear points,
this function could fail and return <code>GSL_EINVAL</code>.
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005fmultifit_005flinear_005fgcv_005finit_0028const"></a>Function: <em>int</em> <strong>gsl_multifit_linear_gcv_init(const</strong> <em>gsl_vector * <var>y</var>, gsl_vector * <var>reg_param</var>, gsl_vector * <var>UTy</var>, double * <var>delta0</var>, gsl_multifit_linear_workspace * <var>work</var>)</em></dt>
<dd><p>This function performs some initialization in preparation for computing
the GCV curve and its minimum. The right hand side vector is provided
in <var>y</var>. On output, <var>reg_param</var> is set to a vector of regularization
parameters in decreasing order and may be of any size. The vector
<var>UTy</var> of size <em>p</em> is set to <em>U^T y</em>. The parameter
<var>delta0</var> is needed for subsequent steps of the GCV calculation.
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005fmultifit_005flinear_005fgcv_005fcurve_0028const"></a>Function: <em>int</em> <strong>gsl_multifit_linear_gcv_curve(const</strong> <em>gsl_vector * <var>reg_param</var>, const gsl_vector * <var>UTy</var>, const double <var>delta0</var>, gsl_vector * <var>G</var>, gsl_multifit_linear_workspace * <var>work</var>)</em></dt>
<dd><p>This funtion calculates the GCV curve <em>G(\lambda)</em> and stores it in
<var>G</var> on output, which must be the same size as <var>reg_param</var>. The
inputs <var>reg_param</var>, <var>UTy</var> and <var>delta0</var> are computed in
<code>gsl_multifit_linear_gcv_init</code>.
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005fmultifit_005flinear_005fgcv_005fmin_0028const"></a>Function: <em>int</em> <strong>gsl_multifit_linear_gcv_min(const</strong> <em>gsl_vector * <var>reg_param</var>, const gsl_vector * <var>UTy</var>, const gsl_vector * <var>G</var>, const double <var>delta0</var>, double * <var>lambda</var>, gsl_multifit_linear_workspace * <var>work</var>)</em></dt>
<dd><p>This function computes the value of the regularization parameter
which minimizes the GCV curve <em>G(\lambda)</em> and stores it in
<var>lambda</var>. The input <var>G</var> is calculated by
<code>gsl_multifit_linear_gcv_curve</code> and the inputs
<var>reg_param</var>, <var>UTy</var> and <var>delta0</var> are computed by
<code>gsl_multifit_linear_gcv_init</code>.
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005fmultifit_005flinear_005fgcv_005fcalc_0028const"></a>Function: <em>double</em> <strong>gsl_multifit_linear_gcv_calc(const</strong> <em>double <var>lambda</var>, const gsl_vector * <var>UTy</var>, const double <var>delta0</var>, gsl_multifit_linear_workspace * <var>work</var>)</em></dt>
<dd><p>This function returns the value of the GCV curve <em>G(\lambda)</em> corresponding
to the input <var>lambda</var>.
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005fmultifit_005flinear_005fgcv_0028const"></a>Function: <em>int</em> <strong>gsl_multifit_linear_gcv(const</strong> <em>gsl_vector * <var>y</var>, gsl_vector * <var>reg_param</var>, gsl_vector * <var>G</var>, double * <var>lambda</var>, double * <var>G_lambda</var>, gsl_multifit_linear_workspace * <var>work</var>)</em></dt>
<dd><p>This function combines the steps <code>gcv_init</code>, <code>gcv_curve</code>,
and <code>gcv_min</code> defined above into a single function. The input
<var>y</var> is the right hand side vector. On output, <var>reg_param</var> and
<var>G</var>, which must be the same size, are set to vectors of
<em>\lambda</em> and <em>G(\lambda)</em> values respectively. The
output <var>lambda</var> is set to the optimal value of <em>\lambda</em>
which minimizes the GCV curve. The minimum value of the GCV curve is
returned in <var>G_lambda</var>.
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005fmultifit_005flinear_005fLk"></a>Function: <em>int</em> <strong>gsl_multifit_linear_Lk</strong> <em>(const size_t <var>p</var>, const size_t <var>k</var>, gsl_matrix * <var>L</var>)</em></dt>
<dd><p>This function computes the discrete approximation to the derivative operator <em>L_k</em> of
order <var>k</var> on a regular grid of <var>p</var> points and stores it in <var>L</var>. The dimensions of <var>L</var> are
<em>(p-k)</em>-by-<em>p</em>.
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005fmultifit_005flinear_005fLsobolev"></a>Function: <em>int</em> <strong>gsl_multifit_linear_Lsobolev</strong> <em>(const size_t <var>p</var>, const size_t <var>kmax</var>, const gsl_vector * <var>alpha</var>, gsl_matrix * <var>L</var>, gsl_multifit_linear_workspace * <var>work</var>)</em></dt>
<dd><p>This function computes the regularization matrix <var>L</var> corresponding to the weighted Sobolov norm
<em>||L c||^2 = \sum_k \alpha_k^2 ||L_k c||^2</em> where <em>L_k</em> approximates the derivative
operator of order <em>k</em>. This regularization norm can be useful in applications where
it is necessary to smooth several derivatives of the solution. <var>p</var> is the number of
model parameters, <var>kmax</var> is the highest derivative to include in the summation above, and
<var>alpha</var> is the vector of weights of size <var>kmax</var> + 1, where <var>alpha</var>[k] = <em>\alpha_k</em>
is the weight assigned to the derivative of order <em>k</em>.  The output matrix <var>L</var> is size
<var>p</var>-by-<var>p</var> and upper triangular.
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005fmultifit_005flinear_005frcond"></a>Function: <em>double</em> <strong>gsl_multifit_linear_rcond</strong> <em>(const gsl_multifit_linear_workspace * <var>work</var>)</em></dt>
<dd><p>This function returns the reciprocal condition number of the least squares matrix <em>X</em>,
defined as the ratio of the smallest and largest singular values, rcond = <em>\sigma_{min}/\sigma_{max}</em>.
The routine <code>gsl_multifit_linear_svd</code> must first be called to compute the SVD of <em>X</em>.
</p></dd></dl>

<hr>
<div class="header">
<p>
Next: <a href="Robust-linear-regression.html#Robust-linear-regression" accesskey="n" rel="next">Robust linear regression</a>, Previous: <a href="Multi_002dparameter-regression.html#Multi_002dparameter-regression" accesskey="p" rel="previous">Multi-parameter regression</a>, Up: <a href="Least_002dSquares-Fitting.html#Least_002dSquares-Fitting" accesskey="u" rel="up">Least-Squares Fitting</a> &nbsp; [<a href="Function-Index.html#Function-Index" title="Index" rel="index">Index</a>]</p>
</div>



</body>
</html>