Least Angle Regression in Tangent Space and LASSO for Generalized Linear Models

TL;DR

This paper introduces a novel approach for sparse estimation in generalized linear models by adapting LARS and LASSO to the tangent space of the model manifold, enabling efficient computation and improved performance.

Contribution

It extends LARS and LASSO algorithms to generalized linear models through tangent space transformation, enhancing efficiency and estimation accuracy.

Findings

Methods perform comparably to $l_1$-regularized maximum likelihood estimation.

Numerical experiments show improved generalization and model selection.

Real-data analysis confirms the effectiveness of the proposed methods.

Abstract

This study proposes sparse estimation methods for the generalized linear models, which run one of least angle regression (LARS) and least absolute shrinkage and selection operator (LASSO) in the tangent space of the manifold of the statistical model. This study approximates the statistical model and subsequently uses exact calculations. LARS was proposed as an efficient algorithm for parameter estimation and variable selection for the normal linear model. The LARS algorithm is described in terms of Euclidean geometry regarding the correlation as the metric of the parameter space. Since the LARS algorithm only works in Euclidean space, we transform a manifold of the statistical model into the tangent space at the origin. In the generalized linear regression, this transformation allows us to run the original LARS algorithm for the generalized linear models. The proposed methods are efficient and perform well. Real-data analysis indicates that the proposed methods output similar results to that of the $l_1$-regularized maximum likelihood estimation for the aforementioned models. Numerical experiments reveal that our methods work well and they may be better than the $l_1$-regularization in generalization, parameter estimation, and model selection.

Figures

Tables

Table 1. Table 1 : The results of the numerical experiments. Generalization: the average prediction error. Model selection: the proportion of the trials where the methods selected the true model. Seq: the proportion of the trials in which each sequence of estimates included the true model. Parameter estimation: the average of the squared error of the selected estimate. Values in bold are the best and better values.

	Method	Generalization  ($\times {10}^{-2}$)				Model selection					Parameter estimation
		AIC1	AIC2	BIC1	BIC2	Seq	AIC1	AIC2	BIC1	BIC2	AIC1	AIC2	BIC1	BIC2
A1	TLARS	10.70	9.80	12.97	10.72	0.7246	0.3969	0.1838	0.4973	0.3672	168.3	178.7	195.5	167.4
	TLASSO1	10.70	9.80	12.97	10.72	0.7247	0.3968	0.1838	0.4974	0.3784	168.3	178.7	195.5	167.4
	TLASSO2	15.73	12.49	18.35	15.04	0.7086	0.4062	0.0662	0.4865	0.2769	249.3	171.4	310.2	232.4
	L1	18.81	9.74	22.60	10.46	0.6897	0.3996	0.0301	0.4824	0.1548	315.7	183.5	404.5	169.1
A2	TLARS	4.04	3.60	5.14	3.96	0.9785	0.4955	0.1252	0.8573	0.4988	58.7	45.6	99.6	56.4
	TLASSO1	4.04	3.58	5.14	3.96	0.9785	0.4955	0.1252	0.8573	0.4988	58.7	45.6	99.6	56.4
	TLASSO2	4.73	3.69	5.91	4.40	0.9787	0.4959	0.0561	0.8575	0.4022	79.6	47.2	126.9	69.0
	L1	8.20	3.59	10.57	4.05	0.9732	0.4968	0.0721	0.8570	0.3810	234.3	45.4	352.9	58.7
B1	TLARS	13.52	13.89	13.42	13.16	0.5500	0.1799	0.1455	0.4293	0.3369	146.7	221.6	104.9	106.5
	TLASSO1	13.33	13.67	13.28	13.00	0.5643	0.1784	0.1508	0.4316	0.3474	144.0	214.2	102.8	102.1
	TLASSO2	14.05	13.94	14.81	14.23	0.5666	0.1820	0.1152	0.4366	0.3257	105.0	133.8	106.2	100.0
	L1	15.98	14.43	19.16	13.59	0.5560	0.1814	0.0785	0.4342	0.2671	131.4	334.1	155.9	101.4
B2	TLARS	4.95	5.20	5.16	4.95	0.5848	0.1926	0.1127	0.5402	0.4157	96.8	140.6	89.9	84.8
	TLASSO1	4.95	5.20	5.16	4.94	0.5852	0.1918	0.1127	0.5400	0.4159	96.8	140.6	89.9	84.8
	TLASSO2	5.00	5.31	5.29	5.05	0.5850	0.1926	0.0978	0.5400	0.4104	95.0	143.7	90.8	85.1
	L1	5.90	5.10	7.62	4.90	0.5793	0.1925	0.0935	0.5367	0.3946	109.1	132.2	158.8	82.0
C1	TLARS	10.18	9.27	14.71	10.97	0.1479	0.0087	0.0009	0.0787	0.0237	373.3	324.8	751.4	435.9
	TLASSO1	10.18	9.27	14.71	10.97	0.1479	0.0087	0.0009	0.0787	0.0237	373.3	324.8	751.4	435.9
	TLASSO2	14.78	11.65	19.09	15.72	0.1399	0.0098	0.0000	0.0742	0.0147	811.9	537.0	1183.8	895.2
	L1	13.48	9.48	19.48	12.05	0.1137	0.0088	0.0000	0.0706	0.0050	690.5	351.6	1215.4	566.4
C2	TLARS	3.98	3.36	6.22	4.17	0.773	0.014	0.000	0.486	0.077	247.4	172.1	608.9	274.2
	TLASSO1	3.98	3.36	6.22	4.17	0.773	0.014	0.000	0.486	0.077	247.4	172.1	608.9	274.2
	TLASSO2	4.45	3.53	6.64	4.57	0.779	0.014	0.000	0.486	0.068	311.6	190.4	687.7	329.0
	L1	4.58	3.40	8.08	4.34	0.736	0.015	0.000	0.486	0.046	330.8	176.0	982.7	297.9

Equations

$$h({\mu }^a)=\sum _{i=1}^d{x}_i^a{\theta }^i,(a=1,2,\dots ,n),$$

$$\mathcal{M}=\{p(\cdot |\mathbf{\theta })|\mathbf{\theta }\in {\mathbb{R}}^d\},$$

$$p(y|\mathbf{\theta })=p(y|X,\mathbf{\theta })=\exp \left\{{\mathbf{y}}^{\top }X\mathbf{\theta }-\psi (\mathbf{\theta })\right\},$$

$${\hat{\mathbf{c}}}_{(k)}:=X^{\top }X{\mathbf{r}}_{(k-1)},{\hat{C}}_{(k)}:=\underset{j}{\max }\{|{\hat{c}}_{(k),i}|\},$$

$$I_{(k)}:=\{i||{\hat{c}}_{(k),i}|={\hat{C}}_{(k)}\}.$$

$$X_{(k)}:=(\dots s_i{x}_i\dots {)}_{i\in I_{(k)}},G_{(k)}:=X_{(k)}^{\top }{X}_{(k)},$$

$$A_{(k)}:=(1_k^{\top }{G}_{(k)}^{-1}{1}_k{)}^{-1/2},{\mathbf{w}}_{(k)}:=A_{(k)}{G}_{(k)}^{-1}{1}_k,$$

$${\mathbf{a}}_{(k)}:=X^{\top }{X}_{(k)}{\mathbf{w}}_{(k)}.$$

$$\{\begin{array}{l}({\hat{\theta }}_{(k)}^i{)}_{i\in I_{(k)}}:=({\hat{\theta }}_{(k-1)}^i{)}_{i\in I_{(k)}}+\gamma \mathrm{diag}(s_i{)}_{i\in I_{(k)}}{\mathbf{w}}_{(k)},\\ ({\hat{\theta }}_{(k)}^i{)}_{i\in I_{(k)}^{\mathrm{c}}}:=0\end{array}$$

$$\hat{\gamma }:=\underset{j\in I_{(k)}^c}{\min }{}^{+}\left\{\frac{{\hat{C}}_{(k)}-{\hat{c}}_{(k),j}}{A_{(k)}-a_{(k),j}},\frac{{\hat{C}}_{(k)}+{\hat{c}}_{(k),j}}{A_{(k)}+a_{(k),j}}\right\}>0,$$

$$\underset{\mathbf{\theta }\in {\mathbb{R}}^d}{\min }\left\{\Vert \mathbf{y}-X\mathbf{\theta }{\Vert }_2^2+\lambda \Vert \mathbf{\theta }{\Vert }_1\right\},$$

$$\underset{\mathbf{\theta }\in {\mathbb{R}}^d}{\min }\left\{-{\mathbf{y}}^{\top }X\mathbf{\theta }+\psi (\mathbf{\theta })+\lambda \Vert \mathbf{\theta }{\Vert }_1\right\}.$$

$$g_{ij}(\mathbf{\theta })=\mathrm{E}\left[{\partial }_{i}l(\mathbf{\theta }){\partial }_{j}l(\mathbf{\theta })\right],$$

$$⟨v_1,v_2⟩=\sum _{i=1}^d\sum _{j=1}^d{v}_1^i{v}_2^j⟨{\partial }_i,{\partial }_j⟩=\sum _{i=1}^d\sum _{j=1}^d{v}_1^i{v}_2^j{g}_{ij}.$$

$$\underset{\mathbf{\theta }\in {\mathbb{R}}^d}{\min }\left\{\Vert X{\hat{\mathbf{\theta }}}_{\mathrm{MLE}}-X\mathbf{\theta }{\Vert }_2^2+\lambda \Vert \mathbf{\theta }{\Vert }_1\right\},$$

$$\underset{\mathbf{\theta }\in {\mathbb{R}}^d}{\min }\left\{\Vert \alpha X\tilde{\mathbf{\theta }}-X\mathbf{\theta }{\Vert }_2^2+\lambda \Vert \mathbf{\theta }{\Vert }_1\right\},$$

$$\log p(y|\mathbf{\theta })\approx -\frac{1}{2\alpha }(\mathbf{\theta }-\alpha \tilde{\mathbf{\theta }}{)}^{\top }{X}^{\top }X(\mathbf{\theta }-\alpha \tilde{\mathbf{\theta }})+\frac{\alpha }{2}{\tilde{\mathbf{\theta }}}^{\top }\tilde{\mathbf{\theta }}-\psi (0).$$

$$\mathrm{AIC}$$

$$\mathrm{BIC}$$

$$\displaystyle\frac{\partial\psi}{\partial\theta^{i}}=\sum_{a=1}^{n}x_{i}^{a}h^{-1}\Big{(}\sum_{j=1}^{d}x_{j}^{a}\theta^{j}\Big{)}.$$

$$\left(\frac{\partial \psi }{\partial {\theta }^1},\frac{\partial \psi }{\partial {\theta }^2},\dots ,\frac{\partial \psi }{\partial {\theta }^d}\right)$$

$$\displaystyle=\Big{(}h^{-1}\Big{(}\sum_{i=1}^{d}x_{i}^{1}\theta^{i}\Big{)},h^{-1}\Big{(}\sum_{i=1}^{d}x_{i}^{2}\theta^{i}\Big{)},$$

$$\displaystyle\hskip 56.9055pt\dots,h^{-1}\Big{(}\sum_{i=1}^{d}x_{i}^{n}\theta^{i}\Big{)}\Big{)}X.$$

$$g_{ij}(\mathbf{\theta })$$

$$\displaystyle=\sum_{a=1}^{n}x_{i}^{a}\partial_{j}\mu^{a}({\bm{\theta}})=\sum_{a=1}^{n}x_{i}^{a}\partial_{j}h^{-1}\Big{(}\sum_{k=1}^{d}x_{k}^{a}\theta^{k}\Big{)}$$

$$\displaystyle=\sum_{a=1}^{n}x_{i}^{a}x_{j}^{a}\tilde{h}\Big{(}\sum_{k=1}^{d}x_{k}^{a}\theta^{k}\Big{)},$$

$$\log p(y|\mathbf{\theta })=-\frac{1}{2\alpha }(\mathbf{\theta }-\alpha \tilde{\mathbf{\theta }}{)}^{\top }{X}^{\top }X(\mathbf{\theta }-\alpha \tilde{\mathbf{\theta }})+\frac{\alpha }{2}{\tilde{\mathbf{\theta }}}^{\top }\tilde{\mathbf{\theta }}-\psi (0)+O(\Vert \mathbf{\theta }{\Vert }^3).$$

$$\psi (\mathbf{\theta })$$

$$=\psi (0)+h^{-1}(0)1^{\top }X\mathbf{\theta }+\frac{1}{2\alpha }{\mathbf{\theta }}^{\top }{X}^{\top }X\mathbf{\theta }+O(\Vert \mathbf{\theta }{\Vert }^3)$$

$$=\psi (0)+\frac{1}{2\alpha }{\mathbf{\theta }}^{\top }{X}^{\top }X\mathbf{\theta }+O(\Vert \mathbf{\theta }{\Vert }^3).$$

$$\log p(y|\mathbf{\theta })$$

$$={\tilde{\mathbf{\theta }}}^{\top }{X}^{\top }X\mathbf{\theta }-\psi (\mathbf{\theta })$$

$$={\tilde{\mathbf{\theta }}}^{\top }{X}^{\top }X\mathbf{\theta }-\left\{\psi (0)+\frac{1}{2\alpha }{\mathbf{\theta }}^{\top }{X}^{\top }X\mathbf{\theta }+O(\Vert \mathbf{\theta }{\Vert }^3)\right\}$$

$$=-\frac{1}{2\alpha }(\mathbf{\theta }-\alpha \tilde{\mathbf{\theta }}{)}^{\top }{X}^{\top }X(\mathbf{\theta }-\alpha \tilde{\mathbf{\theta }})+\frac{\alpha }{2}{\tilde{\mathbf{\theta }}}^{\top }\tilde{\mathbf{\theta }}-\psi (0)+O(\Vert \mathbf{\theta }{\Vert }^3).$$

Full text

Keywords: Exponential family, Generalized linear regression, Information geometry, Sparse modeling

1 Introduction

We propose sparse estimation methods for generalized linear models (GLM). One of the proposed methods is based on least angle regression (LARS) [7] and is described in terms of information geometry. The main features of our approach are i) we use an approximation of a statistical model and do not use the statistical model itself, and ii) the proposed methods are calculated exactly, which allows us to compute the estimators efficiently. A few extensions of LARS, which are based on information, Riemannian, and differential geometry, have been proposed in the literature such as [10] and [4], for example. The existing methods take advantage of a dual structure of a model manifold, which requires computational costs. Our method utilizes a part of the dual structure and uses the original LARS algorithm in the tangent space. The proposed method enables us to compute the estimator easily. Furthermore, we demonstrate that least absolute shrinkage and selection operator (LASSO) [17] for the normal linear model is also available in the tangent space.

Sparse modeling has been extensively investigated in this two decades. As a representative method, LASSO has motivated many researchers in statistics, machine learning, and other fields. LASSO was proposed as an estimation and variable-selection method for the normal linear model. It minimizes the $l_{1}$-regularized least square with a tuning parameter. Various generalizations have been proposed for other problems. For example, [16] and [18] treat the generalized linear regression and Gaussian graphical models, respectively. See also [9].

LARS was proposed for the same problem as LASSO. The LARS algorithm is very efficient, and can also compute the LASSO estimator if a minor change is added. The LARS algorithm uses only correlation coefficients between the response and explanatory variables. The algorithm is therefore described in terms of Euclidean geometry.

Information geometry is a Riemannian-geometric framework for statistics and other fields [1, 2, 3, 5, 13, 15]. In this framework, we treat a statistical model as a Riemannian manifold and take advantage of its geometrical properties for estimation, test, and other tasks. Each probability distribution is treated as a point in the manifold. For example, estimation problem for the generalized linear regression can be described in terms of the geometry. The GLM is treated as a manifold and an estimator assigns a point therein to an observed data. The maximum likelihood estimator (MLE) uses a kind of projection.

Some extensions of LARS have been proposed based on the information geometry of the exponential family of distributions. [10] and [4] proposed different extensions of LARS, that take advantage of the dual structure of the model manifold. Their works are theoretically natural and can be extended to other models than the GLM [11, 12]. However, the existing methods require many iterations of approximation computation, which is inevitable for treating more complicated objects than Euclidean space. For example, [4] treated many tangent spaces that each to an estimate while our methods use only one tangent space. [4] stated “DGLARS method may be computationally more expensive than other customized techniques” for the $l_{1}$-regularization method. We aim to provide as an efficient method as the $l_{1}$-regularization for the GLM. It should be noted that our approach is different from that of existing methods. We approximate the model manifold with only one tangent space and use the exact computation of LARS in the space. This approximation is natural from the viewpoint of information geometry. The usefulness of our idea is validated by numerical experiments. One advantage of our methods is that they do not require additional implementation because we can use existing packages.

The rest of this paper is organized as follows. In section 2, we introduce our problem and the related works. In Section 3, we propose a sparse estimation method based on LARS. Furthermore, LASSO-type estimators are also proposed. We compare our methods with the $l_{1}$-regularization for the GLM by performing numerical experiments in section 4. Finally, our conclusions are presented in section 5. Lemmas and remarks are presented in Appendix A.

2 Problem and Related Method

In subsection 2.1, we formulate the problem and introduce our notation. In subsections 2.2 and 2.3, we briefly describe the LARS algorithm and the LASSO estimators, respectively.

2.1 Problem and Notation

We consider the generalized linear regression, which is an estimation problem of an exponential family of probability distributions [1, 6, 14]. In this 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

[TABLE]

where $a$ is the index indicating $a$-th sample, $h:\mathbb{R}\rightarrow\mathbb{R}$ is a link function, $n$ is the sample size, $d$ is the number of the explanatory variables, and ${\bm{\theta}}=(\theta^{i})$ is the parameter to be estimated. Let $X=(x^{a}_{i})$ be the design matrix, which is an $(n\times d)$-matrix. Let ${\bm{y}}=(y^{a})$ and ${\bm{\mu}}=(\mu^{a})$ 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 current study, we focus on the canonical link function, which results in useful properties of the exponential family. Our method is based on this assumption.

In terms of probability distributions, the aforementioned problem corresponds to estimation for an exponential family of distributions,

[TABLE]

where $\psi:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is a potential function. Our notation takes after that of [1] and [6]. For example, this formulation includes logistic and Poisson regressions.

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

Herein, 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^{a}_{i}=0$ and $\sum_{a=1}^{n}(x^{a}_{i})^{2}=1$ for $i=1,2,\dots,d$. Furthermore, we assume that column vectors of $X$ are linearly independent.

2.2 LARS

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

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 ${\bm{\theta}}=0$ to the MLE $\hat{\bm{\theta}}_{\mathrm{MLE}}$ of the full model. The full model refers to the linear model including all the explanatory variables. The MLE $\hat{\bm{\theta}}_{\mathrm{MLE}}$ is determined by the design matrix $X$ and the response $\bm{y}$. The detailed algorithm of LARS is shown in Algorithm 1, where $\hat{\bm{\theta}}_{(k)}$ is $k$-th estimate the algorithm outputs. After $d$ iterations, LARS outputs a sequence of the estimates $\hat{\bm{\theta}}_{(0)},\hat{\bm{\theta}}_{(1)},\dots,\hat{\bm{\theta}}_{(d)}$.

The LARS algorithm is presented in Figure 1. Figures 1 and 1 indicate the estimator’s move and the residual’s move, respectively, in the parameter space $\mathbb{R}^{d}$ when $d=2$. The LARS estimator i) selects an element of the parameter that forms a least angle between the residual ${\bm{r}}_{(k)}=\hat{\bm{\theta}}_{\mathrm{MLE}}-\hat{\bm{\theta}}_{(k)}$ 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.

2.3 LASSO

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

[TABLE]

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

LASSO can be applied to the GLM as the $l_{1}$-regularized MLE, which is the minimization problem

[TABLE]

For example, see [16].

3 The Proposed Methods

Our main idea, which is very simple, is to run the LARS algorithm in the tangent space of the model manifold. Although the idea appears to be extremely simple, it works well as is illustrated in sections 3 and 4.

In subsection 3.1, we introduce information geometry used herein. In subsection 3.2, we propose LARS in tangent space, which is an extension of the original LARS to the GLM. The proposed method is identical to the original LARS when applied to the normal linear model. In subsection 3.3, we propose other methods that are related with LASSO. Subsection 3.4 explains the difference between the proposed and the existing methods.

3.1 Information Geometry

We briefly introduce some tools from information geometry, including model manifold, tangent space, and exponential map (Figure 2). For details, see [1, 2, 3, 5, 13, 15].

In the generalized linear regression, we need to select one distribution from the exponential family (1). The parameter $\bm{\theta}$ works as a coordinate system in the manifold $\cal M$.

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

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

[TABLE]

where $l({\bm{\theta}})=\log p(y|\,{\bm{\theta}})$ 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

[TABLE]

In the generalized linear regression, the Fisher metric $G$ at $T_{0}{\cal M}$ is proportional to the correlation matrix $X^{\top}X$ of the explanatory variables, that is, $G=cX^{\top}X$ for some $c>0$. This is why we use the tangent space at the origin. For details, see subsection A.1.

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

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

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

3.2 LARS in Tangent Space

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

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

In the estimation step of the proposed method, we run the original LARS algorithm in the tangent space $T_{0}{\cal M}$ as if the response is $\hat{\bm{y}}$. LARS outputs a sequence $\{\hat{\bm{\theta}}_{(0)},\hat{\bm{\theta}}_{(1)},\dots,\hat{\bm{\theta}}_{(d)}\}$ of the model parameter $\bm{\theta}$. As shown in Figure 1, the LARS estimator $\hat{\bm{\theta}}$ can be regarded as moving from the origin to the MLE $\hat{\bm{\theta}}_{\mathrm{MLE}}$ of the full model. At the same time, however, the residual ${\bm{r}}(\hat{\bm{\theta}}):=\hat{\bm{\theta}}_{\mathrm{MLE}}-\hat{\bm{\theta}}$ of the estimator $\hat{\bm{\theta}}$ is moving from the MLE $\hat{\bm{\theta}}_{\mathrm{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. The residual ${\bm{r}}(\hat{\bm{\theta}})$ moves, not the estimator $\hat{\bm{\theta}}$. Note that Algorithm 1 in subsection 2.2 is actually described from the latter perspective.

LARS in Tangent Space (TLARS)

LARS in tangent space (TLARS) is given as follows:

Calculate the MLE $\hat{\bm{\theta}}_{\mathrm{MLE}}$ of the full model. 2. 2.

Run the LARS algorithm for the design matrix $X$ and the response $\hat{\bm{y}}=X\hat{\bm{\theta}}_{\mathrm{MLE}}$. 3. 3.

Using the sequence $\{\hat{\bm{\theta}}_{(0)},\hat{\bm{\theta}}_{(1)},\dots,\hat{\bm{\theta}}_{(d)}\}$ made by LARS, the result is the sequence $\{p(\cdot|\,\hat{\bm{\theta}}_{(0)}),p(\cdot|\,\hat{\bm{\theta}}_{(1)}),\dots,p(\cdot|\,\hat{\bm{\theta}}_{(d)})\}$.

As a special case, the proposed method coincides with the original LARS when we consider the normal linear regression with a known variance. Note that TLARS is as computationally efficient as LARS although it solves the estimation problem of the GLM. Furthermore, we can use existing packages of LARS for the computation of TLARS.

3.3 LASSO in Tangent Space

We propose two estimation methods. One is a LASSO-type modification of TLARS and the other is an approximation of the $l_{1}$-regularization for the GLM (2).

LASSO in Tangent Space 1 (TLASSO1)

By modifying the LARS algorithm so that it outputs the LASSO estimator [7], we can use LASSO in the tangent space $T_{0}{\cal M}$. LASSO in tangent space (TLASSO1) is formally defined as a minimization problem

[TABLE]

which implies that we use the design matrix $X$ and the response $\hat{\bm{y}}=X\hat{\bm{\theta}}_{\mathrm{MLE}}$ in the ordinary LASSO. This corresponds to the LASSO modification of TLARS.

LASSO in Tangent Space 2 (TLASSO2)

Another LASSO-type method is a direct approximation of (2). TLASSO2 is defined as

[TABLE]

where $\alpha=1/(h^{-1})^{\prime}(0)$ and $\tilde{\bm{\theta}}$ satisfies $X^{\top}X\tilde{\bm{\theta}}=X^{\top}{\bm{y}}$. Since the column vectors of the design matrix $X$ are assumed to be linearly independent, $\tilde{\bm{\theta}}$ uniquely exists. Problem (4) is LASSO for the normal linear regression with the design matrix $X$ and the response $\alpha X\tilde{\bm{\theta}}$. TLASSO2 (4) is an approximation of (2). In fact, using $\tilde{\bm{\theta}}$ and $\alpha$, the log-likelihood is approximated as follows (see subsection A.2):

[TABLE]

Note that $\alpha\tilde{\bm{\theta}}$ is an approximation of the MLE $\hat{\bm{\theta}}_{\mathrm{MLE}}$.

3.4 Remarks on Other Information-Geometrical Methods

We briefly compare TLARS with two existing methods that are extensions of LARS based on information geometry. One is bisector regression (BR) by [10] 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 attempts to form 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{\bm{\theta}}_{\mathrm{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{\bm{\theta}}_{\mathrm{MLE}})$ in $\cal M$. Accordingly, the tangent space at the current estimator moves, meaning that we treat the tangent spaces at many points in $\cal 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}{\cal M}$.

4 Numerical Examples

We present results of numerical examples and compare our methods with a related method. In detail, we compare four methods in the logistic regression setting: LARS in Tangent Space (TLARS), LASSO in Tangent Space (TLASSO1 and 2), and the $l_{1}$-regularized maximum likelihood estimation for the GLM (L1).

Our methods do not require an extra implementation since the LARS algorithm has already been implemented in the lars package of the software R. Using R, we only needed glm() for calculating the MLE and the lars package for the proposed methods. For the computation of $l_{1}$-regularization, we used the glmnet package [8].

4.1 Real Data

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

We report 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 the 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 the results for our methods, especially, for TLARS and TLASSO1.

4.2 Numerical Experiments

We performed numerical experiments of logistic regression. The topic is three-fold: generalization, parameter estimation, and model selection. The results are presented in Table 1. Values in bold 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 ${\bm{\theta}}_{0}$ of the parameter ${\bm{\theta}}$, and the sample size $n$. For each of $m$ trials, we made the design matrix $X$ using the rnorm() function in R. Furthermore, we made the response $\bm{y}$ based on $X$ and ${\bm{\theta}}_{0}$, that is, elements of $\bm{y}$ have different Bernoulli distributions. The four methods were applied to $({\bm{y}},X)$.

For selecting one model and one estimate from a sequence of parameter estimates, we used AIC and BIC:

[TABLE]

where $d^{\prime}$ is the dimension of the parameter of the model under consideration. For a sequence $(\hat{\bm{\theta}}_{(k)})$ made by each of the four methods, let $I_{(k)}=\{i|\,\hat{\theta}_{(k)}^{i}\not=0\}$ and $\hat{\bm{\theta}}_{\mathrm{MLE}}^{(k)}$ the MLE of the model ${\cal M}_{(k)}=\{p(\cdot|{\bm{\theta}})|\,\theta^{j}=0\,\,(j\not\in I_{(k)})\}$. We call (5) with $\hat{{\bm{\theta}}}=\hat{\bm{\theta}}_{\mathrm{MLE}}^{(k)}$ AIC1, and (5) with $\hat{{\bm{\theta}}}=\hat{\bm{\theta}}_{(k)}$ AIC2. Similarly, (6) with $\hat{{\bm{\theta}}}=\hat{\bm{\theta}}_{\mathrm{MLE}}^{(k)}$ is BIC1, and (6) with $\hat{{\bm{\theta}}}=\hat{\bm{\theta}}_{(k)}$ is BIC2.

For evaluating the generalization error of the four methods, we newly made $m$ observations $\{({\bm{y}}^{l}_{1},X^{l}_{1}),\dots,({\bm{y}}^{l}_{m},X^{l}_{m})\}$ in $l$-th trial ($l=1,2,\dots,m$). We computed the difference between $({\bm{y}}^{l}_{1},\dots,{\bm{y}}^{l}_{m})$ and $m$ predictions by each of the methods. The “Generalization” columns of Table 1 report the average prediction error over $m$ trials; a smaller value is better.

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

In the “Parameter estimation” columns, each value means the average of $\|\hat{\bm{\theta}}-{\bm{\theta}}_{0}\|_{2}^{2}$ of the selected estimate $\hat{\bm{\theta}}$; a smaller value is better.

In Table 1, we report the results of three cases. We used $m=10,000$ for all cases except case C2, for which we set $m=1,000$.

In case A, we set $d=10$ and ${\bm{\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, ${\bm{\theta}}_{0}=(10,10,-10,-10,0,0,0,0,0,0)^{\top}$.

Case B is the case of $d=10$ and ${\bm{\theta}}_{0}=(10,10,0,0,0,$ $0,0,0,0,0)^{\top}$ with the relation ${\bm{x}}_{3}={\bm{x}}_{2}+{\bm{\epsilon}}$, where ${\bm{x}}_{2}$ and ${\bm{x}}_{3}$ are the second and third columns of the design matrix $X$, respectively, and $\bm{\epsilon}$ 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 ${\bm{\theta}}_{0}$, the vector of the length 50 with ten $10$s, ten $-10$s, and thirty [math]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.

[FIGURE:]

5 Conclusion

We proposed the sparse estimation methods as an extension of LARS for the GLM. The methods take advantage of the tangent space at the origin, which is a natural approximation of the model manifold. The proposed methods are computationally efficient because the problem is approximated by the normal linear regression. The numerical experiments showed that our idea worked well by comparing the proposed methods with the $l_{1}$-regularization for the GLM. One of our future works is to evaluate our methods theoretically. Furthermore, we will apply tools developed for LARS and LASSO to TLARS and TLASSO, for example, screening and post-selection inference.

Appendix A Lemmas and Remarks

Some lemmas and remarks are presented. We use well-known facts on the exponential family, the GLM, and information geometry. For details, see [1, 2, 3, 5, 6, 13, 14, 15].

As introduced in subsection 3.3, $\tilde{\bm{\theta}}$ satisfies $X^{\top}X\tilde{\bm{\theta}}=X^{\top}{\bm{y}}$ and $\alpha=1/\tilde{h}(0)$, where $h$ is the link function and $\tilde{h}(0)=(h^{-1})^{\prime}(0)$. Let $\eta_{j}$ be $j$-th element of the expectation parameter ${\bm{\eta}}({\bm{\theta}})=\mathrm{E}_{\bm{\theta}}[X^{\top}{\bm{y}}]$. Letting ${\bm{\mu}}({\bm{\theta}})=\mathrm{E}_{\bm{\theta}}[{\bm{y}}]$, it holds ${\bm{\eta}}({\bm{\theta}})=X^{\top}{\bm{\mu}}({\bm{\theta}})$ and $\mu^{a}({\bm{\theta}})=h^{-1}(\sum_{i=1}^{d}x_{i}^{a}\theta^{i})\,\,(a=1,2,\dots,n)$.

A.1 Metric at Tangent Space and Correlation Between Explanatory Variables

We show that the Fisher metric $G=(g_{ij})$ at the tangent space $T_{0}{\cal M}$ is proportional to the correlation matrix $X^{\top}X$ of the explanatory variables (Lemma 2). To avoid confusion, in this subsection, we use $G(0)=(g_{ij}(0))$ for the metric in $T_{0}{\cal M}$ and $G({\bm{\theta}})=(g_{ij}({\bm{\theta}}))$ for the metric in the tangent space at $p(\cdot|\,{\bm{\theta}})$.

Lemma 1**.**

It holds that

[TABLE]

Proof.

Since it is known that $\eta_{i}=\partial\psi/\partial\theta^{i}$,

[TABLE]

∎

Lemma 2**.**

$G(0)=cX^{\top}X$* for some $c>0$.*

Proof.

It is known that the metric $g_{ij}$ is derived from the potential function $\psi$: $g_{ij}({\bm{\theta}})=\partial_{i}\partial_{j}\psi({\bm{\theta}})$. Therefore, it holds

[TABLE]

where $\tilde{h}$ is the derivative of $h^{-1}$. Letting ${\bm{\theta}}=0$ and $c=\tilde{h}(0)$, we have $G(0)=cX^{\top}X$. Since both $G(0)$ and $X^{\top}X$ are known to be positive definite, $c$ is a positive constant. ∎

Note that $c$ is common to all $i,j$ and $a$ in the proof. Hence, the tangent space $T_{0}{\cal M}$ at the origin ${\bm{\theta}}=0$ is selected as the space in which LARS runs.

A.2 Approximations of the Likelihood and MLE

We approximate the log-likelihood and the MLE of the GLM (1). Lemma 3 implies that $\alpha\tilde{{\bm{\theta}}}$ is an approximation of the MLE $\hat{\theta}$

Lemma 3**.**

The log-likelihood is expanded as

[TABLE]

Proof.

Using $\alpha=1/\tilde{h}(0)$ and Lemmas 1 and 2, the potential function $\psi$ is expanded as follows:

[TABLE]

At the last equal sign, we used ${\bm{1}}^{\top}X=0$ since each column vector of $X$ is assumed to be normalized. Therefore,

[TABLE]

∎

A.3 e-Exponential Map

In Riemannian geometry, a point in a tangent space is mapped to a manifold via an exponential map, which is defined using a geodesic. A geodesic in a manifold corresponds to a straight line in Euclidean space. When we consider an exponential map, we need to introduce not only a metric but also a connection, which determines flatness and straightness in a manifold. In section 3, we implicitly introduced the e-connection. From the viewpoint of the e-connection, each curve of $\theta^{i}$-axis is an e-geodesic in $\cal M$.

For a manifold $\cal M$ and a point $p\in{\cal M}$, an exponential map $f$ at $p$ is formally defined as follows. First, we consider the geodesic $\gamma_{v}(t)$ for $v\in T_{p}{\cal M}$ which satisfies $\gamma_{v}(0)=p$ and $\mathrm{d}\gamma_{v}(t)/\mathrm{d}t|_{t=0}=v$. Here, the parameter $t$ moves in an interval including [math]. Note that, given a connection, the geodesic $\gamma_{v}$ locally exists and is uniquely determined. The exponential map $f$ is $f:T_{p}{\cal M}\rightarrow{\cal M}$ and $f(v)=\gamma_{v}(1)$ for $v\in D\subset T_{p}{\cal M}$, where $D=\{v\in T_{p}{\cal M}|\,\gamma_{v}(1)\,\,\mathrm{exists}\}$.

In general, an exponential map is not necessarily easy to treat. For example, the domain $D$ of such a map is known as a star-shaped domain and does not coincide with a whole tangent space. However, our exponential map $\mathrm{Exp_{0}}:T_{0}{\cal M}\rightarrow M$ has a useful property. The domain of $\mathrm{Exp_{0}}$ is a whole $T_{0}{\cal M}$ and the range is a whole $\cal M$.

Lemma 4**.**

The map $\mathrm{Exp}_{0}:T_{0}{\cal M}\rightarrow M$ defined in subsection 3.1 is the e-exponential map for a manifold of the GLM. Furthermore, $\mathrm{Exp}_{0}$ is a bijection from the tangent space $T_{0}{\cal M}$ to the manifold $\cal M$.

Proof.

For $v=\sum_{i=1}^{d}v^{i}\partial_{i}\in T_{0}{\cal M}$, the value of the map is $\mathrm{Exp_{0}}(v)=p(\cdot|\,{\bm{v}})$, where ${\bm{v}}=(v^{i})$. It is known that the e-geodesic $\gamma(t)$ satisfying $\gamma(0)=p(\cdot|\,0)$ and $\mathrm{d}\gamma(t)/\mathrm{d}t|_{t=0}=v\in T_{0}{\cal M}$ is represented as $\gamma(t)=p(\cdot|\,t{\bm{v}})$. Therefore, $\mathrm{Exp_{0}}(v)=p(\cdot|\,{\bm{v}})=\gamma(1)$, which means that $\mathrm{Exp}_{0}$ is the e-exponential map.

Since ${\cal M}=\{p(\cdot|\,{\bm{\theta}})|\,{\bm{\theta}}\in\mathbb{R}^{d}\}$, the e-exponential map is defined on a whole $T_{0}{\cal M}$. For ${\bm{\theta}}\in\mathbb{R}^{d}$, $w=\sum_{i=1}^{d}\theta^{i}\partial_{i}$ is in $T_{0}{\cal M}$ and $\mathrm{Exp_{0}}(w)=p(\cdot|\,{\bm{\theta}})$, which implies that the e-exponential map is a surjection. Furthermore, if $v,w\in T_{0}{\cal M}$ are different, $\mathrm{Exp_{0}}(v)\not=\mathrm{Exp_{0}}(w)$ because the column vectors of $X$ are assumed to be linearly independent. ∎