LARS Library: Least Angle Regression laSso Library

Frank Vanden Berghen

The PDF version of this page is available here.
You can download here the LARS library ( v1.06 - released 2/10/2007) (new binaries that works on any windows-matlab version).

Let's assume that we want to perform a simple identification task. Let's assume that you have a set of

measures $(x_{(1)},x_{(2)},\ldots,x_{(n)},y)$ . We want to find a function (or a model) that predicts the measure

in function of the measures $x_{(j)}, \; j=1,\ldots,n$ . The vector $x \in \Re^n = ( x_{(1)},\ldots,x_{(n)})^t$ has several names:

is the regressor vector.
is the input.
is the vector of independent variables.

has several names:

is the target.
is the output.
( is the dependant variable).

The pair

has several names:

is an input $\leftrightarrow$ output pair.
is a sample.

We want to find

such that

. Let's also assume that you know that

belongs to the family of linear models. I.e. $f(x)=\sum_{j=1}^n x_{(j)} \beta_j+p=<x,\beta>+p$ where $<\cdot,\cdot>$ is the dot product of two vectors and

is the constant term. You can also write $f(x)=x^t \beta+p$ where $\beta \in \Re^n$ and

are describing

. We will now assume, without loss of generality, that

. We will see later how to compte

if this is not the case. $\beta$ is the model that we want to find. Let's now assume that we have many input $\leftrightarrow$ output pairs. i.e. we have

$\displaystyle ( x_{(1)},\ldots, x_{(n)} )_{(1)} = X_{(1)}$	$\displaystyle \leftrightarrow$	$\displaystyle y_1$	(1)
$\displaystyle X_{(2)}$	$\displaystyle \leftrightarrow$	$\displaystyle y_2$
$\displaystyle \vdots$		$\displaystyle \vdots$
$\displaystyle X_{(m)}$	$\displaystyle \leftrightarrow$	$\displaystyle y_m$

We want to compute $\beta$ such that $\forall i \;\; X_{(i)}^t \beta = y_i$ . Using matrix notation, $\beta$ is the solution of:

$\displaystyle X \beta = y$

(2)

where $X \in \Re^{m \times n}$ is a matrix containing on each line a different regressor and $y \in\Re^m$ contains all the targets. The columns of

contain the independent variables or, in short, the variables. If $m \leq n$ , we can compute $\beta$ using a simple Gauss-Jordan elimination. That's not very interesting. We are usually more interested in the case where

: when the linear system $X \beta =y$ is over-determined. We want to find $\beta^*$ such that the Total Prediction Error (TPError) is minimum:

TPError $\displaystyle (\beta^*) = \min_{\beta} \left( \mbox{\em TPError}(\beta) \right)$

(3)

If we define the total prediction error TPError as

Sum of Squared Errors $\displaystyle (\beta) = \sum_{i=1}^m ( X_{(i)} \beta - y_i ) ^2 = ( X \beta - y )^2 = \Vert X \beta - y \Vert _2$

where $\Vert \cdot \Vert _2$ is the

-norm of a vector, we obtain the Classical Linear Regression or Least Min-Square(LS):

Classical Linear Regression: $\beta^*$ is the solution to: $\displaystyle \min_{\beta} \left( \Vert X \beta - y \Vert _2 \right) = \min_{\beta} \left( ( X \beta - y )^2 \right)$

Such a definition of the total prediction error TPError can give un-deserved weight to a small number of bad samples. It means that one sample with a large prediction error is enough to change completely $\beta^*$ (this is due to the squaring effect that ``amplifies'' the contribution of the ``outliers''). This is why some researcher use a

instead of a

:

Robust Regression: $\beta^*$ is the solution to: $\displaystyle \min_{\beta} \left( \Vert X \beta - y \Vert _1 \right) = \min_{\beta} \left( \sum_{i=1}^m \vert X_{(i)} \beta - y_i \vert \right)$

where $\vert\cdot\vert$ is the absolute value. We will not discuss here Robust regression algorithms. To find the solution of a classical Linear Regression, we have to solve:

$\displaystyle Error^2(\beta^*) =$	$\displaystyle \min_{\beta} \left\{ ( X \beta - y )^2 \right\}$
$\displaystyle =$	$\displaystyle \min_{\beta} \left\{ (X \beta - y)^t (X \beta - y) \right\}$
$\displaystyle =$	$\displaystyle \min_{\beta} \left\{ (\beta^t X^t - y^t) (X \beta - y) \right\}$
$\displaystyle =$	$\displaystyle \min_{\beta} \left\{ ( \beta^t X^t X \beta - y^t X \beta - \beta^t X^t y + y^t y) \right\} \nonumber$

Since $y^t X \beta$ is a scalar, we can take its transpose without changing its value: $y^t X \beta= (y^t X \beta)^t= \beta^t (y^t X)^t = \beta^t X^t y$ . We obtain:

$\displaystyle Error^2(\beta^*) =$

$\displaystyle \min_{\beta} \left\{ \beta^t X^t X \beta - 2 \beta^t X^t y \right\}$

$\displaystyle Error^2(\beta^*)$ is minimum	$\displaystyle \Leftrightarrow$	$\displaystyle \frac{\partial Error^2(\beta)}{\partial \beta}(\beta^*)=0$
	$\displaystyle \Leftrightarrow$	$\displaystyle 2 X^t X \beta^* - 2 X^t y = 0$
	$\displaystyle \Leftrightarrow$	$\displaystyle X^t X \beta^* = X^t y$	(4)

The equation 4 is the well-known normal equation. On the left-hand-side of this equation we find $X^t X = C \in \Re^{n \times n}$ : the correlation matrix.

is symmetric and semi-positive definite (all the eigenvalues of

are $\geq 0$ ). The element $C_{ij}= \sum_{k=1}^m X_{ik} X_{jk} = <X_i,X_j>= X_i^t X_j$ is the correlation between column

and column

(where

stands for column

of

). The right-hand-side of equation 4 is also interesting: it contains the univariate relation of all the columns of

(the variables) with the target

.

We did previously the assumption that the constant term

of the model is zero.

is null when

and

where $mean(\cdot)$ is the mean operator and

is the $i^{th}$ column of

. Thus, in the general case, when $p \neq 0$ , the procedure to build a model is:

Compute where $std(\cdot)$ is the ``standard deviation'' operator.
Normalize the target: $\displaystyle \tilde{y}=\frac{y-mean(y)}{std(y)}$
Normalize the columns of : $\displaystyle \tilde{X_i}=\frac{X_i-mean(X_i)}{std(X_i)}$
Since the columns of $\tilde{X}$ have zero mean, and since $\tilde{y}$ has zero mean, we have . We can thus use the normal equations: $\tilde{X}^t \tilde{X} \tilde{\beta} = \tilde{X^t} \tilde{y}$ to find $\tilde{\beta}$ . All the columns inside $\tilde{X}$ have a standard deviation of one. This increase numerical stability.
Currently, has been computed on the ``normalized columns'' and . Let's convert back the model so that it can be applied on the un-normalized and :
- $\displaystyle \beta_i= \frac{std(y)}{std(X_i)} \; \tilde{\beta}_i, \;\; i=1,\ldots,n$
- $\displaystyle p=mean(y)-\sum_{i=1}^n mean(X_i) \beta_i$

From now on, in the rest of this paper, we will always assume that

(data has been normalized).

One simple way to solve the normal equations 4 is the following:

Compute a Cholesky factorization of . i.e. Compute such that is a lower triangular matrix and such that .
Let's define $v:= L \beta$ . Then we have: $C \beta = L^t L \beta = L^t v = X^t y$ (the last equality comes from equation 4). This last equation ( )can be solved easily because is triangular. Once we have found , we solve $L \beta = v$ to find $\beta$ (this is also easy because is triangular).

Unfortunately, the above algorithm is not always working properly. Let's consider the case where

. i.e. the variables

and

are 100% linearly correlated together. In this case, we will have $\beta_i X_i+ \beta_j X_j=\beta_i X_i+\beta_j X_i=(\beta_i+\beta_j)X_i$ . $\beta _i$ can be any value at the condition that $\beta_j:=-\beta_i$ (so that we obtain $\beta_i+\beta_j=0$ and none of the columns

and

are meaningful to predict

). Thus, $\beta _i$ can become extremely large and give undeserved weight to the column

. We have one degree of liberty on the $\beta _i$ 's. This difficulty arises when:

A column of is a linear combination of other columns.
The rank of A is lower than .
The matrix is not of full rank .

This problem is named the ``multi-colinearity problem''. It is illustrated in figure 1 and 2.

**Figure 1:** Stable Model
$\begin{figure} \centering\epsfig{figure=stability1.eps, width=10cm, height=5.5cm} \end{figure}$

**Figure 2:** Unstable Models
$\begin{figure} \epsfig{figure=stability2.eps, width=8cm, height=5cm} \epsfig{figure=stability3.eps, width=8cm, height=5cm} \end{figure}$

The blue points inside figures 1 and 2 are the samples. The green points are the lines of $X \in \Re^{9 \times 2}$ (the inputs) and the heights of the blue points are the targets ( $y \in \Re^9$ )(the output). There is a linear correlation between

and

: $X_1 \approx C^{te} X_2$ . In the ``blue direction'' (perpendicular to the black dashed line), we do not have enough reliable information to decide if the output

should increase or decrease. Inside figures 1 and 2, a model $\beta$ is represented by a plane. In this situation, the ideal model that we want to obtain, is drawn in green in figure 1 (the plane/model is ''flat'': there is no un-due variation of

in the ``blue direction''). The equation of the ideal model is either:

or

or

. Note that, inside figure 2 (and also 1), we do NOT have exactly $X_1 = C^{te} X_2$ (the green point are not totally on the dashed line $x_1=\frac{b}{a}x_2$ ). Thus, we theoretically have enough information in the ``blue direction'' to construct a model. However exploiting this information to predict

is not a good idea and leads to unstable models: a small perturbation of the data (

or

) leads to a completely different model: Inside figure 2, we have represented in red two models based on small perturbations of

and

. Usually, unstable models have poor accuracy. The main question now is: How to stabilize a model? Here are two classical solutions:

Inside figure 1, the 'information' that is contained in is duplicated in . We can simply drop one of the two variables. We obtain the models or . The main difficulty here is to select carefully the columns that must be dropped. The solution is the ``LARS/Lasso Algorithm'' described later in this paper.
The other choice is to keep inside the model all the variables ( and ). However, in this case, we must be sure that:
- $\beta_1 \approx \beta_2$
- $\beta_1$ and $\beta_2$ are small values (I remind you that $\beta_1$ and $\beta_2$ can become arbitrarily large).
We must find an algorithm to compute that will ``equilibrate'' and ``reduce'' the 's:
$\displaystyle \beta^*$ is the solution of $\displaystyle = \min_{\beta} \left\{ ( X \beta - y )^2 + \lambda \beta^2 \right\}$

$\displaystyle \hspace{5cm} \Leftrightarrow (X^t X + \lambda I) \beta^* = X^t Y$ (5)

This technique is named ``Ridge regression''. When $\lambda$ increases, all the $\beta _i$ are ``equilibrated'' and ``pushed'' to zero as illustrated in figure 3.

Figure 3: The $\beta _i$ 's sorted from greatest to smallest for different value of $\lambda$ .
$\begin{figure} \centering\epsfig{figure=ridge1.eps, width=10cm, height=5.5cm} \end{figure}$

We will discuss later of a way to compute an optimal regularization parameter $\lambda$ . Inside figure 1, we obtain the model $\beta=(\beta_1 \;\; \beta_2)^t=(c\;\;c)^t$ and we have thus with ''small''.

The ``Ridge regression'' technique is particularly simple and efficient. It has a nice graphical explanation. The ``Ridge
regression'' technique is in fact searching for $\beta^*$ that minimizes $Q(\beta)=$ Squared Error $(\beta)=(X \beta-y)^2$ under the constraint that $\Vert\beta \Vert _2< r$ . This definition is illustrated on the left of figure 4. We seek the solution $\beta^*$ of the minimization problem:

$\displaystyle \min_{\beta \in \Re^n} Q(\beta) \equiv g_k^t \beta + \frac{1}{2} \beta^t H \beta \;\;\;$ subject to $\displaystyle \Vert \beta \Vert _2 < r$

(with

and

) To solve this minimization problem, we first have to rewrite the constraint: $c(\beta)=\frac{1}{2}r^2-\frac{1}{2}\beta^t \beta<0$ . Now, we introduce a Lagrange multiplier $\lambda$ and the Lagrange function:

$\displaystyle \L (\beta,\lambda)=Q(\beta)-\lambda c(\beta)$

Using the Lagrange First order optimality criterion, we know that the solution $(\beta^*,\lambda^*)$ to this optimization problem is given by:

$\displaystyle 0=\nabla_s \L (\beta^*,\lambda^*)= \nabla Q(\beta^*)-\lambda^* \... ...\beta^*)= H \beta^* + g + \lambda^* \beta^* = (H + \lambda^* I)\beta^* + g = 0$

(6)

which is equation 5.

**Figure 4:** On the left: Illustration of Ridge Regression. On the right: Illustration of Lasso
$\begin{figure} \epsfig{figure=ridge2.eps, width=8cm, height=7cm} \epsfig{figure=lasso1.eps, width=8cm, height=7cm} \end{figure}$

The constraint $\Vert\beta \Vert _2< r$ prevents the $\beta _i$ 's from becoming arbitrarily large. As illustrated in figure 3, the ridge regression algorithm ``pushes'' all the $\beta _i$ 's towards zero but none of the $\beta _i$ actually gets a null value. It's interesting to have $\beta_i=0$ because, in this case, we can completely remove from the model the column

. A model that is using less variables has many advantages: it is easier to interpret for a human and it's also more stable. What happens if we replace the constraint $\Vert\beta \Vert _2< r$ with the constraint $\Vert\beta \Vert _1 = \sum_{i=1}^n \vert\beta_i\vert < r$ ? See the right of illustration 4. You can see on the figure that $\beta_{lasso}=(\beta_1^*,\beta_2^*)$ with $\beta_1^*=0$ . You can see on figure 5 that, as you decrease

, more and more $\beta_i^*$ are going to zero.

**Figure 5:** The $\beta _i$ 's sorted from greatest to smallest for different values of
$\begin{figure} \centering\epsfig{figure=lasso2.eps, width=10cm, height=5.5cm} \end{figure}$

The Lasso regression is defined by: $\beta^*$ is the solution of the minimization problem:

$\displaystyle \min_{\beta \in \Re^n} (X \beta -y)^2 \;\;\;$ subject to $\displaystyle \Vert \beta \Vert _1 < r$

(7)

As you can see, computing different $\beta$ 's (with a different number of null $\beta _i$ 's) involves solving several times equation 7 for different values of

. Equation 7 is a QP (Quadratic Programming) problem and is not trivial to solve at all. The Lasso algorithm, as it is described here, is thus very demanding on computer resources and is nearly of no practical interest. Hopefully, there exist another, more indirect, way to compute at a very little cost all the solutions $\beta$ of the Lasso Algorithm for all the values of

. A modification of the LARS algorithm computes all the Lasso solution in approximatively the time needed to compute and solve the simple ``normal equations'' (equation 4) (the algorithmic complexity of the two procedures is the same). LARS stands for ``Least Angle Regression laSso''. We will not describe here the LARS algorithm but we will give an overview of some of its nicest properties.

The LARS algorithm is a refinement of the FOS (Fast Orthogonal Search) algorithm. We will start by describing the FOS algorithm. The FOS (and the LARS algorithm) is a forward stepwise algorithm. i.e. FOS is an iterative algorithm that, at each iteration, includes inside its model a new variable. The set of variables inside the model at iteration

is the ``active set'' $\mbox{$\cal A$}$

. The set of variables outside the model is the ``inactive set'' $\mbox{$\cal I$}$

. $X_{\cal A}$ is a subset of column of

containing only the active variables at iteration

. $X_{\cal I}$ is a subset of column of

containing only the inactive variables at iteration

.

Let's first illustrate the equation 2: $X \beta =y$ . Equation 2 can be rewritten:

$\displaystyle X \beta = \beta_1 \left( \begin{array}{c} \vdots \\ X_1 \\ \v... ...\right) = \left( \begin{array}{c} \vdots \\ y \\ \vdots \end{array} \right)$

(8)

If $X \in \Re^{m \times n}$ , then $X_i \in \Re^m$ is a column of X.

can also be seen as a vector in a

dimensional space and $\beta \in \Re^n$ is the best linear combination of the

vectors to reach the target

. The ``vector view'' of equation 8 is illustrated in figure 6 where

$\begin{displaymath}X \in \Re^{3 \times 2}= \left( \begin{array}{cc} 1 & 1 \\ 1... ...; y=\left( \begin{array}{c} 3 \\ 2 \\ 3 \end{array} \right) \end{displaymath}$

**Figure 6:** Graphical interpretation of $X \beta =y$
$\begin{figure} \centering\epsfig{figure=vectors1.eps, width=7.5cm, height=4.5cm} \end{figure}$

The FOS algorithm does the following:

Let's define the residual prediction error at step as . We have . Normalize the variables and the target to have $\Vert X_i\Vert _2=1$ and $\Vert y\Vert _2=1$ . Set , the iteration counter.

Figure 7: Illustration of the FOS algorithm

$\begin{figure} \centering\epsfig{figure=vectors2.eps, width=15.5cm, height=4.5cm} \end{figure}$
This step is illustrated in figure 7, Part 1. We will add inside $\mbox{$\cal A$}$ the variable that is ``the most in the direction of ''. i.e. the ``angle'' $\alpha_j$ between and is smaller than all the other angles $\alpha_i$ ( $i\neq j$ ) between the 's and . Using the well-know relation:

$\displaystyle cos(\alpha_i)=\frac{<X_i,E_k>}{\Vert X_i\Vert\; \Vert E_k\Vert}$ (9)

Since the variables have been normalized, we obtain:

$\displaystyle cos(\alpha_i)= C^{te} \; X_i^t E_k \;\;\;\;\;\;($ where $\displaystyle C^{te}=\frac{1}{\Vert E_k\Vert})$ (10)

Note that is simply the correlation between variable and the Residual Error .
Thus the algorithm is: ``compute $a \in \Re^n := X_{\cal I}^t E_k$ . Find such that

$\displaystyle a_j=\max_i \{ a_i \}$ (11)

Add the variable to the active set $\mbox{$\cal A$}$ . Remove variable from the inactive set $\mbox{$\cal I$}$ ''.
This step is illustrated in figure 7, Part 2. We now compute $\beta_j$ , the length of the ``step'' in the direction (inside the example j=1). Let's project (orthogonally) inside the space that is ``spanned'' by the active variables. The projection of is named $E_{proj}$ . We obtain $\beta_j=\frac{\Vert E_{proj}\Vert}{\Vert X_j\Vert}=\Vert E_{proj}\Vert$ since is normalized. If the projection operator fails ( is
inside the space spanned by the active set $\mbox{$\cal A$}$ of variables), remove from the active set $\mbox{$\cal A$}$ and go to step 2.
Deflation: update of the Residual Error: $E_{k+1}:=E_k-E_{proj}$ .
If $\mbox{$\cal I$}$ is empty then stop, otherwise increment and go to step 2.

Let's define

, the number of active variable at the end of the algorithm. If you decide to stop prematurely (for example, when $E_{k-1} - E_k \approx 0$ ), you will have

. The FOS algorithm has several advantages:

It's ultra fast! Furthermore, the computation time is mainly proportional to . If , then you will go a lot faster than when using the standard normal equations.
The algorithm will ``avoid'' naturally the correlated columns. I remind you that a model that is using two strongly correlated columns is likely to be unstable. If column and are strongly correlated, they are both ``pointing'' to the same direction (the angle between them is small: see equation 9 and figure 6) . If you include inside the active set the column , you will obtain, after the deflation, a new error $E_{k+1}$ that will prevent you to choose at a later iteration. One interpretation is the following: ``All the information contained inside the direction has been used. We are now searching for other directions to reach ''. It is possible that, at the very end of the procedure, the algorithm still tries to use column . In this case, the projection operator used at step 3 will fail and the variable will be definitively dropped. Models produced with the FOS algorithm will thus be usually very stable.
The selection of the column that enters the active set is based on the ``residual error'' at step (with ). This selection is thus exploiting the information contained inside the target . The ``target information'' is not used by other algorithms that are based on SVD of or QR factorization of . The FOS algorithm will thus performs a better variable selection than SVD- or QR-based algorithm. Many other algorithms do not use any deflation. Deflation is important because it allows us to search for variables that explains the part of the target that is still un-explained (the residual error ). A good variable selection is useful when .
The memory consumption is very small. To be able to use the FOS algorithm, you only need to be able to store in memory the target , one line of , and a dense triangular matrix of dimension Other algorithms based on QR factorization of requires to have the full matrix inside memory because the memory space used to store is used during the QR factorization to store temporary, intermediate results needed for the factorization. Algorithms based on SVD of are even more memory hungry. This is a major drawback for most advanced applications, especially in Econometrics where we often have $X \in \Re^{10^7 \times 10^4}$ . There exists some ``out of core'' QR and SVD factorizations that are able to work even when does not fit into memory. However ``out of core'' algorithms are currently extremely slow and unstable.

There is however one major drawback to the FOS algorithm. It is illustrated in figure 8 where we have:

$\displaystyle X \in \Re^{3 \times 3}= \left( \begin{array}{ccc} 1 & 0 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 1 \end{array} \right)$ and $\displaystyle \;\; \beta=\left( \begin{array}{c} 1 \\ 1 \\ 0 \end{array} \right) \;\;$ and $\displaystyle \;\; y=\left( \begin{array}{c} 1 \\ 1 \\ 0 \end{array} \right)$

(12)

**Figure:** FOS algorithm fails and gives $\beta=(0 \;\; 0 \;\; 1)^t$
$\begin{figure} \centering\epsfig{figure=fos1.eps, width=15.5cm, height=5.5cm} \end{figure}$

At the first iteration, the FOS algorithm selects

because $X_3^t \;y=2$ is greater than $X_1^t \;y=X_2^t \;y=1$ (and

. Thereafter, the FOS algorithm computes the new Residual Error $E_1=(0 \;\;0 \;\;-1)^t$ and get stuck. Indeed, if we add to the active set either

or

, the L2-norm of the Residual Error do not decrease: $\Vert E_2\Vert _2=\Vert E_1\Vert _2$ . There is no way to add a new variable to decrease the Residual Error and thus the algorithm stops. There are mainly two conceptual reasons why the FOS algorithm fails:

The FOS algorithm is a forward stepwise algorithm and, as almost all forward stepwise algorithms, it is not able to detect multivariate concepts. In the example illustrated in figure 8, the target is defined precisely to be =+. You are able to predict the target accurately if your model is $\beta=(1 \;\; 1 \;\; 0)^t$ . The FOS algorithm gives the wrong model $\beta=(0 \;\; 0 \;\; 1)^t$ . In this example, the target is a multivariate concept hidden inside two variables: and . A forward stepwise algorithm that is able to detect a multivariate concept is very difficult to develop. Inside FOS, the selection of the variable that enters the active set is only based on the univariate relation of with the target . On the contrary, Backward stepwise algorithms have no problem to detect multi-variable concepts and should thus be preferred to a Simple,Classical Forward Stepwise Algorithm. However, usually, Backward stepwise algorithms are very time consuming because they need first to compute a full model and thereafter to select carefully which variables they will drop.

The LARS algorithm with Lasso modification is a forward stepwise algorithm that produces all the solutions of the Lasso algorithm in a computing time proportional to ( is the number of active variable at the end of the algorithm). The Lasso algorithm is part of the family of the backward stepwise algorithm (It starts with a full model and remove some variables when decreases). Thus, the Lasso algorithm is able to discover multi-variable concepts. And thus, the LARS algorithm with Lasso modification is also able to discover multi-variable concepts (since LARS with Lasso gives the same models than the Lasso algorithm). From my knowledge, the LARS algorithm with Lasso modification is the only forward stepwise algorithm able to discover multi-variable concepts. Since, it's a forward stepwise algorithm, it's also very fast.

On the other side, conjunctive concepts are easy to detect. Let's assume that you have the following:
$\begin{displaymath}X \in \Re^{3 \times 2}= \left( \begin{array}{ccc} 1 & 0 \\ ... ...; y=\left( \begin{array}{c} 0 \\ 0 \\ 1 \end{array} \right)\end{displaymath}$

If you create a column ( $X_{3,i}=X_{1,i} X_{2,i}, \;\;i=1,\ldots,m$ ), you obtain:
$\displaystyle X3=\left( \begin{array}{c} 1 \times 0 = 0\\ 0 \times 1 = 0 \\ 1 \times 1 = 1 \end{array} \right) \;\;$ and $\displaystyle \;\; \beta=\left( \begin{array}{c} 0 \\ 0 \\ 1 \end{array} \right)$

i.e. You can represent conjunctive concepts with the product of variables. You will thus obtain second order models. Computing second order models usually involves extending the original matrix with columns that are the product of some of the columns of the original . This is usually a very memory expansive operation. Indeed if is the number of column of , then the extended $X_{extended}$ matrix has $\frac{1}{2}n(n+1)$ columns. The LARS library has been built to allow you to easily create ``on the fly'' a regressor (= one line of $X_{extended}$ ). i.e. you only need to be able to store in memory one line of $X_{extended}$ , not the whole $X_{extended}$ matrix. This allows you to compute at a low memory consumption second, third, fourth,... order models.
The FOS algorithm is ``too greedy'': when the FOS algorithm finds a good direction, it is exploiting this direction ``to the maximum'' (the length of the steps is ). Thereafter there is no ``room left'' for further improvements, adding other variables. The FOS has been ``trapped'' in a so-called ``local optimum''. The LARS algorithm is not so greedy. Its steps lengthes are shorter (see step 3. of the FOS algorithm about ``steps''). This means that the different models that are computed at each steps of the LARS algorithm are not 100% efficient. With the same set of active variables, one can compute a better model. The LARS library downloadable on this page allows you to obtain both models:
- the most efficient one.
- the one used inside the LARS algorithm to search for the next variables to add inside the active set.

To resume, the LARS algorithm has all the advantages of the FOS algorithm:

ultra fast.
correlated columns are avoided to obtain stable models.
Selection of the active variables based on the target (Deflation).
small memory consumption.

Furthermore, the LARS algorithm with Lasso modification is able to discover multivariate concepts easily. For example, the LARS algorithm with LASSO modification finds $\beta=(1 \;\; 1 \;\; 0)^t$ for the example illustrated in figure 8 .

There exists one last question to answer when using the LARS algorithm: When to stop? How to choose

? The question is not simple. Several strategies are available inside the LARS library downloadable on this page:

Let's define the maximum correlation of the inactive columns with , the residual error at iteration (see the definition of in equation 11).

if $\displaystyle ((a(k-10)-a(k)) < s_a \; a(0)$ , then stop (13)

where is a user-defined parameter (usually ) and is the iteration index.
if $\displaystyle ((E_{k-10}^2-E_k^2) < s_E \; E_0^2$ , then stop (14)

where is the residual error at iteration and is a user-defined parameter (usually ).
Let's define where is the number of active variable at iteration .

if $\displaystyle (C_p(k)> max\{ C_p(k-s_c), C_p(k-s_c+1),\ldots, C_p(k-1)\}$ , then stop (15)

where is a user-defined parameter (usually ) and is the iteration index.
Let's define as the n-Fold-Cross-Validation-Error of a model $\beta(k)$ that has been regularized using an optimized ridge parameter $\lambda$ .

if $\displaystyle (minCVError(k)>minCVError(k-1))$ , then stop (16)

where is the iteration index (WARNING: very slow! It involves a full derivative-free optimization algorithm (Brent) at each iteration of the LARS algorithm to find the optimum $\lambda$ ).

My advice is :

to choose a very stringent termination test (like the first one with ) to be sure to have all the informative variables. has a nice and easy interpretation: it's the percentage of
the variance inside that will not be explained.
to perform a backward stepwise, using as initial model, the model given by the LARS algorithm. The LARS library has been designed so that such task is very easy to code.

My experience is that stopping criterions based on MDL (maximum description length), Akkaike index or statistics are not reliable.

Methodology

To be able to answer the simple question ``How good is my model?'' you need a good methodology. Let's assume that we have a set of data

and

. You will cut this data in three parts: a creation set $X_{crea}$ , a validation set $X_{val}$ and a test set $X_{test}$ . The procedure is illustrated in figure 9.

**Figure 9:** Data are cut into three parts
$\begin{figure} \centering\epsfig{figure=sets1.eps, width=8cm, height=6cm} \end{figure}$

The creation set will be used to build a first model. You can build the model $\beta$ using the normal equation ( $X_{crea}^t X_{crea} \beta=X_{crea}^t y_{crea}$ ) or using the LARS algorithm. Thereafter the model $\beta$ needs to be stabilized (stabilization is less important if you have built the model with the LARS algorithm but can still improve a little your model). There are mainly two simple ways to stabilize a model: Remove un-informative variables (this is called ''backward stepwise'') and Ridge Regression.

To remind you, the Ridge Regression is based on equation 5:

$\displaystyle (X_{crea}^t X_{crea} + \lambda I) \; \beta(\lambda) \; = X_{crea}^t y_{crea}$

(17)

where $\lambda$ is the regularization parameter that needs to be adjusted, $\beta(\lambda)$ is the model (depending on the value of $\lambda$ ) and $X_{crea}$ is the creation set. If $\lambda$ is high, the model will be strongly stabilized. The optimal value of the $\lambda$ parameter is $\lambda^*$ . $\lambda^*$ is the solution to the optimization problem:

$\displaystyle \min_{\lambda} \left( (X_{val} \;\; \beta(\lambda) - y_{val} )^2 \right)$

where $\beta(\lambda)$ is computed from equation 17. We are searching for the value of $\lambda$ that gives the best model when applied on the validation set. We are improving the ``generalization ability'' of our model. Since we never used any data from the validation set to build our model, the performance of the model $\beta(\lambda)$ is entirely due to its ``generalization ability''.

To stabilize your model, you can also completely remove some variables using a ``backward stepwise'' algorithm. Here is a simple ``backward stepwise'' algorithm:

Let's define $X_{crea}^{(-i)}$ , the creation set without the column . Let's define $sqE_{ref}= (X_{val} \beta - y_{val})^2$ where $\beta$ is the solution to the standard normal equation: $X_{crea}^t X_{crea} \; \beta\; = X_{crea}^t y_{crea}$ . $sqE_{ref}$ is the performance on the validation set of a model $\beta$ using all the columns of the creation set $X_{crea}$ . Set .
We will try to remove the variable . Compute $\beta^{(-i)}$ , the solution to $X_{crea}^{(-i)t} X_{crea}^{(-i)} \; \beta^{(-i)}\; = X_{crea}^{(-i)t} y_{crea}$ . Compute $sqE^{(-i)}= (X_{val} \beta^{(-i)} - y_{val})^2$ : the performance of the model $\beta^{(-i)}$ on the validation set $X_{val}$ when the model is built without taking into account the $i^{th}$ column of the creation set $X_{crea}$ . If $sqE^{(-i)} \approx sqE_{ref}$ , then the column did not contain any information and can be dropped: Set $X_{crea}:=X_{crea}^{(-i)}$ .
Choose another and go back to step 2.

Let's define $\mbox{$\cal A$}$ , the optimal set of column after the backward stepwise.

Once we have found an optimal ridge parameter $\lambda^*$ and/or the optimal set $\mbox{$\cal A$}$ of informative columns (backward stepwise), we can compute the final model $\beta^*$ using:

$\displaystyle (X_{train}^t X_{train} + \lambda^* I) \; \beta^* \; = X_{train}^t y_{train}$

(18)

where $X_{train}$ contains only the columns in $\mbox{$\cal A$}$ .

WARNING: If you try to use both regularization techniques at the same time, they will ``enter in competition''. You should then be very careful.

To get an idea of the quality of your final model $\beta^*$ , you can apply it on the test set:

Squared Error on test set= $\displaystyle ( X_{test} \; \beta^*- y_{test} )^2$

It's crucial to never use the data in the Test Set $X_{test}$ when building the model because otherwise the results will always be too optimistic. Indeed, if we have used a part of the test set to build the model, we are ``cheating'': we are not truly estimating the ``generalization ability'' of our model.

The methodology just described here is a correct approach to a modelization problem. However, this approach does not use a very large validation and test set. Since the validation set is small, the $\lambda$ parameter and the set $\mbox{$\cal A$}$ of active columns will be very approximative. Since the test set is small, the estimation of the quality of the model is also very approximative. To overcome these difficulties, we can use a n-fold-cross-validation technique. If we want a more accurate estimation of the quality of the model, we can build 5 test sets instead of one: see figure 10. It's named a 5-fold-cross-validation. It also means that we have 5 training sets and thus 5 different models. The main point here is that none of these 5 models have seen the data that are inside their corresponding test set.

**Figure 10:** 5 test sets (and 5 training sets))
$\begin{figure} \centering\epsfig{figure=sets2.eps, width=16cm, height=7cm} \end{figure}$

An estimation of the overall quality of our modelization is:

	$\begin{displaymath}\begin{array}{r} \text{overall quality of the modelization=}... ...rray} \sum_{i=1}^5 ( X_{test(i)} \beta_{(i)} - Y_{test(i)} )^2\end{displaymath}$
	where $\beta_{(i)}$ is the solution to $\displaystyle X_{train(i)}^t X_{train(i)} \; \beta_{(i)} \; = X_{train(i)}^t y_{train(i)}$	(19)

where $X_{train(i)}$ and $X_{test(i)}$ are defined on figure 10.

When we perform a n-fold-cross-validation, we must compute

different models $\beta_{(i)}, \;\; i=1,\ldots,n$ using equation 19. The only time-consuming part inside equation 19 is the computation of $X_{train(i)}^t X_{train(i)}$ . Here is a small ``trick'' to be able to compute $X_{train(i)}^t X_{train(i)}, i=1,\ldots,n$ in a very efficient way: If we define

as in figure 10, we can compute $C_i=B_i^t B_i, \;\; C_i \in \Re^{n \times n}$ . Then, we have:

$\begin{displaymath}\begin{array}{rlclllll} X_{train(1)}^t & X_{train(1)} & = & ... ...train} & = &C_1 +& C_2 +& C_3 + & C_4 + & C_5 \\ \end{array}\end{displaymath}$

(20)

Thanks to equations 20, we can compute $X_{train(i)}^t X_{train(i)}, i=1,\ldots,n$ very easily: only a few matrix additions are needed. The same property of addition exist for the right hand side of equation 19: $X_{train(i)}^t y_{train(i)}$ . Since solving equation 19 is instantaneous on modern CPU's (only the computation of $X_{train}^t X_{train}$ is time-consuming), we can easily obtain the

different models $\beta_{(i)}$ and the n-fold-cross-validation is ``given nearly for free''.

The same technique can be applied to increase the size of the validation set: see figure 11.

**Figure 11:** 4 validation sets (and 4 creation sets))
$\begin{figure} \centering\epsfig{figure=sets3.eps, width=12cm, height=7cm} \end{figure}$

When building a model, the first step is the normalization of all the variables (the columns

). The normalization of a variable

consist of two steps:

We subtract from the column its mean
We divide the column by its standard deviation

Thereafter the

are computed: see equation 20 and figure 10. Based on the

's we can compute easily many different models $\beta_{(i)}$ . Let's assume that we want a model that is completely ignoring the block

: the block

will be used as test set. Referring to figure 10, the test set is $X_{test(5)}$ and the training set is $X_{train(5)}$ . The columns in block

have been normalized using

and

. The two parameters

and

are partially based on information contained inside the block

. We are breaking the rule that says:

``never use any information from the test set to build a model''

(21)

The columns in

(part of the creation set) are normalized using information from

(the test set). This small ``infraction'' of the rules is usually negligible and allows us to reduce considerably the computation time. However, for some rare variables, the rule must be respected strictly. These variables are named ``dynamic variables''. The LARS library does a strong distinction between ``dynamic variables'' and ``static variables''. The LARS library handles both types of variables. However only the ''static variables'' are reviewed inside this documentation.

The concept of ``dynamic variables'' is linked with the concept of ``variable recoding''. The normalization of the variables is a (very basic) recoding. If there is no recoding, there is no ``dynamic variables''. The ``normalization of a variable'' is a recoding that never produces any ``dynamic variables'' (especially if all the lines of

have been correctly scrambled). ``Dynamic variables'' appears mostly when using advanced recodings such as, for example, in the paper ``1-Dimensional Splines as Building Blocks for Improving Accuracy of Risk Outcomes Models'' by David S. Vogel and Morgan C. Wang. The LARS library is handling advanced recodings through the use of ``dynamic variables''.

Inside the LARS library, the $\lambda$ parameter of the Ridge Regression is computed using a n-fold-cross-validation on the validation sets. Equations 20 are used to achieve high efficiency. The $C_i \in \Re^{n_a \times n_a}$ matrices are built incrementally each time a new variable enters the active set $\mbox{$\cal A$}$ . The equations 19 are solved using an advanced Cholesky factorization that is using an advanced pivoting strategy and an advanced perturbation strategy to obtain a very high accuracy model even when the matrices $(X_{\cdot(i)}^t X_{\cdot(i)})$ are close to singularity or badly conditioned.

The LARS library has been written to use the less memory possible. For example, all the

matrices are symmetric, thus we only need to store half of these matrices in memory. ``In place'' factorization of matrices have been used when possible.

It can happen inside the LARS forward stepwise algorithm that several variables should enter the active set $\mbox{$\cal A$}$

at the same iteration

(there are several values of

that are possible inside equation 11 ). It can also happens that several variables should leave the active set $\mbox{$\cal A$}$

at the same iteration

. Most algorithms do not treat correctly these special cases. The LARS library is handling correctly these cases.

How to use the library?

The LARS library inside Matlab

The usage is:

[betaBest,betas,errors,cHats]=matlabLARS(X,y,stoppingCriteriaType,stoppingCriteriaValue,lassoModif, verbosity);

matlabLARS is an iterative algorithm that compute different $\beta$ 's $\in \Re^n$ such that $(X \beta- y)^2$ is minimum (classical least-square). The X matrix and the y vector do not need to be normalized.
At the first iteration only one $\beta _i$ is non-null. At each iteration

, we will:

set one more $\beta _i$ to a non-null value.
compute Error where Error $=(X \beta - y)$
compute cHat= (see equation 13 about )

The different $\beta$ 's, computed at each iterations, are stored inside the matrix 'betas'. The best model among all the models according to the selected stopping criteria is stored in 'betaBest'. In a similar way, at each iteration, the squared error Error

and the maximum correlation values cHat, are stored respectively inside the matrices 'errors' and ' cHats'.

If the lasso modification is used (set lassoModif=1), then some $\beta _i$ could be reset to a null value at some rare iterations.

Let's define cHat(k) the value of cHat at iteration

.
Let's define sqEr(k) the value of Error at iteration

.
Let's define n(k) the number of non-null $\beta _i$ at iteration

.
The variable 'stoppingCriteriaType' is a character that can be:

'R': if ((cHat(k-10)-cHat(k))< stoppingCriteriaValue cHat(0)), then stop
'S': if ((sqEr(k-10)-sqEr(k))< stoppingCriteriaValue sqEr(0)), then stop
'C': Cp test: Let's define
- Cp(k) = sqEr(k)-nCol+2*n(k)
- p =stoppingCriteriaValue
if (Cp(k)> max Cp(k-p), Cp(k-p+1),..., Cp(k-1)), then stop
'M': Let's define minCVError(k) as the 6-Fold-Cross-Validation-Error of a model beta(k) that has been regularized using an optimized ridge parameter $\lambda$ .
if (minCVError(k)>minCVError(k-1)), then stop
'E': if (stoppingCriteriaValue=-1) include all columns, otherwise:
if (number of non-null $\beta _i$ =stoppingCriteriaValue), then stop
The variable 'verbosity' is optional. It's a number between 0 and 2. (2=maximum verbosity).

As previously mentioned, if you select a specific set of active variable, you can compute at least two different models $\beta$ :

the most efficient one (stored inside betaBest)
the one used inside the LARS algorithm to search for the next variables to add inside the active set (stored inside betas).

betaBest $_{n+1}$ is the constant term of the model. All the models (in betas and betaBest) are computed to be applied on un-normalized data.

The Matlab interface to the LARS library is very poor compared to the C++ interface. The source code of the matlab interface (in C++) is inside the zip-file downloadable at the top of this webpage. If you want to customize the Matlab version of the LARS library, it's very easy to edit/modify this small C++ code. The C++ interface is documented below.

The LARS library as a stand-alone executable

Inside the zip-file downloadable inside this webpage, you will find the source code of a small C++ executable that solves the `diabetes' example. The `diabetes' data are stored inside an external .csv file. The last column of the .csv file is the target. A very fast .csv-file-reader is provided.

It's very easy to edit/modify the small C++ code to customize it to your needs. The C++ interface is documented below.

The LARS library in C++

The LARS library is available as a .dll file. To use this .dll file in your own project you must:

In the ``debug'' version of you program, add inside the ''Linker/Input/Additional Dependencies'' tab: ``larsDLLDebug.lib''
Inside the directory containing the executable of the ''debug'' version of you program, add the file ''larsDLLDebug.dll''
In the ``release'' version of you program, add inside the ''Linker/Input/Additional Dependencies'' tab: ``larsDLL.lib''
Inside the directory containing the executable of the ''release'' version of you program, add the file ``larsDLL.dll''
Add inside your project the files ``lars.h'' and ''miniSSEL1BLAS.h''

In order to work, the LARS library must access the

matrix and the

vector. The normalization of

and

is done automatically inside the LARS library, if needed. You can provide the data (

and

) in two different ways:

Line by Line: You must create a child of the LARS class. Inside the child class, you must re-define some of the following functions:
- larsFLOAT *getLineX(int i, larsFLOAT *b);
  This function return the $i^{th}$ line of . The required line can be stored inside the b array or inside an other array. A pointer to the array containing the line must be returned.
- larsFLOAT getLineY(int i);
  This function return the $i^{th}$ component of : . Alternatively, you can use the third parameter of the getModel function to specify a target.
- void setSelectedColumns(int n,int *s);
  The re-definition of this function is optional. The default implementation is:
```
void LARS::setSelectedColumns(int n,int *s) { selectedC=s; nSelected=n; }
```
  The setSelectedColumns function announces that subsequent calls to the function getLineXRequired will return only a part of the column of . The indexes of the columns that should be returned are given inside the array s of length n.
- larsFLOAT *getLineXRequired(int i, larsFLOAT *b);
  The re-definition of this function is optional. However, it's strongly recommended that you re-define this function for performance reasons. The default implementation is:
```
larsFLOAT *LARS::getLineXRequired(int i, larsFLOAT *b)
{
    int j,ns=nSelected, *s=selectedC;
    larsFLOAT *f=getLineX(i,b);
    for (j=0; j<ns; j++) b[j]=f[s[j]];
    return b;
}
```
  This function returns only the columns of line i that have been selected using the function setSelectedColumns.
Column by Column: You must create a child of the LARS class. Inside the child class, you must re-define the following functions:
- larsFLOAT *getColX(int i, larsFLOAT *b);
  This function return the $i^{th}$ column of . The required column can be stored inside the b array or inside an other array. A pointer to the array containing the line must be returned.
- larsFLOAT *getColY(larsFLOAT *b);
  This function return the target . The target can be stored inside the b array or inside an other array. A pointer to the array containing must be returned. Alternatively, you can use the third parameter of the getModel function to specify a target.

Let's now assume that you have created a child of the base LARS class named LARSChild. You can now configure the different parameters of the LARS algorithm. Typically, you define these parameters inside the constructor of LARSChild. The parameters are:

nCol: the number of column of .
nRow: the number of rows of .
orientation: if orientation='L', then we access the matrix line by line. Otherwise, we access the matrix column by column.
lasso: if lasso=1, then the lasso modification of the LARS algorithm is used.
fullModelInTwoSteps: if fullModelInTwoSteps=1, then no forward stepwise is performed: we are directly building a model using all the column of . When the number of column nCol is small this is faster than the full LARS algorithm.
ridge: if ridge=1, then an optimization of the $\lambda$ parameter of the ridge regression is performed at the end of the build of the model (this requires ).
memoryAligned: all the linear algebra routines are SSE2 optimized. This means that the functions getCol* and getLine* should return a pointer to an array that is 16-byte-memory-aligned. Such an array ca be obtained using the mallocAligned function of the miniSSEL1BLAS library. If you are not able to provide aligned-pointers, you can set memoryAligned=0, but the computing time will be more or less doubled.
nMaxVarInModel: maximum number of variables inside the model when performing a forward stepwise (this is a stopping criteria).
nFold: the number of block (as defined in figure 10) for the n-fold-cross-validation.
folderBounds: an array of size of integer that contains the line-index of the first line of each (as defined in figure 10). We should also have folderBounds[nFold]=nRow. If you don't initialize this parameter, then folderBounds will be initialized for you automatically so that all the blocks have the same size.
stoppingCriteria: define which kind of stopping criteria will be used inside the forward stepwise:
- stoppingCriteria='R': stop based on the maximum correlation