Least Angle Regression in Tangent Space and LASSO for Generalized Linear Model¹

In this paper, we consider the generalized linear regression, which is an estimation problem of the GLM [12]. In the generalized linear regression, the expectation $\mu$ of a response $y$ is represented by a linear combination of explanatory variables $x_1,x_2,\dots ,x_d$ as

where $h:R\to R$ is called a link function, $n$ is the sample size, $d$ is the number of the explanatory variables, and $\theta =\left({\theta }^i\right)$ is the parameter to be estimated. Let $X=\left(x_i^a\right)$ be the design matrix, which is an $(n\times d)$-matrix. Let $y=\left(y^a\right)$ and $\mu =\left({\mu }^a\right)$ be the response vector and its expectation, respectively, which are column vectors of length $n$. In general, the link function is a function of ${\mu }^a$ and is not determined uniquely. However, in the paper, we only use the canonical link function, which results in useful properties of the GLM.

As a special case, the normal linear regression uses the link function $h\left(y\right)=y$ and a quadratic function as the potential function. Another example is the logistic regression, where the link function is $h\left(y\right)=y/(1-y)$ and the potential function is $\psi \left(\theta \right)={\sum }_{a=1}^n\log \{1+\exp ({\sum }_{i=1}^d{x}_i^a{\theta }^i\left)\right\}$.

Through the paper, we assume that the design matrix $X$ is normalized, that is, each column vector has the mean zero and the $l_2$-norm one: ${\sum }_{a=1}^n{x}_i^a=0$ and ${\sum }_{a=1}^n(x_i^a{)}^2=1$ for $i=1,2,\dots ,d$. Furthermore, we assume that column vectors of $X$ are linearly independent.

We briefly describe the LARS algorithm. In subsection 3.2, we use the LARS algorithm for proposing an estimation method. The detail and more discussions on LARS can be found in, for example, [6] and [8].

LARS was proposed as an algorithm for parameter estimation and variable selection in the normal linear regression. In the LARS algorithm, the estimator moves from the origin to the maximum likelihood estimate (MLE) ${\hat{\theta }}_{MLE}$ of the full model. The full model means the linear model including all the explanatory variables. The MLE ${\hat{\theta }}_{MLE}$ is determined by the design matrix $X$ and the response $y$. The detailed algorithm of LARS is showed in Algorithm 1, where ${\hat{\theta }}_{\left(k\right)}$ is $k$-th estimate the algorithm outputs. After iterations, LARS outputs a sequence of the estimates ${\hat{\theta }}_{\left(0\right)},{\hat{\theta }}_{\left(1\right)},\dots ,{\hat{\theta }}_{\left(d\right)}$.

The idea of the LARS algorithm is showed by Figure 1. Figures 1 and 1 indicate the estimator’s move and the residual’s move, respectively, in the parameter space $R^d$ when $d=2$. The estimator i) selects an element of the parameter which makes a least angle between ${\hat{\theta }}_{MLE}-{\hat{\theta }}_{\left(k\right)}$ and ${\theta }^i$-axis, and ii) uses it as a trajectory in the form of the bisector of an angle. The LARS algorithm is described in terms of Euclidean geometry and can be computed efficiently. Furthermore, $X^{\top }X$ plays an important role in the LARS algorithm, which is one of our motivations for considering the tangent space of a statistical model.

LASSO is an optimization problem for parameter estimation and variable selection in the normal linear regression. LASSO solves the minimization problem

where $\lambda \ge 0$ is a tuning parameter. The path of the LASSO estimator when $\lambda$ varies can be made by the LARS algorithm with a minor modification.

We introduce some tools from information geometry, including model manifold, tangent space, and exponential map (Figure 2). This is a brief introduction. For details, see [2, 3, 1, 5].

In the generalized linear regression, we need to select one distribution from the GLM. A model manifold is a manifold consisting of probability distributions of interest. That is, the model manifold is $M=\left\{p\right(\cdot |\theta )|\theta \in R^d\}$, where $p(\cdot |\theta )$ indicates the probability distribution with the regression coefficient $\theta$. The parameter $\theta$ works as a coordinate system in $M$.

The tangent space $T_{p}M$ at a point $p\in M$ is a linear space consisting of directional derivatives, that is, $T_{p}M=\{v={\sum }_{i=1}^d{v}^i{\partial }_i|v^i\in R\}$, where ${\partial }_i:=\partial /\partial {\theta }^i$. We consider the tangent space $T_{p(\cdot |0)}M$ at $p(\cdot |0)$. For simplicity, we call $p(\cdot |0)$ and $T_{p(\cdot |0)}M$, the origin and the tangent space $T_{0}M$ at the origin, respectively.

Any pair of two vectors in $T_{0}M$ has its inner product. The inner product is determined by the Fisher information matrix $G=G\left(0\right)=\left(g_{ij}\right(0\left)\right)$:

where $l\left(\theta \right)=\log p\left(y|\theta \right)$ is the log-likelihood. Using the Fisher metric $G$, the inner product of $v_1={\sum }_{i=1}^d{v}_1^i{\partial }_i$ and $v_2={\sum }_{i=1}^d{v}_2^j{\partial }_j$ is given by

In the generalized linear regression, the Fisher metric $G$ at $T_{0}M$ is proportional to the correlation matrix $X^{\top }X$ of the explanatory variables. That is, $G=c{X}^{\top }X$ for some $c>0$. For details, see Lemma 2 in subsection A.1.

A point in the tangent space $T_{0}M$ can be identified with a point in $M$ via an exponential map. We introduce the e-exponential map ${Exp}_0:T_{0}M\to M$ defined as follows. For $v=\sum v^i{\partial }_i\in T_{0}M$, let ${Exp}_0\left(v\right)=p(\cdot |v)\in M$ with $v=\left(v^i\right)$. Our problem in this paper is estimation for the GLM and the parameter is a regression coefficient vector $\theta \in R^d$. Therefore, we can avoid technical difficulties of an exponential map. The map ${Exp}_0$ is a bijection from $T_{0}M$ to $M$. For details, see subsection A.3.

For readers familiar with information geometry, we make an additional remark. The model manifold $M$ of the GLM is e-flat and the regression coefficient $\theta$ is an e-affine coordinate system of $M$. $\left\{{\partial }_i\right\}$ is the natural basis of $T_{0}M$ with respect to the coordinate system $\theta$. Each coordinate axis of ${\theta }^i$ in $M$ corresponds to ${\partial }_i$-axis in $T_{0}M$ via the e-exponential map.

In the following, we also use another representation of $T_{0}M$. This representation is useful for our purpose: $T_{0}M=\{X\theta |\theta \in R^d\}$. In our notation, $X\theta$ also indicates $\sum {\theta }^i{\partial }_i$ in the tangent space $T_{0}M$, not only a point $p(\cdot |\theta )\in M$. However, we believe that it is not confusing because a vector in the tangent space and a point in $M$ are identified through the exponential map.

The main idea of the proposed method is to run LARS in the tangent space $T_{0}M$ at the origin. First, we correspond the model manifold to the tangent space $T_{0}M$ by the e-exponential map. After this mapping, our computation is done by the original LARS algorithm. However, we do not use the response $y$ directly. We introduce a virtual response $\hat{y}$. The LARS algorithm outputs a sequence of parameter estimates, whose length is the same as the dimension of the parameter. Finally, the estimates are mapped to the model manifold.

Before running the original LARS algorithm, we introduce the virtual response $\hat{y}$. The virtual response $\hat{y}$ is defined using the design matrix $X$ and the MLE ${\hat{\theta }}_{MLE}$ of the full model: $\hat{y}=X{\hat{\theta }}_{MLE}$. Note that LARS uses only correlation coefficients between the response $y$ and the explanatory variables $X$ in the form of $y^{\top }X\theta$, which is identical with ${\hat{\theta }}_{MLE}^{\top }{X}^{\top }X\theta$. Therefore, introducing the appropriate representation of the response $\hat{y}=X{\hat{\theta }}_{MLE}$, we need only $X^{\top }X$ as ${\hat{y}}^{\top }X\theta ={\hat{\theta }}_{MLE}^{\top }{X}^{\top }X\theta$.

In the estimation step of the proposed method, we run the original LARS algorithm in the tangent space $T_{0}M$ as if the response is $\hat{y}$. LARS outputs a sequence $\{{\hat{\theta }}_{\left(0\right)},{\hat{\theta }}_{\left(1\right)},\dots ,{\hat{\theta }}_{\left(d\right)}\}$ of the model parameter. As is shown in Figure 1, the LARS estimator $\hat{\theta }$ can be regarded as moving from the origin to the MLE ${\hat{\theta }}_{MLE}$ of the full model. At the same time, however, the residual $r\left(\hat{\theta }\right):={\hat{\theta }}_{MLE}-\hat{\theta }$ of the estimator $\hat{\theta }$ is moving from the MLE ${\hat{\theta }}_{MLE}$ to the origin (Figure 1). The latter is useful for our method because it allows us to fix the estimator’s tangent space to the origin. What moves is the residual $r\left(\hat{\theta }\right)$, not the estimator $\hat{\theta }$. Note that Algorithm 1 in subsection 2.2 is actually described from the latter perspective.

The Kullback-Leibler divergence (KL divergence) is a key quantity in information geometry, which is also important in statistics, machine learning and information theory. For the GLM, the KL divergence $D$ is given by

where $\eta =E_{\theta }\left[X^{\top }y\right]$ is the expectation parameter of the exponential family. The KL divergence is approximated up to second order as

where $d\theta =\left(d{\theta }^i\right)$.

In generalized linear regression, the Fisher metric $G=\left(g_{ij}\right)$ is proportional to the correlation matrix $X^{\top }X$, that is, $G=c{X}^{\top }X$ for some $c>0$. (See Appendix A.1.) The KL divergence is approximately related with the correlation matrix as

We propose two estimation methods. One is LASSO modification of TLARS. The other is an approximation of the $l_1$-regularization for the GLM (1).

We briefly compare TLARS with two existing methods which are extensions of LARS based on information geometry. One is Bisector Regression (BR) by [9] and the other is Differential-Geometric LARS (DGLARS) by [4]. Our concern here is about algorithm itself.

First, the BR algorithm is very different from TLARS. BR takes advantage of the dually flat structure of the GLM and tries to make an equiangular curve using the KL divergence. Furthermore, the BR estimator moves from the MLE of the full model to the origin while, in our method, the residual moves from ${\hat{\theta }}_{MLE}$ to the origin.

DGLARS is also different from TLARS. It uses tangent spaces, where the equiangular vector is considered. However, the DGLARS estimator actually moves from $p(\cdot |0)$ to $p(\cdot |{\hat{\theta }}_{MLE})$ in $M$. Accordingly, the tangent space at the current estimator moves, which makes us treat the tangent spaces at many points in $M$. DGLARS treats the model manifold directly. Therefore, it requires many iterations of approximation computation for the algorithm. Note that, on the other hand, the update of the TLARS estimator is described fully in terms of only the tangent space $T_{0}M$.

We applied the proposed methods and the L1 method to a real data. The data is the South Africa heart disease (SAheart) data included by ElemStatLearn package of R. This data contains nine explanatory variables of 462 samples. The response is a binary variable.

We show the results by the four methods. Figures 3 and 3 are the paths by TLARS and TLASSO1, respectively. In this example, they are the same. Figure 3 is the TLASSO2 path, and Figure 3 is the L1 path. The paths by TLARS, TLASSO1, and TLASSO2 are made by lars() function of R, and that of L1 by glmnet().

As Figure 3 shows, the four paths are very similar. The proposed methods are based only on the tangent space, not on the model manifold itself, while L1 directly takes advantage of the likelihood. These results imply that the approximation of the model does not require deterioration of result for our methods, especially, for TLARS and TLASSO1.

We performed numerical experiments of logistic regression. The topic is three-fold: generalization, parameter estimation, and model selection. The result is shown in Table 1. Bold values are the best and better values.

The procedure of the experiments is as follows. We fixed the number of the parameter $d$, the true value of the parameter ${\theta }_0$, and the sample size $n$. For each of $m$ trials, we made the design matrix $X$ by rnorm() function in R. Furthermore, we made the response $y$ based on $X$ and ${\theta }_0$, that is, elements of $y$ have different Bernoulli distributions. The four methods were applied to $(y,X)$.

For evaluating the generalization error of the four methods, we newly made $m$ observations $\left\{\right(y_1^l,X_1^l),\dots ,(y_m^l,X_m^l\left)\right\}$ in $l$-th trial ($l=1,2,\dots ,m$). We computed the difference between $(y_1^l,\dots ,y_m^l)$ and $m$ predictions by each of the methods. The “Generalization” columns of Table 1 show the average prediction error over $m$ trials. Smaller value is better.

The “Model selection” columns show the fraction of the trials (among $m$ trials) where the methods selected the true model. The “Seq” column indicates the fraction of the trials where each sequence of estimates included the true model. Larger value is better.

In the “Parameter estimation” columns, each value means the average of $\parallel \hat{\theta }-{\theta }_0{\parallel }_2^2$ of the selected estimate $\hat{\theta }$. Smaller value is better.

In Table 1, we report the results of three cases. We used $m=10,000$ for all cases but case C2 where $m=1,000$. In case A, we set $d=10$ and ${\theta }_0=(10,10,10,-10,$ $-10,-10,0,0,0,0{)}^{\top }$. We used $n=100$ for case A1 and $n=1,000$ for A2. In generalization, three methods (TLARS, TLASSO1, and L1) with AIC2 were much better than the other combinations of method and information criterion. In model selection, the four methods with BIC1 were much better regardless of the sample size. In parameter estimation, TLARS and TLASSO1 with AIC1 and BIC2 were better in the small sample setting. However, in the larger sample setting, the four methods with AIC2 were better. These tendencies were observed in other cases not reported here; For example, ${\theta }_0=(10,10,-10,-10,0,0,0,0,0,0{)}^{\top }$.

Case B is the case of $d=10$ and ${\theta }_0=(10,10,0,0,0,$ $0,0,0,0,0{)}^{\top }$ with the relation $x_3=x_2+ϵ$, where $x_2$ and $x_3$ are the second and third columns of the design matrix $X$, respectively, and $ϵ$ is distributed according to a multivariate normal distribution. We set $n=100$ and $n=1,000$ for cases B1 and B2, respectively. In generalization, TLARS and TLASSO1 with AIC1, BIC1, and BIC2 were better than the others in Case B1. Three methods (TLARS, TLASSO1, and L1) with AIC1 and BIC2 were better for the larger sample setting. In Case B, our interest is mainly in generalization because estimation of the true model and the parameter value are not very meaningful. However, the four methods with BIC1 were better in model selection.

In case C, we used $d=50$ and, as ${\theta }_0$, the vector of the length 50 with ten $10$s, ten $-10$s, and thirty $0$s. In generalization and parameter estimation, three methods (TLARS, TLASSO1, and L1) with AIC2 were better than the others regardless of the sample size. In model selection, the four methods with BIC1 were much better than the others.

In summary, the proposed methods worked very well. Of course, the L1 method sometimes performs better than our methods. However, the proposed methods, especially TLARS and TLASSO1, are better than L1 in many situations. Furthermore, TLARS and TLASSO1 output the same results in very many trials.