# A Generic Path Algorithm for Regularized Statistical Estimation

Hua Zhou
Department of Statistics
North Carolina State University
Raleigh, NC 27695-8203
E-mail<EMAIL_ADDRESS>
Yichao Wu
Department of Statistics
North Carolina State University
Raleigh, NC 27695-8203
E-mail<EMAIL_ADDRESS>
###### Abstract

Regularization is widely used in statistics and machine learning to prevent
overfitting and gear solution towards prior information. In general, a
regularized estimation problem minimizes the sum of a loss function and a
penalty term. The penalty term is usually weighted by a tuning parameter and
encourages certain constraints on the parameters to be estimated. Particular
choices of constraints lead to the popular lasso, fused-lasso, and other
generalized $l_{1}$ penalized regression methods. Although there has been a
lot of research in this area, developing efficient optimization methods for
many nonseparable penalties remains a challenge. In this article we propose an
exact path solver based on ordinary differential equations (EPSODE) that works
for any convex loss function and can deal with generalized $l_{1}$ penalties
as well as more complicated regularization such as inequality constraints
encountered in shape-restricted regressions and nonparametric density
estimation. In the path following process, the solution path hits, exits, and
slides along the various constraints and vividly illustrates the tradeoffs
between goodness of fit and model parsimony. In practice, the EPSODE can be
coupled with AIC, BIC, $C_{p}$ or cross-validation to select an optimal tuning
parameter. Our applications to generalized $l_{1}$ regularized generalized
linear models, shape-restricted regressions, Gaussian graphical models, and
nonparametric density estimation showcase the potential of the EPSODE
algorithm.
Keywords: Gaussian graphical model, generalized linear model, lasso, log-
concave density estimation, ordinary differential equations, quasi-
likelihoods, regularization, shape restricted regression, solution path

##  1 Introduction

Regularization is a frequently used framework in statistics. Examples include
the lasso regression (Tibshirani, 1996; Chen et al., 2001) and the $l_{1}$
penalized generalized linear models (GLMs) among many others. For both the
lasso and the $l_{1}$ penalized GLMs, efficient solution path algorithms have
been proposed to ease the tuning of the regularization parameter (Osborne et
al., 2000; Efron et al., 2004; Park and Hastie, 2007). Yet extension to other
more general settings is nontrivial and has been an active research area.

In this article, we consider a general regularization framework

$\displaystyle\min_{\boldsymbol{\beta}\in\mathbb{R}^{p}}f(\boldsymbol{\beta})+\rho\|\boldsymbol{V}\boldsymbol{\beta}-\boldsymbol{d}\|_{1}+\rho\|\boldsymbol{W}\boldsymbol{\beta}-\boldsymbol{e}\|_{+},$
(1)

for which we propose an efficient exact path solver based on ordinary
differential equations (EPSODE). Here $f:\mathbb{R}^{p}\mapsto\mathbb{R}$ can
be any convex, smooth loss function of $\boldsymbol{\beta}\in\mathbb{R}^{p}$,
where $p>0$ is the dimensionality of the parameters. For any vector
$\boldsymbol{v}=(v_{i})$, $\|\boldsymbol{v}\|_{1}=\sum_{i}|v_{i}|$ denotes its
$l_{1}$ norm and $\|\boldsymbol{v}\|_{+}=\sum_{i}\max\\{v_{i},0\\}$ is the sum
of positive parts of its components. The EPSODE provides the exact solution
path to (1) as the tuning parameter $\rho$ varies.

###  1.1 Generality of (1)

The generality of (1) is two-fold. First $f$ can by any convex loss function.
For example, it can be the negative log-likelihood function of GLMs, negative
quasi-likelihood, the exponential loss function of the AdaBoost (Friedman et
al., 2000), or many other frequently used loss functions in statistics and
machine learning. Second we allow $\boldsymbol{V}$ and $\boldsymbol{W}$ to be
any regularization matrices of $p$ columns. This leads to broad applications.
In particular, the first regularization term
$\rho\|\boldsymbol{V}\boldsymbol{\beta}-\boldsymbol{d}\|_{1}$ encourages
equality constraints among parameters $\boldsymbol{\beta}$. When $\rho$ is
large enough, the minimizer $\boldsymbol{\beta}(\rho)$ of (1) satisfies
$\boldsymbol{V}\boldsymbol{\beta}(\rho)=\boldsymbol{d}$. For instance, when
$\boldsymbol{V}$ is the identity matrix and $\boldsymbol{d}={\bf 0}$, it
recovers the well-known lasso regression (Tibshirani, 1996; Chen et al.,
2001), which encourages sparsity of the estimates. When

$\displaystyle\boldsymbol{V}$
$\displaystyle=\left(\begin{array}[]{rrrrr}-1&1&\\\ &\ddots&\ddots\\\
&&-1&1\end{array}\right)$

and $\boldsymbol{d}={\bf 0}$, it corresponds to the fused-lasso penalty
(Tibshirani et al., 2005), which leads to smoothness among neighboring
regression coefficients. As we will show later, more complicated equality
constraints can be incorporated with properly designed $\boldsymbol{V}$ and
$\boldsymbol{d}$. On the other hand, the second regularization term
$\rho\|\boldsymbol{W}\boldsymbol{\beta}-\boldsymbol{e}\|_{+}$ enforces
regularization by inequality relations among regression coefficients. For
large enough $\rho$, the minimizer $\boldsymbol{\beta}(\rho)$ satisfies
$\boldsymbol{W}\boldsymbol{\beta}(\rho)\leq\boldsymbol{e}$. For instance,
setting $\boldsymbol{W}$ as the negative identity matrix and
$\boldsymbol{e}={\bf 0}$ encourages nonnegativity of the estimates, as
required in nonnegative least squares problems (Lawson and Hanson, 1987). In
the isotonic regression (Robertson et al., 1988; Silvapulle and Sen, 2005),
the estimates have to be nondecreasing. This can be achieved by the
regularization matrix

$\displaystyle\boldsymbol{W}$
$\displaystyle=\left(\begin{array}[]{rrrrr}1&-1&\\\ &\ddots&\ddots\\\
&&1&-1\end{array}\right)$

and $\boldsymbol{e}={\bf 0}$. More complicated constraints that occur in
shape-restricted regression and nonparametric regressions also can be
incorporated as we demonstrate in later examples.

In certain applications, both equality and inequality regularizations are
required. In that case, as shown in Section 2, at a large but finite $\rho$,
the minimizer $\boldsymbol{\beta}(\rho)$ coincides with the solution to the
following constrained optimization problem

$\displaystyle\min$ $\displaystyle f(\boldsymbol{\beta})$ s.t.
$\displaystyle\boldsymbol{V}\boldsymbol{\beta}=\boldsymbol{d}\mbox{ and
}\boldsymbol{W}\boldsymbol{\beta}\leq\boldsymbol{e}.$

Consequently EPSODE solves the linearly constrained estimation problem (1.1)
as a by-product. In this case, path following commences from the unconstrained
solution $\text{argmin}f(\boldsymbol{\beta})$ and ends at the constrained
solution to (1.1).

###  1.2 Previous Work

Several path algorithms have been devised for special cases of the general
regularization problem (1). For example, the homotopy method (Osborne et al.,
2000) and the least angle regression (LARS) procedure (Efron et al., 2004)
handle lasso penalized least squares problem. The solution path generated is
piecewise linear and illustrates the tradeoffs between goodness of fit and
sparsity. Rosset and Zhu (2007) give sufficient conditions for a solution path
to be piecewise linear and expand its applications to a wider range of loss
and penalty functions. Recently Tibshirani and Taylor (2011) devise a dual
path algorithm for generalized $l_{1}$ penalized least squares problems, which
is problem (1) with $f$ quadratic but without the second inequality
regularization term. Zhou and Lange (2011) consider (1) in full generality for
quadratic $f$. All these work concerns regularized linear regression for which
the solution path is piecewise linear. Several attempts have been made to path
following for regularized GLMs for which the solution path is no longer
piecewise linear. Park and Hastie (2007) propose a predictor-corrector
approach to approximate the lasso path for GLMs. Wu (2011) presents an
ordinary differential equation-based path algorithm which delivers the exact
path for lasso penalized GLMs. Friedman (2008) derives an approximate path
algorithm for any convex loss regularized by a separable, but not necessarily
convex penalty. Here a penalty function is called separable if its Hessian
matrix is diagonal. The separability restriction on the penalty term excludes
many important problems encountered in real applications.

Our proposed approach generalizes previous work in several aspects. First, it
works for any convex loss (or criterion) function. Second, it allows for any
type of regularization in terms of linear functions of parameters, equality or
inequality. Equality constrained regularizations include lasso, fused-lasso
and generalized $l_{1}$ penalty for example. Inequality constrained
regularizations are required in shape-restricted regression and nonparametric
log-concave density estimation. Last but not least, it is an exact path
algorithm.

###  1.3 A Motivating Example

For illustration, we consider a merger and acquisition (M&A) data set studied
in (Fan et al., 2011). This data set constitutes $n=1,371$ US companies with a
binary response variable indicating whether the company becomes a leveraged
buyout (LBO) target ($y_{i}=1$) or not ($y_{i}=0$). Seven covariates (1. cash
flow, 2. cash, 3. long term investment, 4. market to book ratio, 5. log market
equity, 6. tax, 7. return on S&P 500 index) are recorded for each company.
There have been intensive studies on the effects of these factors on the
probability of a company being a target for strategic mergers. Exploratory
analysis using linear logistic regression shows no significance in most
covariates.

To explore the possibly nonlinear effects of these quantitative covariates,
the varying-coefficient model (Hastie and Tibshirani, 1993) can be adopted
here. We discretize each predictor into, say, 10 bins and fit a logistic
regression. The first bin of each predictor is used as the reference level and
effect coding is applied to each discretized covariate. The circles (o) in
Figure 1 denote the estimated coefficients for each bin of each predictor and
hint at some interesting nonlinear effects. For instance, the chance of being
an LBO target seems to monotonically decease with market-to-book ratio and be
quadratic as a function of log market equity. Regularization can be utilized
to borrow strength between neighboring bins and gear solution towards clearer
patterns. To illustrate the flexibility of the regularization scheme (1), we
apply cubic trend filtering to 5 covariates (cash flow, cash, long term
investment, tax, return on S&P 500 index), impose the monotonicity (non-
increasing) constraint on the ‘market-to-book ratio’ covariate, and enforce
the concavity constraint on the ‘log market equity’ covariate. This can be
achieved by minimizing a regularized negative logistic log-likelihood of form

$\displaystyle-l(\boldsymbol{\beta}_{1},\ldots,\boldsymbol{\beta}_{7})+\rho\sum_{j\neq
4,5}\|\boldsymbol{V}_{j}\boldsymbol{\beta}_{j}\|_{1}+\rho\sum_{j=4,5}\|\boldsymbol{W}_{j}\boldsymbol{\beta}_{j}\|_{+},$

where $\boldsymbol{\beta}_{j}$ is the vector of regression coefficients for
the $j$-th discretized covariate. The matrices in the regularization terms are
specified as

$\displaystyle\boldsymbol{V}_{j}$
$\displaystyle=\left(\begin{array}[]{rrrrrrrr}-1&2&-1\\\ 1&-4&6&-4&1\\\
&1&-4&6&-4&1&\\\ &&&\ddots&\ddots&\ddots\\\ &&1&-4&6&-4&1&\\\
&&&&&-1&2&-1\end{array}\right)\text{ for }j=1,2,3,6,7,$
$\displaystyle\boldsymbol{W}_{4}$
$\displaystyle=\left(\begin{array}[]{rrrrrr}-1&1&\\\ &-1&1&\\\
&&\ddots&\ddots\\\ &&&-1&1&\\\ &&&&-1&1\end{array}\right),\mbox{ and }$
$\displaystyle\boldsymbol{W}_{5}$
$\displaystyle=\left(\begin{array}[]{rrrrrrr}1&-2&1\\\ &1&-2&1\\\
&&\ddots&\ddots&\ddots\\\ &&&1&-2&1&\\\ &&&&1&-2&1\end{array}\right).$

The equality constraint regularization matrix $\boldsymbol{V}_{j}$,
$j=1,2,3,6,7$, penalizes the fourth order finite differences between the bin
estimates. Thus, as $\rho$ increases, the coefficient vectors of covariates
1-3,6-7 tend to be piecewise cubic with two ends being linear, mimicking the
natural cubic spline. This is one example of the polynomial trend filtering
(Kim et al., 2009; Tibshirani and Taylor, 2011). Similar to semi-parametric
regressions, regularizations in polynomial trend filtering ‘let the data speak
for themselves’. In contrast, the bandwidth selection in semi-parametric
regressions is replaced by parameter tuning in regularizations. The number and
locations of knots are automatically determined by tuning parameter which is
chosen according to model selection criteria. In a similar fashion, the
coefficient vector gradually becomes monotone for covariate ‘market-to-book
ratio’ and concave for covariate ‘log market equity’. In addition, with $\rho$
large enough, we recover the corresponding constrained solution, which are
shown by the crosses (+) on solid lines in Figure 1. As noted above, our exact
path algorithm delivers the whole solution path bridging from the
unconstrained estimates (denoted by o) to the constrained estimates (denoted
by +). For example, the dotted lines in Figure 1 is a snapshot of the solution
at $\rho=0.6539$. Availability of the whole solution path renders model
selection along the path easy. For instance the regularization parameter
$\rho$ can be chosen by minimizing the cross-validation error or other model
selection criteria such as AIC, BIC, or $C_{p}$. Figure 2 displays the
solution path and the AIC and BIC along the path. It shows that both criteria
favor the fully regularized solution, namely the constrained estimates. The
whole solution path is obtained within seconds on a laptop using a Matlab
implementation of EPSODE.

The patterns revealed by the regularized estimates match some existing finance
theories. For instance, a company with low cash flow is unlikely to be an LBO
target because low cash flow is hard to meet the heavy debt burden associated
with the LBO. On the other hand, company carrying a high cash flow is likely
to possess a new technology. It is risky to acquire such firms because it is
hard to predict their profitability. The tax reason is obvious from the
regularized estimates. The more tax the company is paying, the more tax
benefits from an LBO. Log of market equity is a measure of company size.
Smaller companies are unpredictable in their profitability and extremely large
companies are unlikely to be an LBO target because LBOs are typically financed
with a large proportion of external debts. Interested readers are referred to
(Shivdasani and Wang, 2009) and references therein for related theories on
LBO.

This illustrative example demonstrates the flexibility of our novel path
algorithm. First, it can be applied to any convex loss function. In this
example, the loss function is the negative log-likelihood of a logistic model.
Second, it works for complicated regularizations like polynomial trend
filtering (equality constraints), monotonicity constraint, and concavity
constraint. More applications will be presented in Section 7 to illustrate the
potential of EPSODE.

$\begin{array}[]{cc}\includegraphics[width=433.62pt]{MandA_estimates}\end{array}$
Figure 1: Snapshots of the path solution to the regularized logistic
regression on the M&A data set.
$\begin{array}[]{cc}\includegraphics[width=180.67499pt]{MandA_solpath}&\includegraphics[width=180.67499pt]{MandA_AICBICpath}\end{array}$
Figure 2: Solution and AIC/BIC paths of the regularized logistic regression on
the M&A data set.

The rest of the paper is organized as follows. Section 2 reviews the exact
penalty method for optimization. Here the connections between constrained
optimization and regularization in statistics are made clear. Section 3
derives in detail the EPSODE algorithm for strictly convex loss function $f$.
Its implementation via the sweep operator and ordinary differential equations
are described in Section 4. An extension of EPSODE for $f$ convex but not
necessarily strictly convex is discussed in Section 5. Section 6 concerns
model selection along the path. Section 7 presents various applications of
EPSODE. Finally, Section 8 discusses the limitations of the path algorithm and
hints at future generalizations.

##  2 Exact Penalty Method for Convex Constrained Optimization

Consider the convex program

$\displaystyle\min$ $\displaystyle f(\boldsymbol{x})$ s.t. $\displaystyle
g_{i}(\boldsymbol{x})=0,1\leq i\leq r$ $\displaystyle
h_{j}(\boldsymbol{x})\leq 0,1\leq j\leq s,$

where the objective function $f$ is convex, equality constraint functions
$g_{i}$ are affine, and the inequality constraint functions $h_{j}$ are
convex. We further assume that $f$ and $h_{j}$ are smooth. Specifically we
require that $f$ and $h_{j}$ are continuously twice differentiable. To fix
notation, differential $df(\boldsymbol{x})$ is the row vector of partial
derivatives of $f$ at $\boldsymbol{x}$ and the gradient $\nabla
f(\boldsymbol{x})$ is the transpose of $df(\boldsymbol{x})$. The Hessian
matrix of $f(\cdot)$ is denoted by $d^{2}f(\boldsymbol{x})$.

Exact penalty method minimizes the function

$\displaystyle{\cal E}_{\rho}(\boldsymbol{x})$ $\displaystyle=$ $\displaystyle
f(\boldsymbol{x})+\rho\sum_{i=1}^{r}|g_{i}(\boldsymbol{x})|+\rho\sum_{j=1}^{s}\max\\{0,h_{j}(\boldsymbol{x})\\}$
(4)

for $\rho\geq 0$. Classical results (Ruszczyński, 2006, Theorems 6.9 and 7.21)
state that for $\rho$ large enough, the solution to the optimization problem
(4) coincides with the solution to the original constrained convex program
(2). This justifies the exact penalty method as one way to solve constrained
optimization problems.

According to convex calculus (Ruszczyński, 2006, Theorem 3.5), the optimal
point $\boldsymbol{x}(\rho)$ of the function ${\cal E}_{\rho}(\boldsymbol{x})$
is characterized by the necessary and sufficient condition

$\displaystyle{\bf 0}$ $\displaystyle=$ $\displaystyle\nabla
f(\boldsymbol{x})+\rho\sum_{i=1}^{r}s_{i}\nabla
g_{i}(\boldsymbol{x})+\rho\sum_{j=1}^{s}t_{j}\nabla h_{j}(\boldsymbol{x})$ (5)

with coefficients satisfying

$\displaystyle s_{i}\in\begin{cases}\\{-1\\}&g_{i}(\boldsymbol{x})<0\\\
[-1,1]&g_{i}(\boldsymbol{x})=0\\\
\\{1\\}&g_{i}(\boldsymbol{x})>0\end{cases},\hskip 36.135pt\mbox{ and }\hskip
36.135ptt_{j}\in\begin{cases}\\{0\\}&h_{j}(\boldsymbol{x})<0\\\
[0,1]&h_{j}(\boldsymbol{x})=0\\\ \\{1\\}&h_{j}(\boldsymbol{x})>0\end{cases}.$
(6)

The sets defining possible values of $s_{i}$ and $t_{j}$ are the
subdifferentials of the functions $|x|$ and $x_{+}=\max\\{x,0\\}$. For path
following to make sense, we require uniqueness and continuity of the solution
$\boldsymbol{x}(\rho)$ to (4) as $\rho$ varies. The following lemma concerns
the continuity of the solution path and is the foundation of our path
algorithm.

######  Lemma 2.1.

  1. 1.

(Uniqueness) If ${\cal E}_{\rho}$ is strictly convex, then its minimizer
$\boldsymbol{x}(\rho)$ is unique.

  2. 2.

(Continuity) If ${\cal E}_{\rho}$ is strictly convex and coercive over an open
neighborhood of $\rho$, then the minimizer $\boldsymbol{x}(\rho)$ is
continuous at $\rho$.

  3. 3.

(Continuity of $s_{i}$ and $t_{j}$) Furthermore, if the gradients $\\{\nabla
g_{i}(\boldsymbol{x}):g_{i}(\boldsymbol{x})=0\\}\cup\\{\nabla
h_{j}(\boldsymbol{x}):h_{j}(\boldsymbol{x})=0\\}$ of active constraints are
linearly independent at the solution $\boldsymbol{x}(\rho)$ over an open
neighborhood of $\rho$, then the coefficient paths $s_{i}(\rho)$ and
$t_{j}(\rho)$ are unique and continuous at $\rho$.

###### Proof.

The uniqueness of minimum under strict convexity is well-known (Ruszczyński,
2006). For continuity, suppose that the solution $\boldsymbol{x}(\rho)$ is not
continuous at $\rho$. Then there exists $\epsilon>0$ and a sequence
$\rho_{n}\to\rho$ such that
$\|\boldsymbol{x}(\rho_{n})-\boldsymbol{x}(\rho)\|\geq\epsilon$ for all $n$.
Since ${\cal E}_{\rho}$ is coercive, $\boldsymbol{x}(\rho_{n})$ is bounded and
there exists a subsequence of $\boldsymbol{x}(\rho_{n})$ that converges to
some point $\boldsymbol{y}$. Taking limits in the inequality ${\cal
E}_{\rho_{n}}[\boldsymbol{x}(\rho_{n})]\leq{\cal
E}_{\rho_{n}}(\boldsymbol{x})$ shows that ${\cal
E}_{\rho}(\boldsymbol{y})\leq{\cal E}_{\rho}(\boldsymbol{x})$ for all
$\boldsymbol{x}$, i.e., $\boldsymbol{y}=\boldsymbol{x}(\rho)$. This
contradicts with $\|\boldsymbol{y}-\boldsymbol{x}(\rho)\|\geq\epsilon$.
Therefore $\boldsymbol{x}(\rho)$ is continuous at $\rho$. For the continuity
of coefficients, under the linearly independence assumption, $s_{i}(\rho)$ and
$t_{j}(\rho)$ can be uniquely solved by the stationarity condition (5) given
solution vector $\boldsymbol{x}(\rho)$. Therefore continuity of $s_{i}(\rho)$
and $t_{j}(\rho)$ inherits from continuity of $\boldsymbol{x}(\rho)$. ∎

We remark that strict convexity only gives an easy-to-check sufficient
condition for uniqueness and continuity; it is not a necessary condition. A
convex but not strictly convex function can still have a unique minimum. The
absolute value function $|x|$ is such an example. When the loss function $f$
is strictly convex, then ${\cal E}_{\rho}$ is strictly convex for all
$\rho\geq 0$ and by Lemma 2.1 there exists a unique, continuous solution path
$\\{\boldsymbol{x}(\rho):\rho\geq 0\\}$. In Section 3 and 4, we derive the
path algorithm assuming that $f$ is strictly convex. When $f$ is convex but
not strictly convex, e.g., when $n<p$ in the least squares problems, the
solutions at smaller $\rho$ may not be unique. In that case, it is still
possible to obtain a solution path over the region of large $\rho$ where the
minimum of ${\cal E}_{\rho}$ is unique. In Section 5, we extend EPSODE to the
case $f$ is convex but may not be strictly convex. The third statement of
Lemma 2.1 implies that the active constraints ($g_{i}(\boldsymbol{x})=0$ or
$h_{j}(\boldsymbol{x})=0$) with interior coefficients must stay active until
the coefficients hit the end points of the permissible range, which in turn
implies that the solution path is piecewise smooth. This allows us to develop
a path following algorithm based on ODE.

##  3 The Path Following Algorithm

In this article, we specialize to the case where the constraint functions
$g_{i}$ and $h_{j}$ are affine, i.e., the gradient vectors $\nabla
g_{i}(\boldsymbol{x})$ and $\nabla h_{j}(\boldsymbol{x})$ are constant. This
leads to the regularized optimization problem formulated as (1) by defining
$g_{i}$ and $h_{j}$ as constraint residuals
$g_{i}(\boldsymbol{x})=\boldsymbol{v}_{i}^{t}\boldsymbol{x}-d_{i}$ and
$h_{j}(\boldsymbol{x})=\boldsymbol{w}_{j}^{t}\boldsymbol{x}-e_{j}$. In
principle a similar path algorithm can be developed for the general convex
program where the inequality constraint functions $h_{j}$ are relaxed to be
convex. But that is beyond the scope of the current paper. In Sections 3 and
4, we assume that the loss function $f$ is strictly convex. This assumption is
relaxed in Section 5.

Our path following algorithm EPSODE works in a segment-by-segment fashion.
Along the path we keep track of the following index sets determined by signs
of constraint residuals

$\displaystyle{\cal N}_{\text{E}}$
$\displaystyle=\\{i:g_{i}(\boldsymbol{x})=\boldsymbol{v}_{i}^{t}\boldsymbol{x}-d_{i}<0\\},\hskip
36.135pt{\cal
N}_{\text{I}}=\\{j:h_{j}(\boldsymbol{x})=\boldsymbol{w}_{j}^{t}\boldsymbol{x}-e_{j}<0\\}$
$\displaystyle{\cal Z}_{\text{E}}$
$\displaystyle=\\{i:g_{i}(\boldsymbol{x})=\boldsymbol{v}_{i}^{t}\boldsymbol{x}-d_{i}=0\\},\hskip
36.135pt{\cal
Z}_{\text{I}}=\\{j:h_{j}(\boldsymbol{x})=\boldsymbol{w}_{j}^{t}\boldsymbol{x}-e_{j}=0\\}$
(7) $\displaystyle{\cal P}_{\text{E}}$
$\displaystyle=\\{i:g_{i}(\boldsymbol{x})=\boldsymbol{v}_{i}^{t}\boldsymbol{x}-d_{i}>0\\},\hskip
36.135pt{\cal
P}_{\text{I}}=\\{j:h_{j}(\boldsymbol{x})=\boldsymbol{w}_{j}^{t}\boldsymbol{x}-e_{j}>0\\}.$

Along each segment of the path, the set configuration is fixed. This is
implied by the continuity of both the solution and coefficient paths
established in Lemma 2.1. Throughout this article, we call the constraints in
${\cal Z}_{\text{E}}$ or ${\cal Z}_{\text{I}}$ active and others inactive.

Next we derive the ODE for the solution $\boldsymbol{x}(\rho)$ on a fixed
segment. Suppose we are in the interior of a segment. Let
$\boldsymbol{x}(\rho)$ be the solution of (4) indexed by the penalty parameter
$\rho$ and $\boldsymbol{x}(\rho+\Delta\rho)$ the solution when the penalty is
increased by an infinitesimal amount $\Delta\rho>0$. Then the difference
$\Delta\boldsymbol{x}(\rho)=\boldsymbol{x}(\rho+\Delta\rho)-\boldsymbol{x}(\rho)$
should minimize the increase in optimal objective value. That is, to the
second order, $\Delta\boldsymbol{x}$ is the solution to

$\displaystyle\min_{\Delta\boldsymbol{x}}$ $\displaystyle{\cal
E}_{\rho+\Delta\rho}(\boldsymbol{x}+\Delta\boldsymbol{x})-{\cal
E}_{\rho}(\boldsymbol{x})$ $\displaystyle\approx$ $\displaystyle
df(\boldsymbol{x})\cdot\Delta\boldsymbol{x}+\frac{1}{2}\Delta\boldsymbol{x}^{t}\cdot
d^{2}f(\boldsymbol{x})\cdot\Delta\boldsymbol{x}$
$\displaystyle+(\rho+\Delta\rho)\cdot\left[-\sum_{i\in{\cal
N}_{\text{E}}}\boldsymbol{v}_{i}+\sum_{i\in{\cal
P}_{\text{E}}}\boldsymbol{v}_{i}+\sum_{j\in{\cal
P}_{\text{I}}}\boldsymbol{w}_{j}\right]\cdot\Delta\boldsymbol{x}$
$\displaystyle+\Delta\rho\cdot\left[-\sum_{i\in{\cal
N}_{\text{E}}}g_{i}(\boldsymbol{x})+\sum_{i\in{\cal
P}_{\text{E}}}g_{i}(\boldsymbol{x})+\sum_{j\in{\cal
P}_{\text{I}}}h_{j}(\boldsymbol{x})\right]$ s.t.
$\displaystyle\boldsymbol{v}_{i}^{t}\cdot\Delta\boldsymbol{x}=0,i\in{\cal
Z}_{\text{E}},$
$\displaystyle\boldsymbol{w}_{j}^{t}\cdot\Delta\boldsymbol{x}=0,j\in{\cal
Z}_{\text{I}}.$

Note that the active constraints have to be kept active since the set
configuration is fixed along this segment by Lemma 2.1. This is why we have
these two sets of equality constraints. To ease notational burden, we define

$\displaystyle\boldsymbol{H}(\boldsymbol{x})$ $\displaystyle=$ $\displaystyle
d^{2}f(\boldsymbol{x})$ (9) $\displaystyle\boldsymbol{u}_{\bar{\cal Z}}$
$\displaystyle=$ $\displaystyle-\sum_{i\in{\cal
N}_{\text{E}}}\boldsymbol{v}_{i}+\sum_{i\in{\cal
P}_{\text{E}}}\boldsymbol{v}_{i}+\sum_{j\in{\cal
P}_{\text{I}}}\boldsymbol{w}_{j}.$

This leads to the corresponding Lagrange multiplier problem

$\displaystyle\left(\begin{array}[]{cc}\boldsymbol{H}(\boldsymbol{x})&\boldsymbol{U}^{t}_{{\cal
Z}}\\\ \boldsymbol{U}_{{\cal Z}}&{\bf
0}\end{array}\right)\left(\begin{array}[]{c}\Delta\boldsymbol{x}\\\
\boldsymbol{\lambda}_{{\cal
Z}}\end{array}\right)=\left(\begin{array}[]{c}-\nabla
f(\boldsymbol{x})-(\rho+\Delta\rho)\boldsymbol{u}_{\bar{\cal Z}}\\\ {\bf
0}\end{array}\right),$

where the rows of the matrix $\boldsymbol{U}_{{\cal Z}}$ are the constant
differentials, $\boldsymbol{v}_{i}^{t}$, $i\in{\cal Z}_{\text{E}}$, and
$\boldsymbol{w}_{j}^{t}$, $j\in{\cal Z}_{\text{I}}(\boldsymbol{x})$, of the
active constraint functions. Denoting the inverse of matrix as

$\displaystyle\left(\begin{array}[]{cc}\boldsymbol{H}(\boldsymbol{x})&\boldsymbol{U}^{t}_{{\cal
Z}}\\\ \boldsymbol{U}_{{\cal Z}}&{\bf
0}\end{array}\right)^{-1}=\left(\begin{array}[]{cc}\boldsymbol{P}(\boldsymbol{x})&\boldsymbol{Q}(\boldsymbol{x})\\\
\boldsymbol{Q}^{t}(\boldsymbol{x})&\boldsymbol{R}(\boldsymbol{x})\end{array}\right),$

where

$\displaystyle\boldsymbol{P}(\boldsymbol{x})$ $\displaystyle=$
$\displaystyle\boldsymbol{H}^{-1}(\boldsymbol{x})-\boldsymbol{H}^{-1}(\boldsymbol{x})\boldsymbol{U}_{{\cal
Z}}^{t}\left[\boldsymbol{U}_{{\cal
Z}}\boldsymbol{H}^{-1}(\boldsymbol{x})\boldsymbol{U}_{{\cal
Z}}^{t}\right]^{-1}\boldsymbol{U}_{{\cal
Z}}\boldsymbol{H}^{-1}(\boldsymbol{x})$
$\displaystyle\boldsymbol{Q}(\boldsymbol{x})$ $\displaystyle=$
$\displaystyle\boldsymbol{H}^{-1}(\boldsymbol{x})\boldsymbol{U}_{{\cal
Z}}^{t}\left[\boldsymbol{U}_{{\cal
Z}}\boldsymbol{H}^{-1}(\boldsymbol{x})\boldsymbol{U}_{{\cal
Z}}^{t}\right]^{-1}$ (12) $\displaystyle\boldsymbol{R}(\boldsymbol{x})$
$\displaystyle=$ $\displaystyle-\left[\boldsymbol{U}_{{\cal
Z}}\boldsymbol{H}^{-1}(\boldsymbol{x})\boldsymbol{U}_{{\cal
Z}}^{t}(\boldsymbol{x})\right]^{-1},$

the solution of the difference vector $\Delta\boldsymbol{x}$ is

$\displaystyle\Delta\boldsymbol{x}$ $\displaystyle=$
$\displaystyle-\boldsymbol{P}(\boldsymbol{x})[\nabla
f(\boldsymbol{x})+(\rho+\Delta\rho)\boldsymbol{u}_{\bar{\cal Z}}]$
$\displaystyle=$ $\displaystyle-\boldsymbol{P}(\boldsymbol{x})[\nabla
f(\boldsymbol{x})+\rho\boldsymbol{u}_{\bar{\cal
Z}}(\boldsymbol{x})+\Delta\rho\cdot\boldsymbol{u}_{\bar{\cal Z}}]$
$\displaystyle=$
$\displaystyle-\boldsymbol{P}(\boldsymbol{x})[-\rho\boldsymbol{U}^{t}_{{\cal
Z}}\boldsymbol{r}_{{\cal Z}}+\Delta\rho\cdot\boldsymbol{u}_{\bar{\cal Z}}].$

Note $\boldsymbol{P}(\boldsymbol{x})\boldsymbol{U}^{t}_{{\cal Z}}={\bf 0}$.
Therefore
$\Delta\boldsymbol{x}=-\Delta\rho\cdot\boldsymbol{P}(\boldsymbol{x})\boldsymbol{u}_{\bar{\cal
Z}}$. This gives the direction for the infinitesimal update of solution vector
$\boldsymbol{x}(\rho)$. Taking limit in $\Delta\rho$ leads to the following
key result for developing the path algorithm.

######  Proposition 3.1.

Within interior of a path segment with set configuration (7), the solution
$\boldsymbol{x}(\rho)$ satisfies an ordinary differential equation (ODE)

$\displaystyle\frac{d\boldsymbol{x}(\rho)}{d\rho}=-\boldsymbol{P}(\boldsymbol{x})\boldsymbol{u}_{\bar{\cal
Z}}$ (13)

where the matrix $\boldsymbol{P}(\boldsymbol{x})$ and vector
$\boldsymbol{u}_{\bar{\cal Z}}$ are defined by (12) and (9).

Note that the right hand side of (13) is a constant vector in $\boldsymbol{x}$
when $f$ is quadratic and $g_{i}$ and $h_{j}$ are affine. Thus the
corresponding solution path is piecewise linear. This recovers the case
studied in (Zhou and Lange, 2011). The differential equation (13) holds on the
current segment until one of two types of events happens: an inactive
constraint becomes active or vice versa. The first type of event is easy to
detect – whenever a constraint function, $g_{i}(\boldsymbol{x})$, $i\in{\cal
N}_{\text{E}}\cup{\cal P}_{\text{E}}$, or $h_{j}(\boldsymbol{x})$, $j\in{\cal
N}_{\text{I}}\cup{\cal P}_{\text{I}}$, hits zero, we move that constraint to
the active set ${\cal Z}_{\text{E}}$ or ${\cal Z}_{\text{I}}$ and start
solving a new system of differential equations. To detect when the second type
of event happens, we need to keep track of the coefficients
$s_{i}(\boldsymbol{x})$ and $t_{j}(\boldsymbol{x})$ for active constraints.
Whenever the coefficient of an active constraint hits the boundary of its
permissible range in (6), the constraint has to be relaxed from being active
in next segment. It turns out the coefficients for active constraints admit a
simple representation in terms of current solution vector.

######  Proposition 3.2.

On a path segment with set configuration (7), the coefficients $s_{i}$ and
$t_{j}$ for active constraints are

$\displaystyle\boldsymbol{r}_{{\cal
Z}}(\rho)=\left(\begin{array}[]{c}\boldsymbol{s}_{{\cal
Z}_{\text{E}}}(\rho)\\\ \boldsymbol{t}_{{\cal
Z}_{\text{I}}}(\rho)\end{array}\right)=-\boldsymbol{Q}^{t}(\boldsymbol{x})\left[\frac{1}{\rho}\nabla
f(\boldsymbol{x})+\boldsymbol{u}_{\bar{\cal Z}}\right]$ (16)

where $\boldsymbol{x}=\boldsymbol{x}(\rho)$ is the solution at $\rho$ and the
matrix $\boldsymbol{Q}(\boldsymbol{x})$ is defined by (12).

###### Proof.

Stationarity condition (5) implies

$\displaystyle\boldsymbol{U}_{\cal Z}^{t}\boldsymbol{r}_{\cal
Z}=-\frac{1}{\rho}\nabla f(\boldsymbol{x})-\boldsymbol{U}_{\bar{\cal
Z}}^{t}\boldsymbol{r}_{\bar{\cal Z}}=-\frac{1}{\rho}\nabla
f(\boldsymbol{x})-\boldsymbol{u}_{\bar{\cal Z}}.$

Multiplying both sides by $\boldsymbol{Q}(\boldsymbol{x})$ gives (16). ∎

Given current solution vector $\boldsymbol{x}(\rho)$, the coefficients of the
active constraints are readily obtained from (16). Once a coefficient hits the
end points, we move that constraint from the active set to the inactive set
that matches the end point being hit. In next section, we detail the
implementation of the path algorithm.

##  4 Implementation: ODE and Sweeping Operator

Initialize $\rho=0$,
$\boldsymbol{\beta}(0)=\text{argmin}f(\boldsymbol{\beta})$ and its set
configuration (7).

repeat

Solve ODE (13) until an inactive constraint becomes active or the coefficient
(16) of an active constraint hits boundary.

Update the set configuration (7).

until ${\cal N}_{\text{E}}={\cal P}_{\text{E}}={\cal P}_{\text{I}}=\emptyset$

Algorithm 1 EPSODE: Solution path for regularization problem (1) with strictly
convex $f$.

Algorithm 1 summarizes EPSODE based on Propositions 3.1 and 3.2. It involves
solving ODEs segment by segment and is extremely simple to implement using
softwares with a reliable ODE solver such as the ode45 function in Matlab and
the deSolve package (Karline10RdeSolve) in R. There has been extensive
research in applied mathematics on numerical methods for solving ODEs, notably
the Runge-Kutta, Richardson extrapolation and predictor-corrector methods.
Some path following algorithms developed for specific statistical problems
(Park and Hastie, 2007; Friedman, 2008) turn out to be approximate methods for
solving the corresponding ODE. Wu (2011) first explicitly uses ODE to derive
an exact solution path for the lasso penalized GLM. The connection of path
following to ODE relieves statisticians from the burden of developing specific
path algorithms for a variety of regularization problems.

Any ODE solver repeatedly evaluates the derivative. Suppose the number of
parameters is $p$. Computation of the matrix-vector multiplications in (13)
and (16) has computation cost of order $O(p^{2})+O(p|{\cal Z}|)+O(|{\cal
Z}|^{3})$ if the inverse $H^{-1}$ of Hessian matrix of loss function $f$ is
readily available, where ${\cal Z}={\cal Z}_{\text{E}}\cup{\cal Z}_{\text{I}}$
and $|{\cal Z}|$ denotes its cardinality. Otherwise the computation cost is
$O(p^{3})+O(p|{\cal Z}|)+O(|{\cal Z}|^{3})$.

An alternative implementation avoids repeated matrix inversions by solving an
ODE for the matrices $\boldsymbol{P}$, $\boldsymbol{Q}$ and $\boldsymbol{R}$
themselves. The computations can be conveniently organized around the
classical sweep and inverse sweep operators of regression analysis (Dempster,
1969; Goodnight, 1979; Jennrich, 1977; Little and Rubin, 2002; Lange, 2010).
Suppose $\boldsymbol{A}$ is an $m\times m$ symmetric matrix. Sweeping on the
$k$th diagonal entry $a_{kk}\neq 0$ of $\boldsymbol{A}$ yields a new symmetric
matrix $\widehat{\boldsymbol{A}}$ with entries

$\displaystyle\hat{a}_{kk}$ $\displaystyle=$ $\displaystyle-\frac{1}{a_{kk}},$
$\displaystyle\hat{a}_{ik}$ $\displaystyle=$
$\displaystyle\frac{a_{ik}}{a_{kk}},\quad i\neq k$ $\displaystyle\hat{a}_{kj}$
$\displaystyle=$ $\displaystyle\frac{a_{kj}}{a_{kk}},\quad j\neq k$
$\displaystyle\hat{a}_{ij}$ $\displaystyle=$ $\displaystyle
a_{ij}-\frac{a_{ik}a_{kj}}{a_{kk}},\quad i,j\neq k.$

These arithmetic operations can be undone by inverse sweeping on the same
diagonal entry. Inverse sweeping on the $k$th diagonal entry sends the
symmetric matrix $\boldsymbol{A}$ into the symmetric matrix
$\check{\boldsymbol{A}}$ with entries

$\displaystyle\check{a}_{kk}$ $\displaystyle=$
$\displaystyle-\frac{1}{a_{kk}},$ $\displaystyle\check{a}_{ik}$
$\displaystyle=$ $\displaystyle-\frac{a_{ik}}{a_{kk}},\quad i\neq k$
$\displaystyle\check{a}_{kj}$ $\displaystyle=$
$\displaystyle-\frac{a_{kj}}{a_{kk}},\quad j\neq k$
$\displaystyle\check{a}_{ij}$ $\displaystyle=$ $\displaystyle
a_{ij}-\frac{a_{ik}a_{kj}}{a_{kk}},\quad i,j\neq k.$

Both sweeping and inverse sweeping preserve symmetry. Thus, all operations can
be carried out on either the lower or upper triangle of $\boldsymbol{A}$
alone, saving both computational time and storage. When several sweeps or
inverse sweeps are performed, their order is irrelevant.

At beginning ($\rho=0$) of the path following, we initialize a sweeping
tableau as

$\displaystyle\left(\begin{array}[]{c|c}\boldsymbol{H}^{-1}(\boldsymbol{x})&\boldsymbol{H}^{-1}(\boldsymbol{x})\boldsymbol{U}^{t}\\\
\hline\cr*&\boldsymbol{U}\boldsymbol{H}^{-1}(\boldsymbol{x})\boldsymbol{U}^{t}\end{array}\right),$

where the matrix $\boldsymbol{U}\in\mathbb{R}^{(r+s)\times p}$ holds all
constraint differentials $\boldsymbol{v}_{i}^{t}$ and $\boldsymbol{w}_{j}^{t}$
in rows. Further sweeping of diagonal entries corresponding to the active
constraints yields

$\displaystyle\left(\begin{array}[]{c|cc}\boldsymbol{P}(\boldsymbol{x})&\boldsymbol{Q}(\boldsymbol{x})&\boldsymbol{P}(\boldsymbol{x})\boldsymbol{U}_{\bar{\cal
Z}}^{t}\\\
\hline\cr*&\boldsymbol{R}(\boldsymbol{x})&\boldsymbol{Q}^{t}(\boldsymbol{x})\boldsymbol{U}_{\bar{\cal
Z}}^{t}\\\ &*&\boldsymbol{U}_{\bar{\cal
Z}}\boldsymbol{P}(\boldsymbol{x})\boldsymbol{U}_{\bar{\cal
Z}}^{t}\end{array}\right).$ (21)

Here we conveniently organized the columns of the swept active constraints
before those of un-swept ones. In practice the sweep tableau is not necessary
as in (21) and it is enough to keep an indicator vector recording which
columns are swept. The key elements for the path algorithm magically appear in
the sweep tableau (21)

$\displaystyle\frac{d\boldsymbol{x}(\rho)}{d\rho}$ $\displaystyle=$
$\displaystyle-\boldsymbol{P}(\boldsymbol{x})\boldsymbol{U}_{\bar{\cal
Z}}^{t}\boldsymbol{r}_{\bar{\cal Z}}$ $\displaystyle\boldsymbol{r}_{{\cal
Z}}(\rho)$ $\displaystyle=$
$\displaystyle-\boldsymbol{Q}^{t}(\boldsymbol{x})\boldsymbol{U}_{\bar{\cal
Z}}^{t}\boldsymbol{r}_{\bar{\cal
Z}}-\frac{1}{\rho}\boldsymbol{Q}^{t}(\boldsymbol{x})\nabla f(\boldsymbol{x}).$

Therefore path following procedure only involves solving ODE for the whole
sweep tableau (21) with sweeping or inverse sweeping at kinks between
successive segments. For this purpose we derive the ODE for the sweep tableau
(21). We adopt the convenient notations in (Magnus and Neudecker, 1999). For a
matrix function $F(\boldsymbol{X}):\mathbb{R}^{n\times
q}\to\mathbb{R}^{m\times p}$,

$\displaystyle
DF(\boldsymbol{X})=\frac{\partial\text{vec}F(\boldsymbol{X})}{\partial(\text{vec}\boldsymbol{X})^{t}}$

denotes the $mp\times nq$ Jacobian matrix. For example, Proposition 3.1 states
$D\boldsymbol{x}(\rho)=-\boldsymbol{P}(\boldsymbol{x})\boldsymbol{u}_{\bar{\cal
Z}}$.

######  Proposition 4.1 (ODE for Sweep Tableau).

On a segment of path with fixed set configuration, the matrices
$\boldsymbol{P}(\rho)$, $\boldsymbol{Q}(\rho)$ and $\boldsymbol{R}(\rho)$
satisfy the ordinary differential equations (ODE)

$\displaystyle D\boldsymbol{P}(\rho)$
$\displaystyle=[\boldsymbol{P}(\boldsymbol{x})\otimes\boldsymbol{P}(\boldsymbol{x})]\cdot[D\boldsymbol{H}(\boldsymbol{x})]\cdot\boldsymbol{P}(\boldsymbol{x})\boldsymbol{u}_{\bar{\cal
Z}}$ $\displaystyle D\boldsymbol{Q}(\rho)$
$\displaystyle=[\boldsymbol{Q}^{t}(\boldsymbol{x})\otimes\boldsymbol{P}(\boldsymbol{x})]\cdot[D\boldsymbol{H}(\boldsymbol{x})]\cdot\boldsymbol{P}(\boldsymbol{x})\boldsymbol{u}_{\bar{\cal
Z}}$ $\displaystyle D\boldsymbol{R}(\rho)$
$\displaystyle=[\boldsymbol{Q}^{t}(\boldsymbol{x})\otimes\boldsymbol{Q}^{t}(\boldsymbol{x})]\cdot[D\boldsymbol{H}(\boldsymbol{x})]\cdot\boldsymbol{P}(\boldsymbol{x})\boldsymbol{u}_{\bar{\cal
Z}}.$

###### Proof.

First consider

$\displaystyle\boldsymbol{R}(\boldsymbol{x})=-\left[\boldsymbol{U}_{{\cal
Z}}\boldsymbol{H}^{-1}(\boldsymbol{x})\boldsymbol{U}_{{\cal
Z}}^{t}\right]^{-1}=-\boldsymbol{M}^{-1}(\boldsymbol{x}).$

By chain rule (Magnus and Neudecker, 1999, p91),

$\displaystyle D\boldsymbol{R}(\rho)$ $\displaystyle=$ $\displaystyle
D\boldsymbol{R}(\boldsymbol{M})\cdot D\boldsymbol{M}(\boldsymbol{H})\cdot
D\boldsymbol{H}(\boldsymbol{x})\cdot D\boldsymbol{x}(\rho)$ $\displaystyle=$
$\displaystyle[\boldsymbol{R}(\boldsymbol{x})\otimes\boldsymbol{R}(\boldsymbol{x})]\cdot
D\boldsymbol{M}(\boldsymbol{H})\cdot D\boldsymbol{H}(\boldsymbol{x})\cdot
D\boldsymbol{x}(\rho)$ $\displaystyle=$
$\displaystyle-[\boldsymbol{R}(\boldsymbol{x})\otimes\boldsymbol{R}(\boldsymbol{x})]\cdot\\{[\boldsymbol{U}_{\cal
Z}\boldsymbol{H}^{-1}(\boldsymbol{x})\otimes\boldsymbol{U}_{\cal
Z}\boldsymbol{H}^{-1}(\boldsymbol{x})]\cdot
D\boldsymbol{H}(\boldsymbol{x})\\}\cdot D\boldsymbol{x}(\rho)$
$\displaystyle=$
$\displaystyle[\boldsymbol{Q}^{t}(\boldsymbol{x})\otimes\boldsymbol{Q}^{t}(\boldsymbol{x})]\cdot[D\boldsymbol{H}(\boldsymbol{x})]\cdot\boldsymbol{P}(\boldsymbol{x})\boldsymbol{u}_{\bar{\cal
Z}}.$

Similar calculations yield formula for
$\boldsymbol{Q}(\boldsymbol{x})=-\boldsymbol{H}^{-1}(\boldsymbol{x})\boldsymbol{U}_{{\cal
Z}}^{t}\boldsymbol{R}(\boldsymbol{x})$ and
$\boldsymbol{P}(\boldsymbol{x})=\boldsymbol{H}^{-1}(\boldsymbol{x})-\boldsymbol{Q}(\boldsymbol{x})\boldsymbol{U}_{\cal
Z}\boldsymbol{H}^{-1}(\boldsymbol{x})$. ∎

Solving ODE for these matrices requires the $p^{2}$-by-$p$ Jacobian matrix of
the Hessian matrix $\boldsymbol{H}(\boldsymbol{x})=d^{2}f(\boldsymbol{x})$,

$\displaystyle
D\boldsymbol{H}(\boldsymbol{x})=\frac{\partial[\text{vec}\boldsymbol{H}(\boldsymbol{x})]}{\partial\text{vec}(\boldsymbol{x})^{t}}=\frac{\partial\text{vec}[df^{2}(\boldsymbol{x})]}{\partial\text{vec}(\boldsymbol{x})^{t}},$

which we provide for each example in Section 7 for convenience. When the
number of parameter $p$ is large, $D\boldsymbol{H}$ is a large matrix. However
there is no need to compute and store $D\boldsymbol{H}$ and we are only
required to compute the matrix vector multiplication
$D\boldsymbol{H}\cdot\boldsymbol{v}$ for any vector $\boldsymbol{v}$. In light
of the useful identity
$(\boldsymbol{B}^{t}\otimes\boldsymbol{A})\text{vec}(\boldsymbol{C})=\text{vec}(\boldsymbol{A}\boldsymbol{C}\boldsymbol{B})$,
evaluating the derivative for the whole tableau only involves multiplying
three matrices and incurs computational cost $O(p^{3})+O(p^{2}|{\cal
Z}|)+O(p|{\cal Z}|^{2})$.

Although we have presented the path algorithm as moving from $\rho=0$ to large
$\rho$, it can be applied in either direction. Lasso and fused-lasso usually
start from the constrained solution, while in presence of general equality
constraints, e.g., polynomial trend filtering, and/or inequality constraints,
the constrained solution is not readily available and the path algorithm must
be initiated at $\rho=0$.

##  5 Extension of EPSODE

So far we have assumed strictly convexity of the loss function $f$. This
unfortunately excludes many interesting applications, especially $p>n$ case of
the regression problems. In this section we briefly indicate an extension of
EPSODE to the case $f$ is convex but not necessarily strictly convex. In the
proof of Proposition 3.1, the infinitesimal change of solution
$\Delta\boldsymbol{x}$ is derived via minimizing the equality-constrained
quadratic program (3), the solution to which requires inverse of Hessian
$\boldsymbol{H}^{-1}$ and thus strict convexity of $f$. Alternatively we may
solve (3) via reparameterization. Let $\boldsymbol{U}_{\cal Z}$ hold the
active constraint vectors and $\boldsymbol{Y}\in\mathbb{R}^{p\times(p-|{\cal
Z}|)}$ be a null space matrix of $\boldsymbol{U}_{\cal Z}$, i.e., the columns
of $\boldsymbol{Y}$ are orthogonal to the rows of $\boldsymbol{U}_{\cal Z}$.
Then the infinitesimal change can be represented as
$\Delta\boldsymbol{x}=\boldsymbol{Y}\Delta\boldsymbol{y}$ for some vector
$\Delta\boldsymbol{y}\in\mathbb{R}^{p-|{\cal Z}|}$. Under this
reparameterization, the quadratic program (3) is equivalent to

$\displaystyle\min_{\Delta\boldsymbol{y}}\,\frac{1}{2}\Delta\boldsymbol{y}^{t}[\boldsymbol{Y}^{t}\boldsymbol{H}(\boldsymbol{x})\boldsymbol{Y}]\Delta\boldsymbol{y}+[df(\boldsymbol{x})+(\rho+\Delta\rho)\boldsymbol{u}_{\bar{\cal
Z}}^{t}]\boldsymbol{Y}\cdot\Delta\boldsymbol{y}$

with explicit solution

$\displaystyle\Delta\boldsymbol{y}=-[\boldsymbol{Y}^{t}\boldsymbol{H}(\boldsymbol{x})\boldsymbol{Y}]^{-1}\boldsymbol{Y}^{t}[\nabla
f(\boldsymbol{x})+(\rho+\Delta\rho)\boldsymbol{u}_{\bar{\cal Z}}].$

Hence the infinitesimal change in $\boldsymbol{x}(\rho)$ is

$\displaystyle\Delta\boldsymbol{x}$
$\displaystyle=-\boldsymbol{Y}[\boldsymbol{Y}^{t}\boldsymbol{H}(\boldsymbol{x})\boldsymbol{Y}]^{-1}\boldsymbol{Y}^{t}[\nabla
f(\boldsymbol{x})+(\rho+\Delta\rho)\boldsymbol{u}_{\bar{\cal Z}}]$
$\displaystyle=-\Delta\rho\cdot\boldsymbol{Y}[\boldsymbol{Y}^{t}\boldsymbol{H}(\boldsymbol{x})\boldsymbol{Y}]^{-1}\boldsymbol{Y}^{t}\boldsymbol{u}_{\bar{\cal
Z}}].$

Again taking limit gives the following result in parallel to Proposition 3.1.

######  Proposition 5.1.

Within interior of a path segment with set configuration (7), the solution
$\boldsymbol{x}(\rho)$ satisfies an ordinary differential equation (ODE)

$\displaystyle\frac{d\boldsymbol{x}(\rho)}{d\rho}=-\boldsymbol{Y}[\boldsymbol{Y}^{t}\boldsymbol{H}(\boldsymbol{x})\boldsymbol{Y}]^{-1}\boldsymbol{Y}^{t}\boldsymbol{u}_{\bar{\cal
Z}}$ (22)

where $\boldsymbol{Y}$ is a null space matrix of $\boldsymbol{U}_{\cal Z}$.

An advantage of (22) is that only non-singularity of the matrix
$\boldsymbol{Y}^{t}\boldsymbol{H}(\boldsymbol{x})\boldsymbol{Y}$ is required
which is much weaker than the non-singularity of $\boldsymbol{H}$. The
computational cost of calculating the derivative in (22) is $O((p-|{\cal
Z}|)^{3})+O(p(p-|{\cal Z}|))$, which is more efficient than (13) when
$p-|{\cal Z}|$ is small. However it requires the null space matrix
$\boldsymbol{Y}$, which is nonunique and may be expensive to compute.
Fortunately the null space matrix $\boldsymbol{Y}$ is constant over each path
segment and in practice can be calculated by QR decomposition of the active
constraint matrix $\boldsymbol{U}_{\cal Z}$. At each kink either one
constraint leaves ${\cal Z}$ or one enters ${\cal Z}$. Therefore
$\boldsymbol{Y}$ can be sequentially updated (Lawson and Hanson, 1987) and
need not to be calculated anew for each segment. Which version of (13) and
(22) to use depends on specific application. When the loss function $f$ is not
strictly convex, e.g., $p>n$ case in regression analysis, only (22) applies.
Interested readers are referred to (Nocedal and Wright, 2006) for a similar
dilemma in optimization methods.

##  6 Model Selection Along the Path

In applications such as penalized GLMs, the tuning parameter $\rho$ in the
regularization problem (1) is chosen by a model selection criterion such as
AIC, BIC, $C_{p}$, or cross-validation. The cross validation errors can be
readily computed using the solution path output by EPSODE. Yet the AIC, BIC,
and $C_{p}$ criteria require an estimate of the degrees of freedom of estimate
$\boldsymbol{\beta}(\rho)$. Specifically AIC and BIC are defined by

AIC
$\displaystyle=-\ell(\boldsymbol{\beta}(\rho))+\text{df}(\boldsymbol{\beta}(\rho))$
BIC $\displaystyle=-\ell(\boldsymbol{\beta}(\rho))+\frac{\log
n}{2}\text{df}(\boldsymbol{\beta}(\rho)),$

where $-\ell(\cdot)$ denotes the negative log-likelihood and
$\text{df}(\boldsymbol{\beta}_{\rho})$ is the degrees of freedom for estimate
$\boldsymbol{\beta}_{\rho}$. We propose to use

$\displaystyle\text{df}(\boldsymbol{\beta}(\rho))=p-|{\cal
Z}_{\text{E}}\cup{\cal Z}_{\text{I}}|$ (23)

as a measure of the degrees of freedom under GLMs. It is previously shown that
(23) is an unbiased estimate of the degrees of freedom for lasso penalized
least squares (Efron et al., 2004; Zou et al., 2007), generalized lasso
penalized least squares (Tibshirani and Taylor, 2011), and the least squares
version of the regularized problem (1) (Zhou and Lange, 2011). Using the same
degrees of freedom formula (23) for GLMs is justified by the local
approximation of GLM loglikelihood by weighted least squares. See (Park and
Hastie, 2007) for details.

##  7 Applications

In this section, we collect some representative regularized or constrained
estimation problems and demonstrate how they can be solved by path following.
For all applications, we list the first three derivatives of the loss function
$f$ in (1). In fact, the third derivative is only needed when implementing by
solving the ODE for the sweep tableau.

###  7.1 GLMs and Quasi-Likelihoods with Generalized $l_{1}$ Regularizations

The generalized linear model (GLM) deals with exponential families in which
the sufficient statistics is $Y$ and the conditional mean $\mu$ of $Y$
completely determines its distribution. Conditional on the covariate vector
$\boldsymbol{x}\in\mathbb{R}^{p}$, the response variable $y$ is modeled as

$\displaystyle
p(y|\boldsymbol{x};\boldsymbol{\beta},\sigma)\propto\exp\left\\{\frac{y\langle\boldsymbol{x},\boldsymbol{\beta}\rangle-\psi(\langle\boldsymbol{x},\boldsymbol{\beta}\rangle)}{c(\sigma)}\right\\},$
(24)

where the scalar $\sigma>0$ is a fixed and known scale parameter and the
vector $\boldsymbol{\beta}$ is the parameters to be estimated. The function
$\psi:\mathbb{R}\mapsto\mathbb{R}$ is the link function. When
$y\in\mathbb{R}$, $\psi(u)=u^{2}/2$ and $c(\sigma)=\sigma^{2}$, (24) is the
normal regression model. When $y\in\\{0,1\\}$, $\psi(u)=\ln(1+\exp(u))$ and
$c(\sigma)=1$, (24) is the logistic regression model. When $y\in\mathbb{N}$,
$\psi(u)=\exp(u)$, and $c(\sigma)=1$, (24) is the Poisson regression model.

The quasi-likelihoods generalize GLM without assuming a specific distribution
form of $Y$. Instead only a function relation between the conditional means
$\mu_{i}$ and variances $\sigma_{i}^{2}$, $\sigma_{i}^{2}=V(\mu_{i})$ for some
variance function $V$, is needed. Then the integral

$\displaystyle Q(\mu,y)=\int_{y}^{\mu}\frac{y-t}{\sigma^{2}V(t)}\,dt$

behaves like a log-likelihood function under mild conditions and is called the
quasi-likelihood. The quasi-likelihood includes GLMs as special cases with
appropriately chosen variance function $V(\cdot)$. Readers are referred to the
classical text (McCullagh and Nelder, 1983, Table 9.1) for the commonly used
quasi-likelihoods. By slightly abusing our notation, we assume a known link
function between the conditional mean $\mu_{i}$ and linear predictor
$\boldsymbol{x}_{i}^{T}\boldsymbol{\beta}$,
$\mu=\mu(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta})$ and denote
$Q_{i}(\boldsymbol{\beta})=Q(\mu(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta}),y_{i})$.
Then the quasi-likelihood with generalized $l_{1}$ regularization takes the
form

$\displaystyle-Q(\boldsymbol{\beta})+\rho\|\boldsymbol{V}\boldsymbol{\beta}-\boldsymbol{d}\|_{1}=-\sum_{i=1}^{n}Q_{i}(\boldsymbol{\beta})+\rho\|\boldsymbol{V}\boldsymbol{\beta}-\boldsymbol{d}\|_{1},$
(25)

which is a special case of the general form (1). Specific choices of the
regularization matrix $\boldsymbol{V}$ and constant vector $\boldsymbol{d}$
lead to lasso, fused-lasso, trend filtering, and many other applications.

For the path algorithm, we require the first two or three derivatives of the
complete quasi-likelihood. Denoting
$\boldsymbol{\eta}=\boldsymbol{X}\boldsymbol{\beta}$ with
$\boldsymbol{X}=(\boldsymbol{x}_{1}^{t},\boldsymbol{x}_{2}^{t},\cdots,\boldsymbol{x}_{n}^{t})^{t}$,
we have

$\displaystyle\nabla Q(\boldsymbol{\beta})$ $\displaystyle=$
$\displaystyle[D\boldsymbol{\mu}(\boldsymbol{\eta})]^{t}\boldsymbol{V}^{-1}(\boldsymbol{y}-\boldsymbol{\mu})/\sigma^{2}=\boldsymbol{X}^{t}[D\boldsymbol{\mu}(\boldsymbol{\eta})]\boldsymbol{V}^{-1}(\boldsymbol{y}-\boldsymbol{\mu})/\sigma^{2},$
$\displaystyle\boldsymbol{H}(\boldsymbol{\beta})=d^{2}Q(\boldsymbol{\beta})$
$\displaystyle=$
$\displaystyle[(\boldsymbol{y}-\boldsymbol{\mu})^{t}\boldsymbol{V}^{-1}\otimes\boldsymbol{X}^{t}]\cdot
D^{2}\boldsymbol{\mu}(\boldsymbol{\eta})\cdot\boldsymbol{X}/\sigma^{2},$ (26)
$\displaystyle D\boldsymbol{H}(\boldsymbol{\beta})=d^{3}Q(\boldsymbol{\beta})$
$\displaystyle=$
$\displaystyle[\boldsymbol{X}^{t}\otimes(\boldsymbol{y}-\boldsymbol{\mu})^{t}\boldsymbol{V}^{-1}\otimes\boldsymbol{X}^{t}]\cdot
D^{3}\boldsymbol{\mu}(\boldsymbol{\eta})\cdot\boldsymbol{X}/\sigma^{2},$

where $\boldsymbol{V}$ is a $n$-by-$n$ diagonal matrix with diagonal entries
$V(\mu(\boldsymbol{x}_{i}^{t}\boldsymbol{y}))$,
$D\boldsymbol{\mu}(\boldsymbol{\eta})$ is a $n$-by-$n$ diagonal matrix with
diagonal entries $\mu^{\prime}(\boldsymbol{x}_{i}^{t}\boldsymbol{\beta})$,
$D^{2}\boldsymbol{\mu}(\boldsymbol{\eta})$ is a $n^{2}$-by-$n$ matrix with
$(n(i-1)+i,i)$ entry equal to
$\mu^{\prime\prime}(\boldsymbol{x}_{i}^{t}\boldsymbol{\beta})$ for
$i=1,\ldots,n$ and 0 otherwise, and $D^{3}\boldsymbol{\mu}(\boldsymbol{\eta})$
is a $n^{3}$-by-$n$ matrix with $(n^{2}(i-1)+n(i-1)+i,i)$ entry equal to
$\mu^{\prime\prime\prime}(\boldsymbol{x}_{i}^{t}\boldsymbol{\beta})$ for
$i=1,\ldots,n$ and 0 otherwise. Note for GLM with canonical link, these
formulas simplify (Agresti, 2002, Section 4.6.4).

The most widely used $l_{1}$ regularization is the lasso penalty which imposes
sparsity on the regression coefficients. For numerical demonstration, we
revisit the M&A example introduced in Section 1 without discretizing each
predictor. We standardize each predictor first and consider the lasso
penalized linear logistic regression model. Figure 3 shows the lasso solution
path for each standardized predictor in the left panel and corresponding AIC
and BIC scores in the right panel. The order at which predictors enter the
model matches the more detailed patterns revealed by the varying coefficient
model in Figure 1. The almost monotone effects of the predictors ‘market-to-
book ratio’, ‘cash flow’, ’cash’, and ’tax’ can be captured by the usual
linear logistic regression and these covariates are picked up by lasso first.
The nonlinear effects shown in the other predictors are likely to be missed by
the linear logistic regression. For instance, the quadratic effects of ‘log
market equity’ shown in the regularized estimates in Figure 1 are missed by
both AIC and BIC criteria.

$\begin{array}[]{cc}\includegraphics[width=180.67499pt]{MandA_lassopath}&\includegraphics[width=180.67499pt]{MandA_lassoAICBIC}\end{array}$
Figure 3: M&A example revisited. Lasso solution path on the seven standardized
predictors.

All the generalized lasso problems studied in (Tibshirani and Taylor, 2011)
for Gaussian linear regression naturally generalize to GLMs or quasi-
likelihoods and are subject to the EPSODE path algorithm. This leads to
applications to lasso or fused-lasso penalized GLMs, outlier detections, trend
filtering, and image restoration for GLMs. For instance, cubic trend filtering
is performed on five predictors of the M&A example in Section 1. The graph-
guided penalized linear regression proposed in (Chen et al., 2010) can also be
generalized to GLMs or quasi-likelihoods. Suppose each node $i$ of a graph is
assigned a regression coefficient $\beta_{i}$. In graph penalized regression,
the objective function takes the form

$\displaystyle-\ell(\boldsymbol{\beta})+\lambda_{\text{G}}\sum_{i\sim
j}\left|\frac{\beta_{i}}{\sqrt{d_{i}}}-\text{sgn}(r_{ij})\frac{\beta_{j}}{\sqrt{d_{j}}}\right|+\lambda_{\text{L}}\sum_{j}|\beta_{j}|,$
(27)

where the set of neighboring pairs $i\sim j$ define the graph, $d_{i}$ is the
degree of node $i$, and $r_{ij}$ is the correlation coefficient between $i$
and $j$. This is simply a special case of (25) when the ratio
$\lambda_{\text{G}}/\lambda_{\text{L}}$ is fixed.

###  7.2 Shape-Restricted Regressions

Order-constrained regression has been an important modeling tool (Robertson et
al., 1988; Silvapulle and Sen, 2005). If $\boldsymbol{\beta}$ denotes the
parameter vector, monotone regression imposes isotone constraints
$\beta_{1}\leq\beta_{2}\leq\cdots\leq\beta_{p}$ or antitone constraints
$\beta_{1}\geq\beta_{2}\geq\cdots\geq\beta_{p}$. In partially ordered
regression, subsets of the parameters are subject to isotone or antitone
constraints. In some other problems it is sensible to impose convex or concave
constraints. Note that if locations of regression parameters are at
irregularly spaced time points $t_{1}\leq t_{2}\leq\cdots\leq t_{p}$,
convexity translates into the constraints

$\displaystyle\frac{\beta_{i+2}-\beta_{i+1}}{t_{i+2}-t_{i+1}}\geq\frac{\beta_{i+1}-\beta_{i}}{t_{i+1}-t_{i}}$

for $1\leq i\leq p-2$. When the time intervals are uniform, the constraints
simplify to $\beta_{i+2}-\beta_{i+1}\geq\beta_{i+1}-\beta_{i}$,
$i=1,2,\cdots,p-1$. Concavity translates into the opposite set of
inequalities.

Most of previous work has focused on the linear regression problems because of
the computational and theoretical complexities in the generalized linear model
setting. The recent work (Rufibach, 2010) proposes an active set algorithm for
GLMs with order constraints. The EPSODE algorithm conveniently provides a
solution to the linearly constrained estimation problem (1.1). The relevant
derivatives of loss function are listed in (26). It is noteworthy that EPSODE
not only provides the constrained estimate but also the whole path bridging
the unconstrained estimate to the constrained solution. Availability of the
whole solution path renders model selection between the two extremes simple.

In the illustrative M&A example of Section 1, the bin predictors for the
‘market-to-book ratio’ are regularized by the antitone constraint and those
for the ‘log market equity’ covariate by the concavity constraint.

###  7.3 Gaussian Graphical Models

In recent years several authors (Friedman et al., 2008; Yuan, 2008) proposed
to estimate the sparse undirected graphical model by using lasso
regularizations to the log-likelihood function of the precision matrix, the
inverse of the variance-covariance matrix. Given an observed variance-
covariance matrix $\hat{\boldsymbol{\Sigma}}\in R^{p\times p}$, the negative
log-likelihood of the precision matrix
$\boldsymbol{\Omega}=\boldsymbol{\Sigma}^{-1}$ under normal assumption is

$\displaystyle
f(\boldsymbol{\Omega})=-\log\det\boldsymbol{\Omega}+\text{tr}(\hat{\boldsymbol{\Sigma}}\boldsymbol{\Omega})$
(28)

with the MLE solution $\hat{\boldsymbol{\Sigma}}^{-1}$ when
$\hat{\boldsymbol{\Sigma}}$ is non-degenerate. A zero in the precision matrix
implies conditional independence of the corresponding nodes. Graphical lasso
proposes to solve

$\displaystyle f(\boldsymbol{\Omega})+\rho\sum_{i<j}|\omega_{ij}|,$ (29)

where $\rho\geq 0$ is the tuning constant and $\omega_{ij}$ denotes the
$(i,j)$-element of $\boldsymbol{\Omega}$. It is well-known that the
determinant function is log-concave (Magnus and Neudecker, 1999). Therefore
the loss function $f$ (28) is convex and the EPSODE algorithm applies to (29).
Friedman et al. (2008) proposed an efficient coordinate descent procedure for
solving (29) at a fixed $\rho$. A recent attempt to approximate the whole
solution path is made by Yuan (2008). Again his path algorithm can be deemed
as a primitive predictor-corrector method for approximating the ODE solution.

With symmetry in mind, we parameterize $\boldsymbol{\Omega}$ in terms of its
lower triangular part by a $p(p+1)/2$ column vector $\boldsymbol{x}$ and let
$D\boldsymbol{\Omega}(\boldsymbol{x})=\frac{\partial\text{vec}\boldsymbol{\Omega}}{\partial(\text{vec}\boldsymbol{x})^{t}}$
be the corresponding $p^{2}$-by-$p(p+1)/2$ Jacobian matrix. Note
$D\boldsymbol{\Omega}(\boldsymbol{x})\cdot\boldsymbol{x}=\text{vec}\boldsymbol{\Omega}(\boldsymbol{x})$
and each row of $D\boldsymbol{\Omega}(\boldsymbol{x})$ has exactly one nonzero
entry which equals unity. We list here the first three derivatives of $f$. The
proof is straightforward using matrix calculus and omitted for brevity.

######  Lemma 7.1.

  1. 1.

The derivatives for the Gaussian graphical model (28) with respect to
$\boldsymbol{\Omega}$ are

$\displaystyle Df(\boldsymbol{\Omega})$
$\displaystyle=df(\boldsymbol{\Omega})=[\text{vec}(-\boldsymbol{\Omega}^{-1}+\boldsymbol{\Sigma})]^{t}$
$\displaystyle D^{2}f(\boldsymbol{\Omega})$
$\displaystyle=d^{2}f(\boldsymbol{\Omega})=\boldsymbol{\Omega}^{-1}\otimes\boldsymbol{\Omega}^{-1}$
$\displaystyle D^{3}f(\boldsymbol{\Omega})$
$\displaystyle=-(\boldsymbol{I}_{n}\otimes\boldsymbol{K}_{nn}\otimes\boldsymbol{I}_{n})$
$\displaystyle\hskip
14.45377pt\cdot[\boldsymbol{\Omega}^{-1}\otimes\boldsymbol{\Omega}^{-1}\otimes\text{vec}(\boldsymbol{\Omega}^{-1})+\text{vec}(\boldsymbol{\Omega}^{-1})\otimes\boldsymbol{\Omega}^{-1}\otimes\boldsymbol{\Omega}^{-1}],$

where $\boldsymbol{K}_{nn}$ is the commutation matrix (Magnus and Neudecker,
1999).

  2. 2.

The derivatives for the Gaussian graphical model (28) with respect to
$\boldsymbol{x}$ are

$\displaystyle Df(\boldsymbol{x})$ $\displaystyle=Df(\boldsymbol{\Omega})\cdot
D\boldsymbol{\Omega}(\boldsymbol{x})$
$\displaystyle\boldsymbol{H}(\boldsymbol{x})=D^{2}f(\boldsymbol{x})$
$\displaystyle=[D\boldsymbol{\Omega}(\boldsymbol{x})]^{t}\cdot
D^{2}f(\boldsymbol{\Omega})\cdot D\boldsymbol{\Omega}(\boldsymbol{x})$
$\displaystyle D\boldsymbol{H}(\boldsymbol{x})=D^{3}f(\boldsymbol{x})$
$\displaystyle=\\{[D\boldsymbol{\Omega}(\boldsymbol{x})]^{t}\otimes[D\boldsymbol{\Omega}(\boldsymbol{x})]^{t}\\}\cdot
D^{3}f(\boldsymbol{\Omega})\cdot D\boldsymbol{\Omega}(\boldsymbol{x}).$

When the covariance matrix $\hat{\boldsymbol{\Sigma}}$ is nonsingular, EPSODE
can be initiated either at $\rho=0$ or $\rho=\infty$. When
$\hat{\boldsymbol{\Sigma}}$ is singular, we start from $\rho=\infty$ and the
extended version of EPSODE (22) should be used. If starting at $\rho=0$, the
solution is initialized at $\hat{\boldsymbol{\Sigma}}^{-1}$; If starting at
$\rho=\infty$, the solution is initialized at
$\text{diag}(\hat{\sigma}_{ii}^{-1})$. Minimization of both the unpenalized
and penalized objective function has to be performed over the convex cone of
symmetric, positive semidefinite matrices, which is not explicitly
incorporated in our path following algorithm. The next result ensures the
positive definiteness of the path solution.

######  Lemma 7.2 (Positive definiteness along the path).

The path solution $\boldsymbol{\Omega}(\rho)$ minimizes (29) over the convex
cone of symmetric, positive semidefinite matrices.

###### Proof.

Both symmetry and the stationarity condition (5) are preserved by the path
following. Therefore the path solution $\boldsymbol{\Omega}(\rho)$ constitutes
a minimum of the penalized objective function (29). This implies that
$\text{det}(\boldsymbol{\Omega}(\rho))>0$ otherwise
$f(\boldsymbol{\Omega}(\rho))=\infty$, contradicting with the optimality of
$\boldsymbol{\Omega}(\rho)$. ∎

We illustrate the path algorithm by the classical example of 88 students’
scores on five math courses – mechanics, vector, algebra, analysis, and
statistics (Mardia et al., 1979, Table 1.2.1). Figure 4 displays the solution
path from EPSODE. The top three edges chosen by lasso are analysis-algebra,
statistics-algebra, and algebra-vector, matching the findings in (Yuan, 2008).

$\begin{array}[]{c}\includegraphics[width=252.94499pt]{MKBscore_lassopath}\end{array}$
Figure 4: Solution path of the 10 edges in lasso-regularized Gaussian
graphical model for the math score data. The top three edges chosen by lasso
are labeled.

###  7.4 Nonparametric Density Estimation

As part of the trend of statistics shifting from parametric models to semi- or
non-parametric models, nonparametric density estimation has attracted much
attention in recent years. The maximum likelihood estimation for nonparametric
density estimation often involves a nontrivial, high-dimensional constrained
optimization problem. In this section, we briefly demonstrate the
applicability of EPSODE to the maximum likelihood estimation of univariate
log-concave density. Extensions to multivariate log-concave density estimation
(Cule et al., 2009, 2010) and $k$-monotone density estimation (Balabdaoui and
Wellner, 2007) will be pursued elsewhere. Some algorithms have been
specifically crafted for log-concave density estimation, e.g., the iterative
convex minorant algorithm (ICMA) (Groeneboom and Wellner, 1992; Jongbloed,
1998) and more recently an active-set algorithm (Duembgen et al., 2007). It is
noteworthy that, besides providing an alternative solver for log-concave
density estimation, EPSODE offers the whole solution path between the
unconstrained and constrained solutions. For example, an “almost” log-concave
density estimate in the middle of the path can be chosen that minimizes cross-
validation or prediction error. This adds another dimension to the flexibility
of nonparametric modeling.

The family of log-concave densities is an attractive modeling tool. It
includes most of the commonly used parametric distributions as special cases.
Examples include normal, gamma with shape parameter $\geq 1$, and beta
densities with both parameters $\geq 1$. The survey paper (Walther, 2009)
gives a recent review. A probability density $g(\cdot)$ on $\mathbb{R}$ is
log-concave if its logarithm $\phi(x)=\ln g(x)$ is concave. Given iid
observations, from an unknown distribution of density $g(\cdot)$, with support
at points $x_{1}<\ldots<x_{n}$ with corresponding frequencies
$p_{1},\ldots,p_{n}$, it is well-known (Walther, 2002) that the nonparametric
MLE of $g$ exists, is unique and takes the form $\hat{g}=\exp(\hat{\phi})$
where $\hat{\phi}$ is continuous and piecewise linear on $[x_{1},x_{n}]$, with
the set of knots contained in $\\{x_{1},\ldots,x_{n}\\}$, and
$\hat{\phi}=-\infty$ outside the interval $[x_{1},x_{n}]$. This implies that
the MLE is obtained by minimizing the strictly convex function

$\displaystyle
f(\boldsymbol{\phi})=-\sum_{i=1}^{n}p_{i}\phi_{i}+\sum_{k=1}^{n-1}(x_{k+1}-x_{k})\int_{0}^{1}e^{(1-t)\phi_{k}+t\phi_{k+1}}\,dt.$

over
$\boldsymbol{\phi}=(\phi_{1},\phi_{2},\cdots,\phi_{n})^{t}\in\mathbb{R}^{n}$
subject to constraints

$\displaystyle\frac{\phi_{i+1}-\phi_{i}}{x_{i+1}-x_{i}}\leq\frac{\phi_{i}-\phi_{i-1}}{x_{i}-x_{i-1}},\hskip
7.22743pti=2,\ldots,n-1.$

The consistency of the MLE is proved by Pal et al. (2007) and the pointwise
asymptotic distribution of the MLE studied in (Balabdaoui et al., 2009).

Following Duembgen et al. (2007), we use notations

$\displaystyle\delta_{0}$
$\displaystyle=\delta_{n}=0,\,\delta_{i}=x_{i+1}-x_{i},\hskip
7.22743pti=1,\ldots,n-1$ $\displaystyle J(r,s)$
$\displaystyle=\int_{0}^{1}e^{(1-t)r+ts}\,dt=\begin{cases}\frac{e^{s}-e^{r}}{s-r}&r\neq
s\\\ e^{r}&r=s\end{cases}.$

Then the objective function becomes

$\displaystyle
f(\boldsymbol{\phi})=-\sum_{i=1}^{n}p_{i}\phi_{i}+\sum_{k=1}^{n-1}\delta_{k}J(\phi_{k},\phi_{k+1}).$

The path algorithm requires up to the third derivative of the objective
function $f$

$\displaystyle[\nabla f(\boldsymbol{\phi})]_{i}$
$\displaystyle=-p_{i}+\delta_{i-1}J_{01}(\phi_{i-1},\phi_{i})+\delta_{i}J_{10}(\phi_{i},\phi_{i+1})$
$\displaystyle[H(\boldsymbol{\phi})]_{ij}$
$\displaystyle=[d^{2}f(\boldsymbol{\phi})]_{ij}$
$\displaystyle=\begin{cases}\delta_{i-1}J_{11}(\phi_{i-1},\phi_{i})&j=i-1\\\
\delta_{i-1}J_{02}(\phi_{i-1},\phi_{i})+\delta_{i}J_{20}(\phi_{i},\phi_{i+1})&j=i\\\
\delta_{i}J_{11}(\phi_{i},\phi_{i+1})&j=i+1\\\ 0&\text{otherwise}\end{cases}$
$\displaystyle\frac{\partial[H(\boldsymbol{\phi})]_{i,i-1}}{\partial\phi_{k}}$
$\displaystyle=\begin{cases}\delta_{i-1}J_{21}(\phi_{i-1},\phi_{i})&k=i-1\\\
\delta_{i-1}J_{12}(\phi_{i-1},\phi_{i})&k=i\\\ 0&\text{otherwise}\end{cases}$
$\displaystyle\frac{\partial[H(\boldsymbol{\phi})]_{i,i}}{\partial\phi_{k}}$
$\displaystyle=\begin{cases}\delta_{i-1}J_{12}(\phi_{i-1},\phi_{i})&k=i-1\\\
\delta_{i-1}J_{03}(\phi_{i-1},\phi_{i})+\delta_{i}J_{30}(\phi_{i},\phi_{i+1})&k=i\\\
\delta_{i}J_{21}(\phi_{i},\phi_{i+1})&k=i+1\\\ 0&\text{otherwise}\end{cases}$
$\displaystyle\frac{\partial[H(\boldsymbol{\phi})]_{i,i+1}}{\partial\phi_{k}}$
$\displaystyle=\begin{cases}\delta_{i}J_{21}(\phi_{i},\phi_{i+1})&k=i\\\
\delta_{i}J_{12}(\phi_{i},\phi_{i+1})&k=i+1\\\ 0&\text{otherwise}\end{cases}.$

Interchanging the derivative and integral operators, justified by the
dominated convergence theorem, gives a useful representation for the partial
derivatives of $J$

$\displaystyle J_{ab}(r,s)$ $\displaystyle=\frac{\partial^{a+b}}{\partial
r^{a}\partial s^{b}}J(r,s)=\int_{0}^{1}(1-t)^{a}t^{b}e^{(1-t)r+ts}\,dt.$

We derive a recurrence relation for $J_{ab}(r,s)$ to facilitate its
computation.

######  Lemma 7.3.

$J_{ab}(r,s)$ satisfy following recurrence

  1. 1.

For $r\neq s$,

$\displaystyle J_{00}(r,s)$ $\displaystyle=\frac{e^{s}-e^{r}}{s-r}$
$\displaystyle J_{10}(r,s)$
$\displaystyle=-\frac{e^{r}}{s-r}+\frac{e^{s}-e^{r}}{(s-r)^{2}}$
$\displaystyle J_{01}(r,s)$
$\displaystyle=\frac{e^{s}}{s-r}-\frac{e^{s}-e^{r}}{(s-r)^{2}}$ $\displaystyle
J_{11}(r,s)$
$\displaystyle=\frac{e^{s}+e^{r}}{(s-r)^{2}}-\frac{2(e^{s}-e^{r})}{(s-r)^{3}}$
$\displaystyle J_{ab}(r,s)$
$\displaystyle=\frac{a+b+s-r}{s-r}J_{a-1,b}(r,s)-\frac{a-1}{s-r}J_{a-2,b}(r,s)$
$\displaystyle J_{ab}(r,s)$
$\displaystyle=-\frac{a+b-s+r}{s-r}J_{a,b-1}(r,s)+\frac{b-1}{s-r}J_{a,b-2}(r,s).$

  2. 2.

For $r=s$,

$\displaystyle
J_{ab}(r,s)=\frac{e^{r}a!b!}{(a+b+1)!}=\frac{a}{a+b+1}J_{a-1,b}=\frac{b}{a+b+1}J_{a,b-1}.$

###### Proof.

We recognize $J_{ab}(r,s)$ as the Kummer confluent hypergeometric function
multiplied by a constant

$\displaystyle J_{ab}(r,s)$
$\displaystyle=e^{r}B(a+1,b+1)\sum_{k=0}^{\infty}\frac{(b+1)_{(k)}}{(a+b+2)_{(k)}}\frac{(s-r)^{k}}{k!}$
$\displaystyle=e^{r}B(a+1,b+1)\,_{1}F_{1}(b+1,a+b+2\mid s-r)$
$\displaystyle=e^{s}B(a+1,b+1)\,_{1}F_{1}(a+1,a+b+2\mid r-s)$
$\displaystyle=\frac{e^{s}a!b!}{(a+b+1)!}\,_{1}F_{1}(a+1,a+b+2\mid r-s).$

Then the results follow from the well-known recurrence relation for Kummer
hypergeometric function

$\,{}_{1}F_{1}(x,y\mid z)=\frac{(1-y)(y+z-2)}{(x-y+1)z}\,_{1}F_{1}(x,y-1\mid
z)+\frac{(1-y)(2-y)}{(x-y+1)z}\,_{1}F_{1}(x,y-2\mid z).$

and symmetry $J_{ab}(r,s)=J_{ba}(s,r)$. ∎

To illustrate the path algorithm for this problem, we simulate $n=25$ points
from the extremal distribution Gumbel(0,1). Figure 5 displays the constrained
and unconstrained estimates of $\phi_{i}$ and the solution path bridging the
two.

$\begin{array}[]{cc}\includegraphics[width=180.67499pt]{gumbel_estimates}&\includegraphics[width=180.67499pt]{gumbel_solpath}\\\
\includegraphics[width=180.67499pt]{gumbel_cdf}\end{array}$ Figure 5: Log-
concave density estimation. $n=25$ points are generated from Gumbel(0,1)
distribution. Top left: Unconstrained and concavity-constrained estimates
$\phi$. Top right: Solution path. Bottom left: Empirical cdf and the cdf of
MLE density.

##  8 Conclusions

In this article we propose a generic path following algorithm EPSODE that
works for any regularization problems of form (1). The advantages are its
simplicity and generality. Path following only involves solving ODEs segment
by segment and is simple to implement using popular softwares such as R and
Matlab. Besides providing the whole regularization path, it also gives a
solver for linearly constrained optimization problems that frequently arise in
statistics. Our applications to shape-restricted regressions and nonparametric
density estimation are special cases in particular.

At least two extensions deserve further study. Current algorithm requires
sufficient smoothness (twice differentiable) in the loss function. This
precludes certain applications with non-smooth objective function, e.g., the
Huber loss in robust estimation and the loss function in quantile regression.
Generalization of our path algorithm to regularization of these loss functions
requires further research. Another restriction in our formulation is the
linearity in the regularization terms. In sparse regressions, several authors
have proposed nonlinear and non-convex penalties. The bridge regression (Frank
and Friedman, 1993; Fu, 1998) and SCAD penalties (Fan and Li, 2001) fall into
this category. As observed in (Friedman, 2008), when the penalty is not
convex, the solution path may not be continuous and poses difficulty in path
following, which strongly depends on the continuity and smoothness of the
solution path. Fortunately, in these problems, the discontinuities only occur
when new variables enter or leave the model. A promising strategy is to
initialize the starting point of next segment by solving an equality
constrained optimization problem. This again invites further investigation.

## References

  * Agresti (2002) Agresti, A. (2002).  Categorical Data Analysis (Second ed.).  Wiley Series in Probability and Statistics. New York: Wiley-Interscience [John Wiley & Sons].
  * Balabdaoui et al. (2009) Balabdaoui, F., K. Rufibach, and J. A. Wellner (2009).  Limit distribution theory for maximum likelihood estimation of a log-concave density.  Ann. Statist. 37(3), 1299–1331.
  * Balabdaoui and Wellner (2007) Balabdaoui, F. and J. A. Wellner (2007).  Estimation of a $k$-monotone density: limit distribution theory and the spline connection.  Ann. Statist. 35(6), 2536–2564.
  * Chen et al. (2001) Chen, S. S., D. L. Donoho, and M. A. Saunders (2001).  Atomic decomposition by basis pursuit.  SIAM Rev. 43(1), 129–159.
  * Chen et al. (2010) Chen, X., Q. Lin, S. Kim, J. Carbonell, and E. Xing (2010).  An efficient proximal gradient method for general structured sparse learning.  arXiv:1005.4717.
  * Cule et al. (2009) Cule, M., R. B. Gramacy, and R. Samworth (2009, 1).  LogConcDEAD: An R package for maximum likelihood estimation of a multivariate log-concave density.  Journal of Statistical Software 29(2), 1–20.
  * Cule et al. (2010) Cule, M., R. Samworth, and M. Stewart (2010).  Maximum likelihood estimation of a multi-dimensional log-concave density.  Journal of the Royal Statistical Society Series B 72(5), 545–607.
  * Dempster (1969) Dempster, A. P. (1969).  Elements of Continuous Multivariate Analysis.  Addison-Wesley series in behavioral sciences. Reading, MA: Addison-Wesley.
  * Duembgen et al. (2007) Duembgen, L., K. Rufibach, and A. Huesler (2007).  Active set and EM algorithms for log-concave densities based on complete and censored data.
  * Efron et al. (2004) Efron, B., T. Hastie, I. Johnstone, and R. Tibshirani (2004).  Least angle regression.  Ann. Statist. 32(2), 407–499.  With discussion, and a rejoinder by the authors.
  * Fan and Li (2001) Fan, J. and R. Li (2001).  Variable selection via nonconcave penalized likelihood and its oracle properties.  J. Amer. Statist. Assoc. 96(456), 1348–1360.
  * Fan et al. (2011) Fan, J., A. Maity, Y. Wang, and Y. Wu (2011).  Analyzing mergers and acquisition data using parametrically guided generalized additive models.  submitted.
  * Frank and Friedman (1993) Frank, I. E. and J. H. Friedman (1993).  A statistical view of some chemometrics regression tools.  Technometrics 35(2), 109–135.
  * Friedman (2008) Friedman, J. (2008).  Fast sparse regression and classification.  http://www-stat.stanford.edu/ jhf/ftp/GPSpaper.pdf.
  * Friedman et al. (2000) Friedman, J., T. Hastie, and R. Tibshirani (2000).  Additive logistic regression: a statistical view of boosting.  The Annals of Statistics 28(2), 337–407.
  * Friedman et al. (2008) Friedman, J., T. Hastie, and R. Tibshirani (2008, July).  Sparse inverse covariance estimation with the graphical lasso.  Biostatistics 9(3), 432–441.
  * Fu (1998) Fu, W. J. (1998).  Penalized regressions: the bridge versus the lasso.  Journal of Computational and Graphical Statistics 7(3), 397–416.
  * Goodnight (1979) Goodnight, J. H. (1979).  A tutorial on the sweep operator.  Amer. Statist. 33(3), pp. 149–158.
  * Groeneboom and Wellner (1992) Groeneboom, P. and J. A. Wellner (1992).  Information bounds and nonparametric maximum likelihood estimation, Volume 19 of DMV Seminar.  Basel: Birkhäuser Verlag.
  * Hastie and Tibshirani (1993) Hastie, T. and R. Tibshirani (1993).  Varying-coefficient models.  J. Roy. Statist. Soc. Ser. B 55(4), 757–796.  With discussion and a reply by the authors.
  * Jennrich (1977) Jennrich, R. (1977).  Stepwise regression.  In Statistical Methods for Digital Computers, pp. 58–75. New York: Wiley-Interscience.
  * Jongbloed (1998) Jongbloed, G. (1998).  The iterative convex minorant algorithm for nonparametric estimation.  J. Comput. Graph. Statist. 7(3), 310–321.
  * Kim et al. (2009) Kim, S.-J., K. Koh, S. Boyd, and D. Gorinevsky (2009).  $l_{1}$ trend filtering.  SIAM Rev. 51(2), 339–360.
  * Lange (2010) Lange, K. (2010).  Numerical Analysis for Statisticians (Second ed.).  Statistics and Computing. New York: Springer.
  * Lawson and Hanson (1987) Lawson, C. L. and R. J. Hanson (1987).  Solving Least Squares Problems (New edition ed.).  Classics in Applied Mathematics. Society for Industrial Mathematics.
  * Little and Rubin (2002) Little, R. J. A. and D. B. Rubin (2002).  Statistical Analysis with Missing Data (Second ed.).  Wiley Series in Probability and Statistics. Hoboken, NJ: Wiley-Interscience [John Wiley & Sons].
  * Magnus and Neudecker (1999) Magnus, J. R. and H. Neudecker (1999).  Matrix Differential Calculus with Applications in Statistics and Econometrics.  Wiley Series in Probability and Statistics. Chichester: John Wiley & Sons Ltd.
  * Mardia et al. (1979) Mardia, K. V., J. T. Kent, and J. M. Bibby (1979).  Multivariate Analysis.  London: Academic Press [Harcourt Brace Jovanovich Publishers].  Probability and Mathematical Statistics: A Series of Monographs and Textbooks.
  * McCullagh and Nelder (1983) McCullagh, P. and J. A. Nelder (1983).  Generalized Linear Models.  Monographs on Statistics and Applied Probability. London: Chapman & Hall.
  * Nocedal and Wright (2006) Nocedal, J. and S. J. Wright (2006).  Numerical Optimization (Second ed.).  Springer Series in Operations Research and Financial Engineering. New York: Springer.
  * Osborne et al. (2000) Osborne, M. R., B. Presnell, and B. A. Turlach (2000).  A new approach to variable selection in least squares problems.  IMA J. Numer. Anal. 20(3), 389–403.
  * Pal et al. (2007) Pal, J. K., M. Woodroofe, and M. Meyer (2007).  Estimating a Polya frequency function2.  In Complex datasets and inverse problems, Volume 54 of IMS Lecture Notes Monogr. Ser., pp. 239–249. Beachwood, OH: Inst. Math. Statist.
  * Park and Hastie (2007) Park, M. Y. and T. Hastie (2007).  $L_{1}$-regularization path algorithm for generalized linear models.  J. R. Stat. Soc. Ser. B Stat. Methodol. 69(4), 659–677.
  * Robertson et al. (1988) Robertson, T., F. T. Wright, and R. L. Dykstra (1988).  Order Restricted Statistical Inference.  Wiley Series in Probability and Mathematical Statistics: Probability and Mathematical Statistics. Chichester: John Wiley & Sons Ltd.
  * Rosset and Zhu (2007) Rosset, S. and J. Zhu (2007).  Piecewise linear regularized solution paths.  Ann. Statist. 35(3), 1012–1030.
  * Rufibach (2010) Rufibach, K. (2010, June).  An active set algorithm to estimate parameters in generalized linear models with ordered predictors.  Computational Statistics & Data Analysis 54(6), 1442–1456.
  * Ruszczyński (2006) Ruszczyński, A. (2006).  Nonlinear Optimization.  Princeton, NJ: Princeton University Press.
  * Shivdasani and Wang (2009) Shivdasani, A. and Y. Wang (2009).  Did structured credit fuel and LBO boom?  http://ssrn.com/abstract=1394421.
  * Silvapulle and Sen (2005) Silvapulle, M. J. and P. K. Sen (2005).  Constrained Statistical Inference: Inequality, Order, and Shape Restrictions.  Wiley Series in Probability and Statistics. Hoboken, NJ: Wiley-Interscience [John Wiley & Sons].
  * Tibshirani (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  J. Roy. Statist. Soc. Ser. B 58(1), 267–288.
  * Tibshirani et al. (2005) Tibshirani, R., M. Saunders, S. Rosset, J. Zhu, and K. Knight (2005).  Sparsity and smoothness via the fused lasso.  J. R. Stat. Soc. Ser. B Stat. Methodol. 67(1), 91–108.
  * Tibshirani and Taylor (2011) Tibshirani, R. and J. Taylor (2011).  The solution path of the generalized lasso.  Ann. Statist. to appear.
  * Walther (2002) Walther, G. (2002).  Detecting the presence of mixing with multiscale maximum likelihood.  J. Amer. Statist. Assoc. 97(458), 508–513.
  * Walther (2009) Walther, G. (2009).  Inference and modeling with log-concave distributions.  Statist. Sci. 24(3), 319–327.
  * Wu (2011) Wu, Y. (2011).  An ordinary differential equation-based solution path algorithm.  Journal of Nonparametric Statistics 23, 185–199.
  * Yuan (2008) Yuan, M. (2008).  Efficient computation of $\ell_{1}$ regularized estimates in Gaussian graphical models.  J. Comput. Graph. Statist. 17(4), 809–826.
  * Zhou and Lange (2011) Zhou, H. and K. Lange (2011).  A path algorithm for constrained estimation.  arXiv:1103.3738.
  * Zou et al. (2007) Zou, H., T. Hastie, and R. Tibshirani (2007).  On the “degrees of freedom” of the lasso.  Ann. Statist. 35(5), 2173–2192.