A generalized lasso formulation of weighted and regularized (nearly-)isotonic regression on graphs

Introduction

In the area of constrained inference the basic and simplest problem is isotonic regression, i.e. the $\ell_2$-projection of a given vector of data points $\mathbf{y} \in\mathbb{R}^n$ onto the cone of non-decreasing vectors in $\mathbb{R}^n$,

A general introduction can be found in [14]. The nearly-isotonic regression, introduced in [20] and studied in detail in [11], is the less restrictive version

where $(t)_+ = \max\{t,0\}$, in which the constraint is replaced by a penalty on its violations. Its combination with the fused lasso and the extension to a general partial order were studied in [12], where in particular the dual of the problem was shown to be a quadratic program with box constraints.

The subject of penalized estimation starts with the $\ell_1$-penalisation, i.e. the lasso [18], and the $\ell_2$-penalisation, i.e. ridge regression [9]. The $\ell_1$-penalisation of the successive differences of the coefficients gives the fusion approximator [13], and its combination with the lasso is the fused lasso [19]. These estimators were studied and solved for a general graph in [4, 8, 21]. In [21] they appear as special cases of the generalized lasso, in which the $\ell_1$-penalty is imposed on $D\beta$ for a general penalty matrix $D$, and trend filtering on graphs [22] is another example. The lasso with linear constraints on the coefficients, called the constrained lasso, and its relation to the generalized lasso were studied in [5]. There the constraint matrices are assumed to have full row rank, and for isotonic regression this holds only if the graph of the partial order is a tree, which is not the case, for example, for the bimonotone order on a grid. Next, the statistical properties of isotonic regression in general dimensions were studied in [3, 7], and algorithms for it can be found, for example, in [2, 16]. Several regularized versions of isotonic regression have been proposed. Generalizations of nearly-isotonic regression in the one-dimensional case were studied in [23], where a dynamic programming algorithm is proposed, and in [6], where the squares of the positive parts of the differences and a lasso penalty are added to the penalty in (2).

In this paper we study the nearly-isotonic regression with respect to a general partial order in the following modified form: the sum of squares is weighted, with some of the weights possibly equal to zero, and a quadratic penalty is added. This modification is motivated by [1], where the authors study least squares estimation under bimonotonicity constraints for the data observed on an incomplete grid. There the weights are the numbers of observations at the grid points, and they are equal to zero at the empty cells. In this case the problem is not strictly convex, and the authors add a small quadratic penalty, which they call light regularization. The quadratic penalty is also the one of the smoothed monotonic regression of [17] in the one-dimensional case. The main result of the paper is that the modified problem can be written as a generalized lasso problem in the sense of [21], with the same penalty matrix and with a design matrix of full column rank, and can therefore be solved by the algorithms developed for the generalized lasso. We also derive the dual problem.

Notation

We start with basic definitions of partial order and isotonic regression. Let $I = \{i_1,\ldots,i_n\}$ be some index set. A binary relation $\preceq$ on $I$ is called partial order if

• it is reflexive, i.e. $j \preceq j$ for all $j \in I$,

• it is transitive, i.e. $j_1,j_2,j_3 \in I$, $j_1 \preceq j_2$ and $j_2 \preceq j_3$ imply $j_1 \preceq j_3$,

• it is antisymmetric, i.e. $j_1,j_2 \in I$, $j_1 \preceq j_2$ and $j_2 \preceq j_1$ imply $j_1 = j_2$.

Further, a vector $\beta\in\mathbb{R}^n$ indexed by $I$ is called isotonic with respect to the partial order $\preceq$ on $I$ if $j_1 \preceq j_2$ implies $\beta_{j_1} \le\beta_{j_2}$. We denote the set of all isotonic vectors in $\mathbb{R}^n$ with respect to $\preceq$ by $\mathcal{B}^I$. This set is a closed convex cone in $\mathbb{R}^n$, and it is also called the isotonic cone. Next, the isotonic regression of a vector $y \in\mathbb{R}^n$ indexed by the partially ordered index set $I$ is given by

which for the natural order on $I = \{1,\ldots,n\}$ is (1).

For any partial order relation $\preceq$ on $I$ there exists a directed graph $G=(V,E)$, with $V=I$ and $E$ the minimal set of edges $(j_1,j_2)$ with $j_1 \preceq j_2$, such that $\beta$ is isotonic with respect to $\preceq$ if and only if $\beta_{l_1} \le\beta_{l_2}$ whenever $E$ contains a chain of edges from $l_1$ to $l_2$. We label the vertices by $\{1,\ldots,n\}$ and the edges by $\{1,\ldots,m\}$, and let $D \in\mathbb{R}^{m \times n}$ be the oriented incidence matrix of $G$,

so that $D\beta$ is the vector of the differences of $\beta$ along the edges of $G$. Therefore, the isotonic cone is

where the inequality is meant coordinatewise. Further, for $z \in\mathbb{R}^{m}$ we write

For the natural order on $I=\{1,\ldots,n\}$ the edges join consecutive indices, and $\lVert D\beta\rVert_{+}$ is the penalty in (2). For the bimonotone order on a grid, introduced in [1], the edges join horizontally and vertically adjacent grid points.

General statement of the problem

Let $y \in\mathbb{R}^{n}$ be a signal indexed by $I$, let $W$ be a diagonal matrix with nonnegative weights $w_1,\ldots,w_n$ on the diagonal, not all equal to zero, let $S$ be a given matrix with $n$ columns and $\ell$ rows, and let $\lambda_{NI} \ge0$. We consider the estimator

under the assumption that

is strictly convex.

Theorem 2 and Propositions 3 and 6 below hold with $D$ replaced by $D_{\omega}$, since their proofs do not use the structure of $D$.

The structure of (5) is that of the “incomplete layout with light regularization” in Setting 2b of [1]. There, the observations of a bimonotone regression function are reduced to a rectangular grid: $w_u$ is the number of observations at the grid point $u$ and $Z_u$ their average, and the least squares criterion becomes

When the layout is incomplete, i.e. some grid points have no observations and $w_u=0$, this function is not strictly convex, and the least squares fit is not uniquely determined at the empty cells. For this reason the authors consider the minimisation of

over the bimonotone cone, where $P$ is a convex quadratic function chosen so that the sum is strictly convex, and $\lambda_S>0$ is small, its role being only to make the Hessian well conditioned.

The authors discuss two choices of $P$, the Tychonov regularisation with a reference value $\theta_0$ and the sum of squared differences along the edges of the grid,

They prefer the latter, because it interpolates more smoothly into the empty cells. In the notation of (5) the latter penalty is $\lVert S\beta\rVert_2^2$ with

and their requirement on $P$ is the assumption (6). In Lemma 4 below we verify (6) for this choice. Problem (5) is thus the nearly-isotonic relaxation of the problem of [1], in which the isotonic constraint is replaced by the penalty, and the bimonotone order is replaced by an arbitrary partial order. On the chain graph, with unit weights and with the constraint kept, the same quadratic part gives the smoothed monotonic regression of [17] with equal penalty parameters.

Organisation of the paper

The rest of the paper is organised as follows. In Section 2 we relate the estimator to the constrained estimator in which the penalty $\lambda_{NI}\lVert D\beta\rVert_+$ is replaced by the constraint $D\beta\le0$. In Section 3 we show that (5) is a generalized lasso problem. Section 4 contains the dual problem. In Section 5 we verify assumption (6) for the penalty (7). Section 6 is devoted to the computational aspects. The paper closes with a discussion in Section 7.

Relation to regularized isotonic regression

In this section we compare the estimator (5) with the estimator obtained when the penalty $\lambda_{NI}\lVert D\beta\rVert_+$ in (5) is removed and the isotonic constraint $D\beta\le0$ is imposed instead,

For the grid and $S=\sqrt{\lambda_S}D$ this is the estimator of [1], and for the chain it is the smoothed monotonic regression of [17]. For nearly-isotonic regression, i.e. for $W=\frac{1}{2}I_n$ and $S=0$, the property that the estimator coincides with isotonic regression for all sufficiently large $\lambda_{NI}$ is mentioned in [11] and [20]. The same holds for the estimator (5). Under (6) the objective functions in (5) and (8) are strictly convex and coercive, so both problems have unique solutions. Further, the objective function in (5) is the exact penalty function of the problem (8) with the penalty parameter $\lambda_{NI}$, and the constraints in (8) are linear. Therefore, by the classical exact penalty result, see [10] and the references therein, there exists $\lambda_0\ge0$ such that $\widehat{\beta}=\widehat{\beta}^{CIS}$ for all $\lambda_{NI}\ge\lambda_0$, in particular, for nearly-isotonic regression with respect to a general partial order.

Thus, the estimator (5) is the regularized isotonic regression (8) for all sufficiently large values of $\lambda_{NI}$, and the results of the following sections apply to the estimator (8) as well.

Solution via the generalized lasso

First, we rewrite the objective function in (5).

Proof. The function in (6) is a quadratic function with constant Hessian $2Q$, and a quadratic function is strictly convex if and only if its Hessian is positive definite. This gives the last statement of the lemma.

For (9), since $W$ is symmetric, $(\mathbf{y} - \boldsymbol{\beta})^{\mathsf{T}}W(\mathbf{y} - \boldsymbol{\beta}) = \boldsymbol{\beta}^{\mathsf{T}}W\boldsymbol{\beta} - 2(W\mathbf{y})^{\mathsf{T}}\boldsymbol{\beta} + \mathbf{y}^{\mathsf{T}}W\mathbf{y}$, and $\lVert S\boldsymbol{\beta}\rVert_{2}^{2} = \boldsymbol{\beta}^{\mathsf{T}}S^{\mathsf{T}}S\boldsymbol{\beta}$. Together these give $\boldsymbol{\beta}^{\mathsf{T}}Q\boldsymbol{\beta} - 2(W\mathbf{y})^{\mathsf{T}}\boldsymbol{\beta} + \mathbf{y}^{\mathsf{T}}W\mathbf{y}$. For $t \in\mathbb{R}$ one has $(t)_{+} = \frac{1}{2}|t| + \frac{1}{2}t$. Therefore,

which finishes the proof. $\square$

Next, we provide the solution to the problem (5).

The pair $X$, $\widetilde{y}$ in (10) is not the only possible choice. As the proof shows, the statement of Theorem 2 holds for any $X \in\mathbb{R}^{p \times n}$ and $\widetilde{y} \in\mathbb{R}^{p}$ such that

For example, (13) holds for the square matrix $X$ and the vector $\widetilde{y}$ given by

Finally, for nearly-isotonic regression, i.e. for $W = \frac{1}{2}I_n$ and $S = 0$, (13) holds for

and (12) becomes the relation

between nearly-isotonic regression and the fusion approximator, proved in [12].

The dual problem

For a fixed $\lambda_{NI}$ we can also solve the problem (5) numerically through its dual, which is a quadratic program with box constraints.

Proof. First, following [21], we write the problem (5) in the following way

Further, the Lagrangian is given by

where $u \in\mathbb{R}^{m}$ is a dual variable. Note that

and, since $Q$ is positive definite,

Next, the dual function is given by

and, therefore, the dual problem is

which is equivalent to (15). Lastly, since the objective function is finite and convex on $\mathbb{R}^{n}\times\mathbb{R}^{m}$ and the constraint $D\beta=z$ is affine and feasible, strong duality holds, and taking first derivative of the Lagrangian $L(\beta,z,\widehat{u})$ with respect to $\beta$ we get the following relation between $\widehat{\beta}$ and $\widehat{u}$

If $D$ does not have full row rank, the minimiser in (15) need not be unique. However, the objective function in (15) is a strictly convex function of $D^{\mathsf{T}}u$, therefore, $D^{\mathsf{T}}\widehat{u}$ is unique, and (13) does not depend on the choice of $\widehat{u}$. Note also that $\widetilde{y}_d$ is the solution to (5) for $\lambda_{NI}=0$, so (13) represents the estimator as the quadratic smoother minus a correction determined by the dual problem. For nearly-isotonic regression we have $\widetilde{y}_d=y$ and $\widetilde{D}=D$, and Proposition 3 becomes Theorem 1 of [12] without the fusion penalty.

The strict convexity assumption

In this section we verify the assumption (6) for the penalty (6), i.e. for $S=\sqrt{\lambda_S}D$.

Therefore, for a connected graph the standing assumption that not all weights are equal to zero is necessary and sufficient. In particular, for the bimonotone grid, which is connected, the assumption (6) holds for every incomplete layout with at least one observation and for every $\lambda_S>0$. Note that without the quadratic penalty, i.e. for $S=0$, the Hessian of the quadratic part of (5) is $2W$, which is singular as soon as $w_i=0$ for some $i$.

If the graph is not connected, then for $S=\sqrt{\lambda_S}D$ all terms in (5) are sums over the connected components, because no edge joins two different components. Therefore, the problem (5) splits into independent problems, one for each component, and Lemma 4 can be applied to each of them.

Proof. Let $C_1,\ldots,C_K$ be the connected components of the undirected graph underlying $G$. Since no edge joins two different components, after a permutation of the vertices the matrix $Q=W+\lambda_SD^{\mathsf T}D$ is block diagonal with the blocks $Q_k=W_k+\lambda_SD_k^{\mathsf T}D_k$, $k=1,\ldots,K$, where $W_k$ is the diagonal matrix of the weights $w_i$, $i\in C_k$, and $D_k$ is the incidence matrix of the subgraph of $G$ with the vertex set $C_k$. Therefore, $Q$ is positive definite if and only if every $Q_k$ is positive definite, and the statement follows from Lemma 4 applied to each component. $\square$

Finally, the Tychonov penalty of [1] shrinks the estimator towards the reference value $\theta_0$, and for $\theta_0\ne0$ it is not of the form $\lVert S\beta\rVert_2^2$. Nevertheless, the corresponding problem

reduces to (5). Indeed, let $\boldsymbol{\beta}'=\boldsymbol{\beta}-\theta_0\boldsymbol{1}_n$ and $\boldsymbol{y}'=\boldsymbol{y}-\theta_0\boldsymbol{1}_n$. Then

where the last equality holds because $D\boldsymbol{1}_n=0$. Therefore, the solution is

where $\widehat{\boldsymbol{\beta}}'$ is the solution to (5) with the data vector $\boldsymbol{y}'$ and $S=\sqrt{\lambda_S}I_n$. For this choice of $S$ the matrix $Q$ is positive definite for any weights, so the assumption (6) always holds.

Computational aspects

In this section we discuss the computation of the estimator (5). We do not provide a simulation study, because the algorithms below are well established, and their performance has been studied extensively, see [12] for the nearly-isotonic case, and also [15], [21].

First, by Theorem 2 the estimator is the solution to the generalized lasso problem (11). Therefore, we can use the algorithm of [21] for the problem (5), for an arbitrary partial order, and the condition of full column rank of the design matrix, required there, is satisfied automatically. Note, that the data vector $\widetilde{\boldsymbol{y}}$ depends on $\lambda_{NI}$, therefore, the path of [21] in the lasso parameter with fixed $\widetilde{\boldsymbol{y}}$ is not the path of (5) in $\lambda_{NI}$. Further, among all pairs $X,\widetilde{\boldsymbol{y}}$ which satisfy (12), the choice (10) is preferable in computations, because the design matrix in (10) is sparse whenever $S$ is sparse, while $Q^{1/2}$ is dense in general.

Second, the estimator can be computed through the dual problem (15). However, the matrix $Q^{-1}$ in (14) is dense in general, even if $Q$ is sparse. To avoid computing $Q^{-1}$, we keep $\boldsymbol{\beta}$ as a variable in the dual problem.

Proof. The constraint $Q\beta+\frac{1}{2}D^\mathsf{T}u=Wy$ gives $\beta=\widetilde{y}_d-\widetilde{D}^\mathsf{T}u$, and then $\beta^\mathsf{T}Q\beta=\lVert\widetilde{y}_d-\widetilde{D}^\mathsf{T}u\rVert_Q^2$. Therefore, $u^*$ is a solution to (15) and $\beta^*=\widetilde{y}_d-\widetilde{D}^\mathsf{T}u^*$, which is equal to $\hat{\beta}$ by Proposition 3.

The problem (18) contains only the matrices $Q$ and $D$. Every row of the incidence matrix $D$ has two nonzero elements, and $Q$ is sparse whenever $S^\mathsf{T}S$ is sparse. In particular, for the penalty (7) the matrix $Q=W+\lambda_SD^\mathsf{T}D$ has at most $n+2m$ nonzero elements, so it is sparse when the number of edges $m$ is of the order of $n$, as for the chain and for the grid. In this case (18) can be solved by quadratic programming solvers for sparse problems, for example by OSQP [15], as it was done for the nearly-isotonic case in [12].

Lastly, if the constrained estimator (8) is of interest, then, by § 2, it can be obtained from (5) with a sufficiently large value of $\lambda_{NI}$.

Conclusion and discussion

In this paper we studied the nearly-isotonic regression with respect to a general partial order, with a weighted sum of squares, where some of the weights can be equal to zero, and with an additional quadratic penalty. In our opinion, the main result of the paper is Theorem 2, where we proved that this problem can be written as a generalized lasso problem with a design matrix of full column rank and the same penalty matrix. Therefore, the estimator can be computed by the algorithms developed for the generalized lasso. Further, we derived the dual problem and wrote it as a sparse quadratic program, which can be solved by standard solvers for sparse problems. Also, we showed that for a sufficiently large penalty parameter the estimator coincides with the regularized isotonic regression.

One of the possible directions is to study the degrees of freedom of the estimator. In Theorem 2 of [21], for a design matrix of full column rank and Gaussian data with a covariance matrix proportional to the identity matrix, an unbiased estimate of the degrees of freedom of the fitted values of the generalized lasso is derived. In our Theorem 2 the design matrix has full column rank, but the fitted values $X\hat{\beta}$ are not the estimate of interest, and the data vector $\widetilde{y}$ depends on $\lambda_{NI}$. Moreover, even for Gaussian observations with a covariance matrix proportional to the identity matrix, the covariance matrix of $\widetilde{y}$ is not of this form. Therefore, the result of [21] does not apply directly, and it is interesting to derive the degrees of freedom of the estimator $\hat{\beta}$ itself, in the spirit of Theorem 6 in [12].

Another direction is to develop a path algorithm with respect to $\lambda_{NI}$. Since the data vector $\widetilde{y}$ in Theorem 2 depends on $\lambda_{NI}$, the path algorithm of [21] does not give it directly. In particular, it is interesting to find the smallest value of $\lambda_0$ in § 2, for which the estimator coincides with the regularized isotonic regression.