Perspective Functions: Proximal Calculus and Applications in High-Dimensional Statistics

Introduction

Perspective functions appear, often implicitly, in various problems in areas as diverse as statistics, control, computer vision, mechanics, game theory, information theory, signal recovery, transportation theory, machine learning, disjunctive optimization, and physics (see the companion paper for a detailed account). In the setting of a real Hilbert space ${\mathcal{G}}$, the most useful form of a perspective function, first investigated in Euclidean spaces in , is the following.

Let $\varphi\colon{\mathcal{G}}\to\left]-\infty,+\infty\right]$ be a proper lower semicontinuous convex function and let $\text{\rm rec}\,\varphi$ be its recession function. The perspective of $\varphi$ is

which hinges on the perspective function of the squared Euclidean norm (see for further discussion).

In the literature, problems involving perspective functions are typically solved with a wide range of ad-hoc methods. Despite the ubiquity of perspective functions, no systematic structuring framework has been available to approach these problems. The goal of this paper is to fill this gap by showing that they are amenable to solution by proximal methods, which offer a broad array of splitting algorithms to solve complex nonsmooth problems with attractive convergence guarantees . The central element in the successful implementation of a proximal algorithm is the ability to compute the proximity operator of the functions present in the optimization problem. We therefore propose a systematic investigation of proximity operators for perspective functions and show that the proximal framework can efficiently solve perspective-function based problems, unveiling in particular new applications in high-dimensional statistics.

In Section 2, we introduce basic concepts from convex analysis and review essential properties of perspective function. We then study the proximity operator of perspective functions in Sections 3. We establish a characterization of the proximity operator and then provide examples of computation for concrete instances. Section 4 unveils new applications of perspective functions in high-dimensional statistics and demonstrates the flexibility and potency of the proposed framework to both model and solve complex problems in statistical data analysis.

Notation and background

Throughout, ${\mathcal{H}}$, ${\mathcal{G}}$, and ${\mathcal{K}}$ are real Hilbert spaces and ${\mathcal{H}}\oplus{\mathcal{G}}$ denotes their Hilbert direct sum. The symbol $\|\cdot\|$ denotes the norm of a Hilbert space and ${\langle{{\cdot}\mid{\cdot}}\rangle}$ the associated scalar product. The closed ball with center $x\in{\mathcal{K}}$ and radius $\rho\in\left]0,+\infty\right[$ is denoted by $B(x;\rho)$.

A function $f\colon{\mathcal{K}}\to\left]-\infty,+\infty\right]$ is proper if \text{\rm dom}\,f=\big{\{}{x\in{\mathcal{K}}}~{}\big{|}~{}{f(x)<{+\infty}}\big{\}}\neq{\varnothing}, coercive if $\lim_{\|x\|\to{+\infty}}f(x)={+\infty}$, and supercoercive if $\lim_{\|x\|\to{+\infty}}f(x)/\|x\|={+\infty}$. Denote by $\Gamma_{0}({\mathcal{K}})$ the class of proper lower semicontinuous convex functions from ${\mathcal{K}}$ to $\left]-\infty,+\infty\right]$, and let $f\in\Gamma_{0}({\mathcal{K}})$. The conjugate of $f$ is the function

It also belongs to $\Gamma_{0}({\mathcal{K}})$ and $f^{**}=f$. The subdifferential of $f$ is the set-valued operator

If $f$ is Gâteaux differentiable at $x\in\text{\rm dom}\,f$ with gradient $\nabla f(x)$, then

Let $z\in\text{\rm dom}\,f$. The recession function of $f$ is

The infimal convolution operation is denoted by $\,\square\,$. Now let $C$ be a subset of ${\mathcal{K}}$. Then

is the support function of $C$. If $C$ is nonempty, closed, and convex then, for every $x\in{\mathcal{K}}$, there exists a unique point $P_{C}x\in C$, called the projection of $x$ onto $C$, such that $\|x-P_{C}x\|=d_{C}(x)$. We have

For further background on convex analysis, see .

2 Proximity operators

The proximity operator of $f\in\Gamma_{0}({\mathcal{K}})$ is

This operator was introduced by Moreau in 1962 to model problems in unilateral mechanics. In , it was shown to play an important role in the investigation of various data processing problems, and it has become increasingly prominent in the general area of data analysis . We review basic properties and refer the reader to for a more complete account.

Let $f\in\Gamma_{0}({\mathcal{K}})$. Then

If $C$ is a nonempty closed convex subset of ${\mathcal{K}}$, then

Let $\gamma\in\left]0,+\infty\right[$. The Moreau decomposition of $x\in{\mathcal{K}}$ is

Let $(\Omega,{\mathcal{F}},\mu)$ be a complete $\sigma$-finite measure space, let $\mathsf{K}$ be a separable real Hilbert space, and let $\psi\in\Gamma_{0}({\mathsf{K}})$. Suppose that ${\mathcal{K}}=L^{2}((\Omega,{\mathcal{F}},\mu);{\mathsf{K}})$ and that $\mu(\Omega)<{+\infty}$ or $\psi\geqslant\psi(0)=0$. Set

Let $x\in{\mathcal{K}}$ and define, for $\mu$-almost every $\omega\in\Omega$, $p(\omega)=\text{\rm prox}_{\psi}x(\omega)$. Then $p=\text{\rm prox}_{\Phi}x$.

Proof. By [1, Proposition 9.32], $\Phi\in\Gamma_{0}({\mathcal{K}})$. Now take $x$ and $p$ in ${\mathcal{K}}$. Then it follows from (2.14) and [1, Proposition 16.50] that $p(\omega)=\text{\rm prox}_{\Phi}x(\omega)$ $\mu$-a.e. $\Leftrightarrow$ $x(\omega)-p(\omega)\in\partial\psi(p(\omega))$ $\mu$-a.e. $\Leftrightarrow$ $x-p\in\partial\Phi(p)$. $\Leftrightarrow$ $p=\text{\rm prox}_{\Phi}x$.

Let $D\neq\{0\}$ be a nonempty closed convex subset of ${\mathcal{K}}$, let $x\in{\mathcal{K}}$, and let $\gamma\in\left]0,+\infty\right[$. Set $f=\|\cdot\|+\sigma_{D}$ and $C=\gamma D$. Then

If, in addition, $D$ is a cone and $K$ denotes its polar cone, then $f=\|\cdot\|+\iota_{K}$ and

Proof. Using elementary convex analysis, we obtain

Hence, it follows from (2.16) and (2.15) that

However by [1, Propositions 28.1(ii) and 28.10],

Upon combining (2.21) and (2.22), we arrive at (2.18). Now suppose that, in addition, $D$ is a cone. Then $C=D$, $\sigma_{D}=\iota_{K}$, and (2.16) yields $\operatorname{Id}\,-P_{D}=P_{K}$. Altogether, (2.18) reduces to (2.19).

3 Perspective functions

We review here some essential properties of perspective functions.

Let $\varphi\in\Gamma_{0}({\mathcal{G}})$. Then the following hold:

We refer to the companion paper for further properties of perspective functions as well as examples. Here are two important instances of (composite) perspective functions that will play a central role in Section 4.

Proof. This is a special case of [7, Example 4.2].

Proximity operator of a perspective function

We start with a characterization of the proximity operator of a perspective function when $\text{\rm dom}\,\varphi^{*}$ is open.

Suppose that $\eta+\gamma\varphi^{*}(y/\gamma)\leqslant 0$. Then $\text{\rm prox}_{\gamma\widetilde{\varphi}}(\eta,y)=(0,0)$.

Suppose that $\text{\rm dom}\,\varphi^{*}$ is open and that $\eta+\gamma\varphi^{*}(y/\gamma)>0$. Then

where $p$ is the unique solution to the inclusion

If $\varphi^{*}$ is differentiable at $p$, then $p$ is characterized by $y=\gamma p+(\eta+\gamma\varphi^{*}(p))\nabla\varphi^{*}(p)$.

Proof. It follows from Lemma 2.3(ii) that

Since $\varphi\in\Gamma_{0}({\mathcal{G}})$, we have $\varphi^{*}\in\Gamma_{0}({\mathcal{G}})$. Therefore, $C$ is a nonempty closed convex set. In turn, we derive from [9, Proposition 3.2] that $\text{\rm prox}_{\gamma\widetilde{\varphi}}=\text{\rm prox}_{\sigma_{\gamma C}}$ is a proximal thresholder on $\gamma C$ in the sense that

(ii): Set $(\chi,q)=\text{\rm prox}_{\gamma\widetilde{\varphi}}(\eta,y)$ and $p=(y-q)/\gamma$. It follows from (2.14) that $(\chi,q)\in\text{\rm dom}\,(\gamma\partial\widetilde{\varphi})$ and from (3.4) that $(\chi,q)\neq(0,0)$. Hence, we deduce from Lemma 2.3(iv) that $\chi>0$. Furthermore, we derive from (2.14) and Lemma 2.3(iii) that $(\chi,q)$ is characterized by

Hence, we derive from (3.6) that $\varphi^{*}(p)=(\chi-\eta)/\gamma$, i.e.,

Altogether, we have established the characterization (3.1)–(3.2), while the assertion concerning the differentiable case follows from (2.6).

Here is an alternative proof of Theorem 3.1. It follows from Lemma 2.3(ii) that

is a nonempty closed convex set. Hence, using (2.16) and (2.15), we obtain

Now set $(\pi,p)=P_{C}(\eta/\gamma,y/\gamma)$. We deduce from (2.15), (2.16), and (2.12) that $(\pi,p)$ is characterized by

(i): We have $(\eta/\gamma,y/\gamma)\in C$. Hence, $(\pi,p)=(\eta/\gamma,y/\gamma)$ and (3.11) yields $\text{\rm prox}_{\gamma\widetilde{\varphi}}(\eta,y)=(0,0)$.

Now let $z\in\text{\rm dom}\,\varphi^{*}$ and let $\zeta\in\left]{-\infty},-\varphi^{*}(z)\right[$. Then $h(\zeta,z)<0$. Therefore, we derive from [1, Lemma 26.17 and Proposition 16.8] and (3.13) that

Hence, if $\pi<-\varphi^{*}(p)$, then (3.12) yields $(\eta/\gamma-\pi,y/\gamma-p)=(0,0)$ and therefore $(\eta/\gamma,y/\gamma)=(\pi,p)\in C$, which is impossible since $(\eta/\gamma,y/\gamma)\notin C$. Thus, the characterization (3.12) becomes

that is, $y\in\gamma p+(\eta+\gamma\varphi^{*}(p))\partial\varphi^{*}(p)$.

The next example is based on distance functions.

We have $\text{\rm dom}\,\varphi^{*}={\mathcal{G}}$ and, in view of Theorem 3.1(ii), we need only assume that $\eta+\gamma\varphi^{*}(y/\gamma)>0$, i.e.,

In view of Remark 3.2, the normal cone to the set $C$ of (3.10) at $(0,0)$ is

So, for every $(\eta,y)\in K$, $P_{C}(\eta/\gamma,y/\gamma)=(0,0)$ and $\text{\rm prox}_{\gamma\widetilde{\varphi}}(\eta,y)=(\eta,y)$. Now suppose that $(\eta,y)\notin K$. Then $p\neq 0$ and, taking the norm in the upper line of (3.20), we obtain

Since $\phi^{*}$ is convex, $\theta$ is strongly convex and it therefore admits a unique minimizer $t$. Therefore $\psi(t)=\theta^{\prime}(t)=0$ and $\|p\|=t=\psi^{-1}(\|y\|/\gamma)$ is the unique solution to (3.22). In turn, (3.20) yields

and we obtain $\text{\rm prox}_{\gamma\widetilde{\varphi}}(\eta,y)$ via (3.1).

Next, we compute the proximity operator of a special case of the perspective function introduced in Lemma 2.5.

Then $\psi$ is invertible. Moreover, if $\eta+\gamma\phi^{*}(\|y/\gamma-v\|)>\gamma\delta$, set

In view of Theorem 3.1, it remains to assume that $\eta+\gamma\varphi^{*}(y/\gamma)>0$, i.e., $\eta+\phi^{*}(\|y/\gamma-v\|)>\gamma\delta$, and to show that the point $(t,p)$ provided by (3.28) satisfies

Since $\phi^{*}(0)=0$, we recover (3.29).

$y\neq\gamma v$: As seen in (3.33), $p\neq v$. Using (3.30) and (3.31), (3.32) can be rewritten as

Upon taking the norm on both sides of the second equality, we obtain

We note that, since $\phi^{*}$ is convex, $\psi$ is the derivative of the strongly convex function

Consequently, $\psi$ is strictly increasing [1, Proposition 17.13], hence invertible. It follows that $t=\psi^{-1}(\|y/\gamma-v\|)$. In turn, (3.36) yields (3.28).

If $\eta+\sqrt{\gamma^{2}+\|y\|^{2}}>0$, set

Proof. This is a special case of Corollary 3.5 with with $\delta=0$, $v=0$, and

It follows from [1, Example 13.2(vi) and Corollary 13.33] that $\phi^{*}\colon s\mapsto\sqrt{1+s^{2}}$. Hence, $\phi^{*^{\prime}}\colon s\mapsto s/\sqrt{1+s^{2}}$ and we derive (3.42) from (3.29).

Proof. This is a special case of Corollary 3.5 with $\phi=|\cdot|^{q}/\alpha$. Indeed, we derive from [1, Example 13.2(i) and Proposition 13.20(i)] that $\phi^{*}=\varrho|\cdot|^{q^{*}}/q^{*}$, which implies that (3.46)–(3.47) follow from (3.29).

Note that (3.49) can be solved explicitly via Cardano’s formula [4, Chapter 4] to obtain $t$.

We conclude this subsection by investigating integral functions constructed from integrands that are perspective functions.

Now let $x\in{\mathcal{H}}$ and $y\in{\mathcal{G}}$, and set, for $\mu$-almost every $\omega\in\Omega$, $(p(\omega),q(\omega))=\text{\rm prox}_{\widetilde{\varphi}}(x(\omega),y(\omega))$. Then $\text{\rm prox}_{\Phi}(x,y)=(p,q)$.

2 Further results

A convenient assumption in Theorem 3.1(ii) is that $\text{\rm dom}\,\varphi^{*}$ is open, as it allowed us to rule out the case when

and to reduce (3.14) to (3.15) using (3.13). In general, (3.13) has the form $N_{\text{\rm dom}\,h}(\pi,p)=\{0\}\times N_{\text{\rm dom}\,\varphi^{*}}p$ and, if $\text{\rm dom}\,\varphi^{*}$ is simple enough, explicit expressions can still be obtained. To shed more light on the case (3.53), consider the scenario in which $q\neq 0$ and $\text{\rm dom}\,\varphi^{*}$ is closed, and set $p=(y-q)/\gamma$. Then, in view of (2.14), (3.53) yields $(\eta/\gamma,p)\in\partial\widetilde{\varphi}(0,q)$. In turn, we derive from (2.23) that

and we infer from (2.11) that $p=P_{\text{\rm dom}\,\varphi^{*}}(y/\eta)$. Therefore,

and we note that the condition $q\neq 0$ means that $y\notin\gamma\,{\text{\rm dom}\,\varphi^{*}}$. We provide below examples in which $\text{\rm dom}\,\varphi^{*}$ is a simple proper closed subset of ${\mathcal{G}}$ and the proximity operator of the perspective function of $\varphi$ can be computed explicitly.

Suppose that $D\neq\{0\}$ is a nonempty closed convex cone in ${\mathcal{G}}$ and define

Since $\text{\rm dom}\,\vartheta={\mathcal{G}}$, we have \varphi^{*}=(\vartheta+\iota_{D})^{*}=\vartheta^{*}\mbox{\small\,\square\,}\iota_{D^{\ominus}}, where $D^{\ominus}$ is the polar cone of $D$ and (combine [1, Examples 13.2(vi) and 13.7])

Thus, \text{\rm dom}\,\varphi^{*}=\text{\rm dom}\,(\vartheta^{*}\mbox{\small\,\square\,}\iota_{D^{\ominus}})=\text{\rm dom}\,\vartheta^{*}+\text{\rm dom}\,\iota_{D^{\ominus}}=B(0;1)+D^{\ominus} is closed as the sum of two closed convex sets, one of which is bounded. As a result, since $D^{\ominus}\neq{\mathcal{G}}$,

where $\eta_{+}=\max\{0,\eta\}$ and $y_{+}$ is defined likewise componentwise.

The second example provides the proximity operator of the perspective function of the Huber function.

Following [7, Example 3.2], let $\rho\in\left]0,+\infty\right[$ and consider the perspective function

If $\eta+|y|^{2}/(2\gamma)\leqslant 0$ and $|y|\leqslant\gamma\rho$, then Theorem 3.1(i) yields $(\chi,q)=(0,0)$.

We have $\chi=0$ $\Leftrightarrow$ $\eta/\gamma\leqslant-\rho^{2}/2$. Hence, if $\eta\leqslant-\gamma\rho^{2}/2$ and $|y|>\gamma\rho$, (3.56) yields $(\chi,q)=(0,y-P_{[-\gamma\rho,\gamma\rho]}y)=(0,y-\gamma\rho\,\text{\rm sign}(y))$.

If $\eta>-\gamma\rho^{2}/2$ and $|y|>\rho\eta+\gamma\rho(1+\rho^{2}/2)$, then $(\eta/\gamma,y/\gamma)\in(-\rho^{2}/2,\rho\,\text{\rm sign}(y))+N_{C}(-\rho^{2}/2,\rho\,\text{\rm sign}(y))$ and therefore $P_{C}(\eta/\gamma,y/\gamma)=(-\rho^{2}/2,\rho\,\text{\rm sign}(y))$. Hence, (3.11) yields $(\chi,q)=(\eta+\gamma\rho^{2}/2,y-\gamma\rho\,\text{\rm sign}(y))$.

If $\eta>-\gamma\rho^{2}/2$ and $|y|\leqslant\rho\eta+\gamma\rho(1+\rho^{2}/2)$, then $(\chi,q)=\text{\rm prox}_{\gamma[|\cdot|^{2}/2]^{\sim}}(\eta,y)$ is obtained by setting $v=0$, $\delta=0$, and $\alpha=2$ in Example 3.8.

The last example concerns the Vapnik loss function.

Following [7, Example 3.4], let $\varepsilon\in\left]0,+\infty\right[$ and consider the perspective function

of the Vapnik $\varepsilon$-insensitive loss function

We have \varphi=d_{[-\varepsilon,\varepsilon]}=\iota_{[-\varepsilon,\varepsilon]}\mbox{\small\,\square\,}|\cdot| and therefore $\varphi^{*}=\varepsilon|\cdot|+\iota_{}$. Furthermore, (3.10) becomes

If $\eta+\varepsilon|y|\leqslant 0$ and $|y|\leqslant\gamma$, then Theorem 3.1(i) yields $(\chi,q)=(0,0)$.

We have $\chi=0$ $\Leftrightarrow$ $\eta/\gamma\leqslant-\varepsilon$. Hence, if $\eta\leqslant-\gamma\varepsilon$ and $|y|>\gamma$, (3.56) yields $(\chi,q)=(0,y-P_{[-\gamma,\gamma]}y)=(0,y-\gamma\,\text{\rm sign}(y))$.

If $\eta>-\gamma\varepsilon$ and $|y|>\varepsilon\eta+\gamma(1+\varepsilon^{2})$, then $(\eta/\gamma,y/\gamma)\in(-\varepsilon,\text{\rm sign}(y))+N_{C}(-\varepsilon,\text{\rm sign}(y))$ and therefore $P_{C}(\eta/\gamma,y/\gamma)=(-\varepsilon,\text{\rm sign}(y))$. Hence, (3.11) yields $(\chi,q)=(\eta+\gamma\varepsilon,y-\gamma\,\text{\rm sign}(y))$.

If $|y|>-\eta/\varepsilon$ and $\varepsilon\eta\leqslant|y|\leqslant\varepsilon\eta+\gamma(1+\varepsilon^{2})$, then $P_{C}(\eta/\gamma,y/\gamma)$ coincides with the projection of $(\eta/\gamma,y/\gamma)$ onto the half-space with outer normal vector $(1,\varepsilon\,\text{\rm sign}(y))$ and which has the origin on its boundary. As a result, (3.11) yields $(\chi,q)=((\eta+\varepsilon|y|)/(1+\varepsilon^{2}),\varepsilon(\eta+\varepsilon|y|)\text{\rm sign}(y)/(1+\varepsilon^{2}))$.

If $\eta\geqslant 0$ and $|y|\leqslant\varepsilon\eta$, then $P_{C}(\eta/\gamma,y/\gamma)=(0,0)$ and (3.11) yields $(\chi,q)=(\eta,y)$.

Applications in high-dimensional statistics

Sections 2 and 3 provide a unifying framework to model a variety of problems around the notion of a perspective function. By applying the results of Section 3 in existing proximal algorithms, we obtain efficient methods to solve complex problems. To illustrate this point, we focus on a specific application area: high-dimensional regression in the statistical linear model.

We consider the standard statistical linear model

where the fixed weights $w_{j}\in\left]0,+\infty\right[$ are estimated from data. In , it was shown that, for suitable choices of $w_{j}$, the adaptive Lasso produces (asymptotically) unbiased estimates of $b$. One of the first methods to alleviate the $\sigma$-dependency of the Lasso has been the Sqrt-Lasso . The Sqrt-Lasso problem is based on the formulation

This optimization problem can be cast as second order cone program (SOCP) . The modification of the objective function can be interpreted as an (implicit) scaling of the Lasso objective function by an estimate $\|Xb-z\|_{2}/\sqrt{n}$ of $\sigma$ , leading to

In , it was shown that the tuning parameter $\lambda$ does not depend on $\sigma$ in Sqrt-Lasso.

Alternative approaches rely on the idea of simultaneously and explicitly estimating $b$ and $\sigma$ from the data. The scaled Lasso , a robust hybrid of ridge and Lasso regression , and the TREX are important instances. In the following, we will show that these estimators are based on perspective functions under the unifying statistical framework of concomitant estimation. We will introduce a novel family of estimators and show how the corresponding optimization problems can be solved using proximal algorithms. In particular, we will derive novel proximal algorithms for solving both the standard TREX and a novel generalized version of the TREX which includes the Sqrt-Lasso as special case.

2 Penalized concomitant M-estimators

In statistics, the task of simultaneously estimating a regression vector $b$ and an additional model parameter is referred to as concomitant estimation. In , Huber introduced a generic method for formulating “maximum likelihood-type” estimators (or M-estimators) with a concomitant parameter from a convex criterion. Using our perspective function framework, we can extend this framework and introduce the class of penalized concomitant M-estimators defined through the convex optimization problem

3 Proximal algorithms for the TREX

The TREX extends Sqrt-Lasso and scaled Lasso by taking into account the unknown noise distribution of $e$. Recalling that a theoretically desirable tuning parameter for the Lasso is $\lambda\propto\sigma\|X^{\top}e\|_{\infty}$, the TREX scales the Lasso objective by an estimate of this quantity, namely,

The parameter $\alpha>0$ can be set to a constant value ($\alpha=1/2$ being the default choice). In , promising statistical results were reported where an approximate version of the TREX, with no tuning of $\alpha$, has been shown to be a valid alternative to the Lasso. A major technical challenge in the TREX formulation is the non-convexity of the optimization problem. In , this difficulty is overcome by showing that the TREX problem, although non-convex, can be solved by observing that problem (4.8) can be equivalently expressed as finding the best solution to $2p$ convex problems of the form

and the corresponding TREX subproblem is to

Then $f_{j}=g_{j}\circ M_{j}$. Upon setting $h=\|\cdot\|_{1}$, we see that (4.11) is of the form

and where $t$ is the unique solution in $\left]0,+\infty\right[$ to the depressed cubic equation

3.2 Proximal operators for generalized TREX estimators

Thus far, we have shown that the data-fitting function in the TREX subproblem (4.9) is a special case of (2.25). However, the full potential of (2.25) is revealed by taking a general $q\in\left]1,{+\infty}\right[$, leading to the composite perspective function

This function is the data fitting term of a generalized TREX subproblem for the corresponding global generalized TREX objective

we arrive at $f_{j,q}=g_{j,q}\circ M_{j}$. Setting $h=\|\cdot\|_{1}$ the corresponding problem is to

3.3 Douglas-Rachford for generalized TREX subproblems

4 Numerical illustrations

We illustrate the convergence behavior of the Douglas-Rachford algorithm for TREX problems and the statistical performance of generalized TREX estimators using numerical experiments. All presented algorithms and experimental evaluations are implemented in MATLAB and are available at http://github.com/muellsen/TREX. All algorithms are run in MATLAB 2015a on a MacBook Pro with 2.8 GHz Intel Core i7 and 16 GB 1600 MHz DDR3 memory.

We first examine the scaling behavior of the Douglas-Rachford scheme for the TREX subproblem on linear regression tasks. We simulate synthetic data according to the linear model (4.1) with $m=20$ nonzero variables, regression vector $b^{*}=[-1,1,-1,\ldots,0_{p-m}^{\top}]^{\top}$, and feature vectors $X_{i:}\sim N(0,\Sigma)$ with $\Sigma_{ii}=1$ and $\Sigma_{ij}=0.3$, and Gaussian noise $\varepsilon_{i}\sim N(0,\sigma^{2})$ with $\sigma=1$. Each column $X_{:j}$ is normalized to have norm $\sqrt{n}$. We fix the sample size $n=200$ and consider the dimension $p\in\{20,50,100,200,500,1000,2000\}$. We solve one standard TREX subproblem (for $s\in\{-1,1\}$, $X_{:1}$, $\alpha=0.5$) over $d=20$ random realizations of $X$ and $e$. For the TREX subproblem we consider the proximal Douglas-Rachford algorithm 4.29 with parameters $\mu_{k}\equiv 1.95$ and $\gamma=70$. We declare that the Douglas-Rachford algorithm has converged at iteration $K$ if $\min\{\|b_{K+1}-b_{K}\|,\|y_{K+1}-y_{K}\|\}\leqslant 10^{-10}$, resulting in the final estimate $b_{K}$.

In practice, the Douglas-Rachford algorithm for the TREX subproblem can be enhanced by an online sign selection rule (DR-Sel). When a TREX subproblem for fixed $X_{:j}$ is considered, we can solve the problem for $s\in\{-1,1\}$ concurrently for a small number $k_{0}$ of iterations (standard setting $k_{0}=50$) and select the signed optimization problem with best progress in terms of objective function value.

We compare the run time scaling and solution quality of Douglas-Rachford and DR-Sel with those of the state-of-the-art Splitting Conic Solver (SCS). SCS is a general-purpose first-order proximal method that provides numerical solutions to several standard classes of optimization problems, including SOCPs and Semidefinite Programs (SDPs). We use SCS in indirect mode to solve the SOCP formulation of the TREX subproblem with convergence tolerance $10^{-4}$.

The run time scaling results are shown in Figure 1. We emphasize that the scaling experiments are not meant to measure absolute algorithmic performance but rather efficiency with respect to optimization formulations that are subsequently solved by proximal algorithms. We observe that SCS with the SOCP formulation of TREX compares favorably with Douglas-Rachford and DR-Sel in low dimensions while, for $p>200$, both Douglas-Rachford variants perform better. DR-Sel outperforms Douglas-Rachford by a factor of $2$ to $4$ and always selects the correct signed subproblem (data not shown). The TREX solutions found by SCS and Douglas-Rachford are close in terms of $\|b^{(DR)}-b^{(SCS)}\|$, with DR typically reaching slightly lower function values than SCS. Values for the first 40 dimensions of a typical solution $b_{K}$ in $p=2000$ dimensions are shown in Figure 1 (right panels).

4.2 Behavior of generalized TREX estimators

We next study the effect of the exponent $q$ on the statistical behavior of the generalized TREX estimator. We use the synthetic setting outlined in to study the phase transition behavior of the different generalized TREX estimators. We generate data from the linear model (4.1) with $p=64$ and $m=\lceil 0.4p^{3/4}\rceil$ nonzero variables, regression vector $b^{*}=[-1,1,-1,\ldots,0_{p-m}^{\top}]^{\top}$, and feature vectors $X_{i:}\sim N(0,\Sigma)$ with $\Sigma_{ii}=1$ and $\Sigma_{ij}=0$ and Gaussian noise $e$ with $\sigma=0.5$. Each column $X_{:j}$ is normalized to have norm $\sqrt{n}$. We define the rescaled sample size according to $\theta(n,p,m)=n/(2m\log{(p-m)})$ and consider $\theta(n,p,m)\in\{0.2,0.4,\ldots,1.6\}$. At $\theta(n,p,m)=1$, the probability of exact recovery of the support of $b^{*}$ is $0.5$ for the (Sqrt)-Lasso with oracle regularization parameter . We consider the generalized TREX with different exponents $q\in\{9/8,7/6,3/2,2\}$ and the Sqrt-Lasso as limiting case $q=1$. For all generalized TREX estimators we consider regularization parameters $\alpha\in\{0.1,0.15,\ldots,2\}$. For Sqrt-Lasso we consider the standard regularization path setting outlined in . We solve all generalized TREX problems with the Douglas-Rachford scheme using the previously described parameter and convergence settings. We measure the probability of exact support recovery and Hamming distance to the true support over $d=12$ repetitions. We threshold all “numerical zeros” in the generalized TREX solutions vectors at level $0.05$. For all solutions closest to the true support in terms of Hamming distance, we also calculate estimation error $\|b_{K}-b^{*}\|_{2}^{2}/n$ and prediction error $\|Xb_{K}-Xb^{*}\|_{2}^{2}/n$. Figure 2 shows average performance results across all repetitions.

We observe several interesting phenomena for the family of generalized TREX estimators. In terms of exact recovery, the performance is slightly better than predicted by theory (see gray dashed line in Figure 2 top left panel), with decrease in performance for increasing $q$. This is also consistent with average Hamming distance measurements (top right panel). We observe that generalized TREX oracle solutions (according to the minimum Hamming distance criterion) show best performance in terms of estimation and prediction error for exponents $q\in\{9/8,7/6\}$, followed by $q\in\{3/2,2\}$.

The present numerical experiments highlight the usefulness of the family of generalized TREX estimators for sparse linear regression problems. Further theoretical research is needed to derive asymptotic properties of generalized TREX. A central prerequisite for establishing generalized TREX as statistical estimator is to solve the underlying optimization problem with provable guarantees. We have shown that our perspective function framework along with efficient computation of proximity operators enables this important task in a seamless way.

We thank Dr. Jacob Bien for valuable discussions. The Simons Foundation is acknowledged for partial financial support of this research. The work of P. L. Combettes was also partially supported by the CNRS MASTODONS project under grant 2016TABASCO.

References