Estimating Shape Parameters of Piecewise Linear-Quadratic Problems

TL;DR

This paper introduces a statistical framework for automatically tuning shape parameters of Piecewise Linear-Quadratic penalties, enhancing model flexibility and accuracy in statistical inference, signal processing, and machine learning.

Contribution

It develops a joint optimization approach for shape and primary parameters, utilizing normalization constants to inform parameter tuning, and designs efficient interior point methods for solving these problems.

Findings

Effective joint optimization of shape and primary parameters.

Specialized interior point methods improve computational efficiency.

Framework demonstrated on synthetic and real-world examples.

Abstract

Piecewise Linear-Quadratic (PLQ) penalties are widely used to develop models in statistical inference, signal processing, and machine learning. Common examples of PLQ penalties include least squares, Huber, Vapnik, 1-norm, and their asymmetric generalizations. Properties of these estimators depend on the choice of penalty and its shape parameters, such as degree of asymmetry for the quantile loss, and transition point between linear and quadratic pieces for the Huber function. In this paper, we develop a statistical framework that can help the modeler to automatically tune shape parameters once the shape of the penalty has been chosen. The choice of the parameter is informed by the basic notion that each QS penalty should correspond to a true statistical density. The normalization constant inherent in this requirement helps to inform the optimization over shape parameters, giving a joint optimization problem over these as well as primary parameters of interest. A second contribution is to consider optimization methods for these joint problems. We show that basic first-order methods can be immediately brought to bear, and design specialized extensions of interior point (IP) methods for PLQ problems that can quickly and efficiently solve the joint problem. Synthetic problems and larger-scale practical examples illustrate the potential of the approach.

Tables

Table 1. Table 1: Common PLQ Functions and their Conjugate Representations

Name	Conjugate Representation	Figure
Huber	$h_{\kappa }\left(x\right)={sup}_{u\in [-\kappa ,\kappa ]}\left\{ux-\frac{1}{2}{u}^2\right\}$	Fig. 1(a)
Quantile	$q_{\tau }\left(x\right)={sup}_{u\in [-\tau ,\left(1-\tau \right)]}\left\{ux\right\}$	Fig. 1(b)
Quantile Huber	$h_{\tau ,\kappa }\left(x\right)={sup}_{u\in [-\kappa \tau ,\kappa \left(1-\tau \right)]}\left\{ux-\frac{1}{2}{u}^2\right\}$	Fig. 1(c)
Least squares	$\frac{1}{2}{x}^2={sup}_u\left\{ux-\frac{1}{2}{u}^2\right\}$	Fig. 3(a)
Hinge	$h_{ϵ}\left(x\right)={sup}_{u\in [-\tau ,\left(1-\tau \right)]}\left\{u\left(x-ϵ\right)\right\}$	Fig. 3(b)
Vapnik	${\rho }_{ϵ}\left(x\right)={sup}_{u\in {\left[\mathrm{0,1}\right]}^2}\left\{⟨\left[\begin{array}{c}1\\ -1\end{array}\right]x-\left[\begin{array}{c}ϵ\\ ϵ\end{array}\right],u⟩\right\}$	Fig. 3(c)
Smooth insensitive loss	${\rho }_{ϵ}^h\left(x\right)={sup}_{u\in {\left[\mathrm{0,1}\right]}^2}\left\{⟨\left[\begin{array}{c}1\\ -1\end{array}\right]x-\left[\begin{array}{c}ϵ\\ ϵ\end{array}\right],u⟩-\frac{1}{2}{u}^{T}u\right\}$	Fig. 3(d)
Elastic net		Fig. 3(e)

Table 2. Table 2: Joint inference of the shape and model parameters for the quantile Huber loss. Columns 2 and 3 show results for proposed method. r ​ ( x ) = ‖ x − x t ‖ / ‖ x t ‖ 𝑟 𝑥 norm 𝑥 subscript 𝑥 𝑡 norm subscript 𝑥 𝑡 r(x)=\|x-x_{t}\|/\|x_{t}\| denotes relative error. Column 2 contains τ , κ 𝜏 𝜅 \tau,\kappa estimated using joint optimization (compare to ground truth simulation parameters in column 1) solved by the proposed IP algorithm with KKT tolerance set to 10 − 6 superscript 10 6 10^{-6} . Column 3 shows relative error of the new estimate; compare to Columns 4 and 5, which are relative errors for LS and minimum 1-norm estimates.

$[{\tau }_t,{\kappa }_t]$	$[{\tau }^{\ast },{\kappa }^{\ast }]$	$r\left(x^{\ast }\right)$	$r\left(x_{LS}\right)$	$r\left(x_M\right)$
[0.1,1.0]	[0.09,1.17]	0.14	0.41	0.26
[0.2,1.0]	[0.20,1.07]	0.10	0.16	0.13
[0.5,1.0]	[0.50,0.95]	0.08	0.12	0.09
[0.8,1.0]	[0.81,1.04]	0.09	0.19	0.11
[0.9,1.0]	[0.91,1.17]	0.12	0.38	0.36

Table 3. Table 3: Timing comparison (in seconds) of PALM and IP for the quantile Huber family. Columns 2 and 3 show total run time and number of iterations of PALM and IP; column 3 plots difference in final objective values. IP finds lower objectives, but PALM is faster for larger problems. A large-scale implementation of IP that uses iterative methods in each iteration would help with the scaling issues.

$m$	$n$	PALM: Time(s)/#iter	IP: Time(s)/#iter	$f_{\text{PALM}}^{\ast }-f_{\text{IP}}^{\ast }$
100	50	4.12/96	3.11/15	7.25e-12
500	50	5.24/197	4.86/16	1.38e-08
1000	50	5.97/112	11.68/14	1.52e-08
2000	100	11.41/129	69.85/16	1.30e-08
2000	200	22.72/225	72.89/16	9.61e-08
2000	500	66.30/394	87.85/18	3.86e-08

Table 4. Table 4: Column 3 shows time (s)/iterations for gradient descent on a single instance of the quantile huber problem. Columns 4 and 5 show the time/iterations for PALM for the joint problem over ( x , θ ) 𝑥 𝜃 (x,\theta) and resulting estimates of ( κ , τ ) 𝜅 𝜏 (\kappa,\tau) . Columns 6 and 7 shows the time and estimates obtained by a 5 × 5 5 5 5\times 5 grid search over ( κ , τ ) 𝜅 𝜏 (\kappa,\tau) . Penalty optimization finds good estimates very quickly relative to the baseline. Grid search takes an order of magnitude more time.

$m$	$n$	$x$ only: Time (s)/#iter	Joint: Time(s)/#iter	$[{\kappa }^{\ast },{\tau }^{\ast }]$	Cross-Val Time(s)	$[{\kappa }_{CV},{\tau }_{CV}]$
100	50	24.57/4008	16.98/2198	[0.50,1.28]	585.12	[0.5,1.5]
500	50	4.13/430	7.31/386	[0.53,1.03]	309.96	[0.5,1.0]
1000	50	3.78/202	7.58/224	[0.49,1.00]	349.00	[0.5,1.0]
2000	100	14.22/219	19.79/255	[0.49,1.02]	1026.92	[0.5,1.0]
2000	200	29.20/388	42.67/513	[0.49,1.04]	2315.02	[0.5,1.0]
2000	500	158.32/1025	202.23/1261	[0.52,1.20]	7477.29	[0.5,1.0]

Table 5. Table 5: Runtime comparison for self-tuning Huber RPCA (Joint) v.s. single Huber RPCA with fixed ( κ , σ ) 𝜅 𝜎 (\kappa,\sigma) parameters (Simple). We fix number of iterations to be 200.

Joint: Time (s)/#iter	Simple: Time (s)/#iter
203.11/200	126.72/200
210.74/200	107.24/200

Equations

$$y=Ax+ϵ,ϵ\sim N(0,{\sigma }^{2}I),i=1,\dots ,n,$$

$$p(x,{\sigma }^2|y)\propto p(y|x,\sigma 2)=\prod _{i=1}^n\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp ((y_i-⟨a_i,x⟩{)}^2/2{\sigma }^2)$$

$$\underset{x,{\sigma }^2}{\min }\frac{1}{2{\sigma }^2}\Vert y-Ax{\Vert }^2+n\log ({\sigma }^2).$$

$$\hat{x}=(A^{T}A{)}^{-1}{A}^{T}b.$$

$${\hat{\sigma }}^2=\Vert y-A\hat{x}{\Vert }^2/n,$$

$$\underset{x}{\min }\frac{1}{2}\Vert Ax-y{\Vert }^2$$

$$y_i=⟨a_i,x⟩+{ϵ}_i,{ϵ}_i\sim \exp (-{\rho }_{\tau }(\cdot )),$$

$$\underset{x,\tau \in [0,1]}{\min }{q}_{\tau }(Ax-b)+m\log \left(\frac{1}{\tau }+\frac{1}{1-\tau }\right).$$

$$f^{\ast }(v)=\underset{u}{\sup }{u}^{\mathsf{T}}v-f(u).$$

$$\rho (r)=\underset{u}{\sup }\left\{u^{\mathsf{T}}(Br-\bar{b})-\frac{1}{2}{u}^{\mathsf{T}}Mu:C^{\mathsf{T}}u\le \bar{c}\right\},$$

$$\lim \underset{|x|\to \infty }{\inf }f(x)=\infty .$$

$$p(r;\theta )$$

$$\underset{x}{\min }\sum _{i=1}^m\rho (y_i-⟨a_i,x⟩;\theta )$$

$$\underset{x,\theta \in \mathcal{D}}{\min }\sum _{i=1}^m\rho (y_i-⟨a_i,x⟩;\theta )+m\log [n_c(\theta )],$$

$${\rho }_{\tau }(x)=\{\begin{array}{ll}(1-\tau )x& x\ge 0\\ -\tau x& x\le 0\end{array},p(r;\tau )=\{\begin{array}{ll}\exp (-(1-\tau )x)& x\ge 0\\ \exp (-\tau x)& x\le 0\end{array},n_c(\tau )=\frac{1}{1-\tau }+\frac{1}{\tau }.$$

$$\nabla n_c({\theta }_0)$$

$${\nabla }^2{n}_c({\theta }_0)$$

$$n_c({\theta }_0+h)-n_c({\theta }_0)$$

$$\Rightarrow$$

$$\underset{h\to 0}{\lim }\frac{n_c({\theta }_0+h)-n_c({\theta }_0)}{\Vert h\Vert }=\underset{h\to 0}{\lim }{\int }_{\mathbb{R}}⟨{\nabla }_{\theta }\exp [-\rho (r;\bar{\theta })],\frac{h}{\Vert h\Vert }⟩dr={\int }_{\mathbb{R}}\underset{h\to 0}{\lim }⟨{\nabla }_{\theta }\exp [-\rho (r;\bar{\theta })],\frac{h}{\Vert h\Vert }⟩dr$$

$$={\int }_{\mathbb{R}}⟨{\nabla }_{\theta }\exp [-\rho (r;{\theta }_0)],v⟩dr=⟨{\int }_{\mathbb{R}}{\nabla }_{\theta }\exp [-\rho (r;{\theta }_0)]dr,v⟩$$

$$\nabla n_c({\theta }_0)={\int }_{\mathbb{R}}{\nabla }_{\theta }\exp [-\rho (r;{\theta }_0)]dr.$$

$${\nabla }^2{n}_c({\theta }_0)={\int }_{\mathbb{R}}{\nabla }_{\theta }^2\exp [-\rho (r;{\theta }_0)]dr.$$

$$v(\theta )=\underset{x}{\min }\sum _{i=1}^m\rho (y_i-⟨a_i,x⟩;\theta )+m\log [n_c(\theta )].$$

$$\underset{x,\theta }{\min }H(x,\theta )+r_1(x)+r_2(\theta ),$$

$$H(x,\theta )$$

$${\delta }_{\mathcal{D}}(\theta )=\{\begin{array}{ll}0& \text{mbox}if\theta \in \mathcal{D}\\ \infty & \text{mbox}if\theta \notin \mathcal{D}\end{array}.$$

$${\text{prox}}_{\frac{1}{d_k}{r}_2}(\varphi )=\arg \underset{\theta \in \mathcal{D}}{\min }\frac{1}{2d_k}\Vert \theta -\varphi {\Vert }^2+\log [n_c(\theta )].$$

$$\underset{x,\theta \in \mathcal{D}}{\min }\rho (x;\theta )+\log [n_c(\theta )].$$

$$\bar{b}=G^{\mathsf{T}}\theta +b,\bar{c}=H^{\mathsf{T}}\theta +c,\mathcal{D}=\{\theta :S^{\mathsf{T}}\theta \le s\}.$$

$$\underset{x,S^{\mathsf{T}}\theta \le s}{\min }$$

$$\underset{x,\tau \in {\mathcal{C}}_1}{\min }$$

Full text

In this paper, we develop a statistical framework that can help the modeler to automatically tune shape parameters once the shape of the penalty has been chosen. The choice of the parameter is informed by the basic notion that each QS penalty should correspond to a true statistical density. The normalization constant inherent in this requirement helps to inform the optimization over shape parameters, giving a joint optimization problem over these as well as primary parameters of interest. A second contribution is to consider optimization methods for these joint problems. We show that basic first-order methods can be immediately brought to bear, and design specialized extensions of interior point (IP) methods for PLQ problems that can quickly and efficiently solve the joint problem. Synthetic problems and larger-scale practical examples illustrate the potential of the approach.

1 Introduction

Piecewise Linear-Quadratic (PLQ) functions and their superclass, Quadratic Support functions [1] form a rich class of penalties that are used for a range of applications. The choice of penalty plays a key role in the properties of the final estimate. Using the Huber penalty (Figure 1(a)) rather than a quadratic to measure data misfit is a typical approach to make an estimate robust to outliers [16]. The 1-norm applied directly to parameters makes the final solution sparse, and is used in Lasso regression [27] and compressed sensing [14]. The asymmetric quantile loss (see Figure 1(b)) is used for many regression applications [19].

The correspondence of the shape of a PLQ function to its role in estimation is well-understood. For example, it is the linear tails of the Huber penalty that limit the effects of large residuals on the final estimate when used as a misfit, and it is the nonsmooth behavior at the origin of the 1-norm that promotes sparse solutions when used as a regularizer. Likewise, the asymmetry of the quantile loss that make it useful for financial applications, where loss and gain are treated asymmetrically.

Here, we consider the choice of parameters that fully specify this shape. At what value should the Huber transition between quadratic and linear? How asymmetric do we need the quantile loss to be? These questions are easily answered once we have residuals and estimates in hand, but they seem as difficult a priori as any other parameter selection problem.

Unknown parameters in are usually estimated using cross-validation or grid search. Standard methods require multiple solutions of any given learning problem, where a held-out dataset is used to evaluate each configuration [5]. More recently, Bayesian optimization [26, 17, 9, 18] and random search [8, 21] have come to the forefront as two popular techniques used to obtain meta-parameters in a very wide range of contexts.

While these techniques can be used to for the problems we consider, as well as a much more general problem class, their implementation is more expensive than solving a single problem instance — cross-validation, random search and Bayesian optimization require many instance evaluations.

Here we take a different tack, and build on the statistical perspective developed in [1] to develop a simple modification of the PLQ estimation problem to infer shape parameters simultaneously with variables of interest. The most relevant works related to this paper focus on the relation between the quantile penalty and the asymmetric Laplace distribution(ALD) [31, 28, 7]. In particular, [7] jointly estimate the model and the shape parameters for quantile penalty, and [28] infer the joint posterior distribution of these parameters.

To explain the general idea, we look at two simple examples.

1.1 Classic variance estimation in linear regression.

As a warm-up, we consider the classic problem of estimating regression variables jointly with noise variance parameter in a linear Gaussian model. Take the basic statistical model

[TABLE]

we can find the maximum likelihood estimator for $(x,\sigma^{2})$ by maximizing

[TABLE]

or equivalently by minimizing the negative log of this expression:

[TABLE]

In this simple case, the problem separates, and regardless of $\sigma^{2}$ we have

[TABLE]

Minimizing (2) with respect to $\sigma^{2}$, we also get

[TABLE]

the average of empirical residuals squared. The main point to here is that the likelihood arising from the statistical model (1) enables us to estimate $\sigma^{2}$. While the two variables separate in this case, in general the picture is more interesting, as we see in the next example.

1.2 Simultaneous quantile estimation and regression.

Consider Figure 2, where linear observations have been contaminated with asymmetric errors, i.e. they are more likely to be positive than negative. The data generating mechanism is shown in black. The linear regression model for the data $\{y_{i},a_{i}\}$ is equivalent to the least square problem,

[TABLE]

where $x$ contains slope and intercept. As expected, classic regression fails to capture the true mechanism because the errors are asymmetric, see blue dash in Figure 2.

In this case, we can actually recover the ‘true’ line by assuming the errors arise from an asymmetric generalization of the Laplace distribution:

[TABLE]

where $\tau\in[0,1]$ controlling the relative slopes to the left and right of [math], see Figure 1(b). Solving for a joint maximum likelihood estimator for $(x,\tau)$ , we get the fit in red dash in Figure 2. The estimator is derived by considering the density corresponding to (3), and has the simple expression

[TABLE]

where the last term arises from considering the normalization constant required to make $\exp(-\rho_{\tau}(\cdot))$ a true statistical density, and is derived in Section 3, where the general approach is developed. In this simple example, the right solution is found by solving the single problem (4) once, rather than considering multiple problem instances as we expect from classic approaches. The optimization problem itself is nonsmooth and nonconvex, but as we will see straightforward. The main thrust of the paper is to get a systematic way of formulating estimators such as (4), understand their properties, and consider different algorithms for solving these problems.

1.3 Contributions and Roadmap

Building on the toolkit of PLQ penalties and corresponding densities developed in [1], we use the bridge between penalties and densities to develop an extended likelihood formulation over both shape parameters $\theta$ and regression variables $x$. The key idea is encoded in the normalization constant, that arises from the statistical interpretation and essentially acts as a balance against the data to ensure the final shape parameters still yield a statistically valid model. We consider properties of the extended likelihood formulations over $(x,\theta)$. In many cases, we show that standard techniques (including proximal alternating minimization and variable projection) can be applied. We also build on the interior point method of [1] to develop a new extended interior point approach tailored to these joint estimators.

The paper proceeds as follows. In Section 2 we review PLQ functions, along with their conjugate representations, and examples. Conjugate representations are essential for solving PLQ models using interior point methods. In Section 3, we focus on the statistical interpretation, and derive the normalization constant needed to create joint estimators such as (4). In Section 4, we consider optimization methods that can be applied to the new class of estimators, and derive a new customized interior point method for the class. In Section 5 we provide synthetic verifications that show we can recover shape parameters and regression variables in simple examples. Section 6 contains some applications of these ideas to real datasets, particularly focusing on self-tuning penalties for robust PCA.

2 Quadratic Support Functions and Their Conjugate Representations

Piecewise linear-quadratic (PLQ) functions [25, Definition 10.20] are convex functions whose domain can be represented as the union of finitely many polyhedral sets, relative to each of which the function can be written as a convex quadratic. These functions have a convenient representation using their convex conjugates.

Definition 1** (Convex conjugate)**

The convex conjugate of a function $f$ is given by

[TABLE]

A PLQ function is defined as the conjugate to a quadratic over a polyhedral set [25, Example 11.18]. Examples of common penalties used in statistical modeling, machine learning, and inverse problems can all be written using simple conjugates, as shown in Table 1. Quadratic support (QS) functions [1] generalize PLQ functions generalize this class by removing the polyhedral set restriction [1].

PLQ functions are closed under sums and affine compositions [1], and we can explicitly define PLQ functions using their conjugate representation.

Definition 2** (PLQ Functions and penalties.)**

A piecewise linear-quadratic (PLQ) function is given by

[TABLE]

where $M\succeq 0$, $B$ is an injection, and $U:=\{u:C^{\mathsf{T}}u\leq\bar{c}\}$ is a polyhedral.

When $0\in U$, we may call $\rho$ a penalty to emphasize that it is non-negative.

The notion of coercivity is essential to the statistical interpretation of PLQs in the next section.

Definition 3** (Coercivity)**

A function $f(x)$ is coercive if

[TABLE]

Geometric conditions for coercivity of PLQs are given in [1, Theorem 10]. In the next section, we use PLQs and their properties to build a bridge to densities, which enable the joint likelihood approach.

3 Statistical Model and Properties of Joint Objective

To formulate a learning problem over spaces of penalties, we first review the relationship between penalties and corresponding residual distributions. We then use this relationship to develop a joint formulation for model and shape parameter inference, and characterize properties of the resulting objective function.

The relationship between penalties and associated densities can be made precise. Every coercive PLQ penalty is integrable. Given a coercive penalty $\rho(r;\theta)$ whose shape is parametrized by $\theta$, we define an associated density

[TABLE]

The term $n_{c}(\theta)$ is a normalization constant that ensures that $\rho(r;\theta)$ in (6) is a true statistical density, i.e. it integrates to $1$. These observations yield a simple recipe for extending a negative log likelihood in $x$ (based on optimizing a PLQ penalty) to inform shape parameters $\theta$. Specifically, given an optimization problem of the form

[TABLE]

we consider the extended problem

[TABLE]

where $\theta$ may be restricted to a domain $\mathcal{D}$. For example, in the quantile regression problem in Section 1.2, we have $\mathcal{D}=[0,1]$ and

[TABLE]

The $\log(n_{c})$ term acts as a barrier (see Figure 2), pushing $\tau$ away from the boundary points [math] and $1$ into the interior $(0,1)$, and favoring $\tau=0.5$.

Interestingly, the log of the normalization constant for the quantile regression problem is smooth and strongly convex, but its gradient is not Lipschitz continuous. In the remainder of this section, we characterize theoretical properties of the general objective (7).

Assumption 4

To ensure the validity of the statistical viewpoint, we require $\rho$ to satisfy:

$\rho(r;\theta)\geq 0$* is a PLQ penalty for every $\theta\in\mathcal{D}$, so we have non-negativity.* 2. 2.

$\rho$* is coercive for any $\theta\in\mathcal{D}$, so we have integrability.* 3. 3.

For any $\theta_{0}\in\mathcal{D}$, $\rho(r;\theta)$ is $C^{2}$ around $\theta_{0}$ for almost every $r\in\mathbb{R}$ (smoothness in $\theta$)

Under these assumptions, we can obtain formulas for the first and second derivatives of $n_{c}(\theta)$. {theo}[smoothness of $n_{c}(\theta)$]

For $n_{c}(\theta)$in (6), suppose Assumption 4 holds and for $\theta_{0}\in\mathcal{D}$, there exist functions $g_{k}(r)$, $k=1,2$, such that,

for any unit vector $v$, $|\left\langle\nabla_{\theta}\exp[-\rho(r;\theta)],v\right\rangle|\leq g_{1}(r)$ for any $\theta$ around $\theta_{0}$, 2. 2.

for any unit vector $v$, $\left|\left\langle\nabla_{\theta}^{2}\exp[-\rho(r;\theta)]v,v\right\rangle\right|\leq g_{2}(r)$ for any $\theta$ around $\theta_{0}$, 3. 3.

$\int_{\mathbb{R}}g_{k}(r)\,dr<\infty$, $k=1,2$.

then $n_{c}(\theta)$ is $C^{2}$ continuous around $\theta_{0}$ and,

[TABLE]

The proof is elementary, and is included below for completeness.

Proof 3.5**.**

From Assumption 4, we know that for any $\theta_{0}\in\mathcal{D}$, $\nabla_{\theta}\exp[\rho(r;\theta_{0})]$ and $\nabla_{\theta}^{2}\exp[\rho(r;\theta_{0})]$ exist for almost every $r\in\mathbb{R}$. For any $h$ such that $\|h\|$ is small enough to make $\theta_{0}+h$ stay in the neighborhood of $\theta_{0}$. By applying mean value theorem, we have,

[TABLE]

where $\bar{\theta}$ lie in segment with the end points $\theta_{0}$ and $\theta_{0}+h$. The first and third assumptions allow us to apply the dominant convergence theorem and get,

[TABLE]

where we set $h=\alpha v$ and let $\alpha\rightarrow 0^{+}$ and keep $v$ fix as an unit vector. From the definition of the gradient we know that,

[TABLE]

Following the same steps we could also show $\nabla^{2}c(\theta_{0})$ exists and satisfies,

[TABLE]

The derivative formulas (8) are used for first- and second-order methods to infer $x$ and $\theta$. The parametrization conditions in $\theta$ are satisfied by all piecewise linear-quadratic (PLQ) penalties. The theorem can also be applied to densities that are not log-concave. For instance, the Student’s $t$ density and associated penalty satisfy all assumptions of Theorem 4.

In the quantile case (4), the term $\log[n_{c}(\theta)]$ is convex. We characterize sufficient conditions for convexity of $\log[n_{c}(\theta)]$ for a general class of penalties $\rho$. The results can be derived using [11, Chapter 3.5]. {theo}[convexity of $\log[n_{c}(\theta)]$]

Let $n_{c}(\theta)$ be defined as in Theorem 4, and suppose Assumption 4 holds. We have the following results:

If $\rho(r;\theta)$ is jointly convex in $r$ and $\theta$, then $\log[n_{c}(\theta)]$ is a concave function of $\theta$. 2. 2.

If $\rho(r;\theta)$ is concave with respect to $\theta$ for every $r$, then $\log[n_{c}(\theta)]$ is a convex function.

Theorems 4 and 4 tell an interesting story. The log-normalization constant $\log[n_{c}(\theta)]$ is nearly always smooth; even when the loss $\rho$ is nonsmooth in $x$. The inference problem (7) is never guaranteed to be jointly convex in $(x,\theta)$: if $\rho(x;\theta)$ is jointly convex, then $\log[n_{c}(\theta)]$ will be concave. This is intuitive, as we are attempting to learn both the model and error structure at the same time. Objective (7) is generally nonsmooth and nonconvex. In the next section, we show how to optimize it using first and second order methods.

3.1 Level sets of the objective function and maximum likelihood estimate

Although (7) is non-convex, in many cases we can still find the global minimum in $\theta$. To illustrate, we generate the samples $\epsilon_{i}$ from distribution defined by quantile Huber function with $\kappa=1$ and $\tau=0.05$, select a set of weights $x$, generate data $y_{i}=\left\langle a_{i},x\right\rangle+\epsilon_{i}$, and plot the so called value function for (7):

[TABLE]

Evaluating $v(\theta)$ requires solving a smooth convex problem in $x$, and thus problem (7) reduces to minimizing (9) in $\theta$. Results for the simple 2D case are shown in Figure 4. $v(\theta)$ is non-convex (note the level sets), but there is a unique minimum that is close to the true parameters, and was found in every run by a local search. In the next section, we design methods that are far more efficient than optimizing $v(\theta)$.

4 Optimization Methods for Joint PLQ and Shape Parameter Estimation.

In this section, we discuss algorithms for (7), and develop a new interior point method building on the PLQ optimization strategy of [1] When $\rho$ is smooth in $x$ and $\theta$, we show how to apply Proximal Alternating Linearized Minimization (PALM) [10] and Proximal Alternating Minimization (PAM) [6]. We also discuss variable projection (VP) type algorithms [4, 3]. These development is straightforward, but the log normalization constant $\log[n_{c}(\theta)]$ must be treated carefully, as its gradient does not have a global Lipschitz constant.

Requiring smoothness in $\rho$ is restrictive, and the quantile regression example (LABEL:eq:quantile) has a nonsmooth $\rho$. The quantile penalty is not smooth, but it is PLQ. We use conjugate representations of PLQs in Section 2 to extend interior point methods developed in [1] to problem (7) in Section 4.2.

4.1 PALM, PAM, and Variable Projection

Several algorithms from existing literature can be brought to bear on problem (7) when it is of the form

[TABLE]

under some mild assumptions on $H$. Here we briefly review PALM [10], and PAM [6] algorithms, as well as similar algorithms that partially minimize the objective in $\theta$, known as partial minimization or variable projection [4, 3].

The PALM and PAM algorithms [10] can be used to minimize (10) when $H$ is $C^{1}$, with globally Lipschitz partial gradients, while the functions $r_{1}$ and $r_{2}$ are proper lower semicontinuous; in particular they are not required to be convex, finite-valued, or smooth. The term $r_{1}(x)$ is useful if we need simple regularization and constraints on $x$, such as sparsity or non-negativity. Even though $\log[n_{c}(\theta)]$ is smooth (see Theorem 4), it must be relegated to $r_{2}(\theta)$, since otherwise it violates the Lipschitz assumptions on the partial gradients of $H$. Therefore, to apply PALM or PAM to (7), we take

[TABLE]

Here $\delta_{\mathcal{D}}$ is the indicator function for the set $\mathcal{D}$,

[TABLE]

The PALM algorithm is detailed in Algorithm 1. The steps $c_{k}$ and $d_{k}$ are obtained from Lipschitz constants of the (partial) gradients of $H$. The PAM algorithm is given in Algorithm 2, with the same step sizes to facilitate an easier comparison.

PAM essentially requires the ability to partially minimize a quadratically regularized problem $H$ with respect to both $x$ and $\theta$. In many cases, one of these variables may be easier to solve for than the other; in Example 3.1, $\theta$ has dimension $2$. Variable projection algorithms [4, 3] exploit the ability to fully minimize in one variable while applying an iterative approach in the other. Algorithm 3 shows an example where we have chosen to optimize out $\theta$ for every iteration in $x.$

Detail: The prox operator of $\log[n_{c}(\theta)]$ is not available in closed form for any examples of interest. Instead, it can be efficiently computed using the results of Theorem 4:

[TABLE]

In all examples of interest, $\theta$ is low dimensional; and we compute (12) using Newton’s method or an interior point method. This requires $\nabla\log[n_{c}(\theta)]$ and $\nabla^{2}\log[n_{c}(\theta)]$, which are calculated numerically using formulas (8).

The PALM algorithm works well for large-scale shape inference problems with smooth $\rho$. We use it for the self-tuning RPCA experiments in Section 6.

4.2 Interior Point method for nonsmoothly coupled nonconvex joint QS inference

In this section, we use conjugate representations of PLQ penalties to develop an interior point method for the extended inference problem (7):

[TABLE]

The approach converges superlinearly in practice, but each iteration requires solving a linear system. We solve these systems directly, and scaling issues are explored in Table 3. In practice, large-scale linear systems are solved iteratively, often using pre-conditioners [23]; we leave these developments for future work.

To optimize PLQ penalties parametrized by $\theta$, we allow $\bar{b}$ and $\bar{c}$ to be affine functions of $\theta$, and assume $\mathcal{D}$ is also polyhedral:

[TABLE]

We then solve a saddle point for primal variables $x$, conjugate variables $u$, and shape parameters $\theta$:

[TABLE]

For example, the self-tuning quantile penalty (4) is written as

[TABLE]

Interior point (IP) methods apply damped Newton to a relaxation of the optimality conditions of (13), see [20, 22, 30]. Let $F(z)$ denote of optimality conditions for problem (13), with $z=(x,u,\theta,\lambda)$, where $\lambda$ are dual variables for inequality constraints $C^{\mathsf{T}}u\leq H^{\mathsf{T}}\theta+c$ and $S^{\mathsf{T}}\theta\leq s$, and let $F_{\mu}(z)$ denote the relaxed system obtained by using a logarithmic barrier with parameter $\mu$. The IP method is summarized in Algorithm 4.

We introduce a log-barrier for the conjugate variables $u$ and shape parameters $\theta$ as follows.

[TABLE]

As $\mu\downarrow 0$, the barriers approach true indicator functions for the $\mathcal{D}$ and $U$. The parameter $\mu$ is decreased to a specified tolerance as the optimization proceeds. For fixed $\mu$, the approximate subproblem with fixed $\mu$ can itself be written as a saddle point system:

[TABLE]

We introduce dual variable $q$ as a function of the slacks $d$:

[TABLE]

and form the KKT system,

[TABLE]

The Jacobian matrix $F_{\mu}^{(1)}$ of the system is given by

[TABLE]

Algorithm 4 is a damped Newton iteration to find the stationary point of $F_{\mu}$. At each iteration, we shrink the $\mu$ to be a fraction of the current average complementarity conditions, just as in the implementation used in [1, 2].

We can apply block Gaussian elimination to obtain conditions that make the algorithm implementable. We state the result as a siimple theorem.

{theo}

[IP implementability] Suppose that the PLQ penalty is nondegenerate, that is $\mbox{null}(C)\cap\mbox{null}(M)=\{0\}$, and that the linear model $A$ is full-rank. Then the interior point iteration

[TABLE]

is implementable when a certain square symmetric system is invertible:

[TABLE]

$T_{3}$ is a symmetric square matrix with dimension equal to the length of the parameter vector $\theta$.

Proof 4.6**.**

Implementing the IP iteration is equivalent to applying block-Gaussian elimination to the system (16). The set of operations is given by

[TABLE]

where we define

[TABLE]

The invertibility of $D$ is enforced by Algorithm 4, which keeps slack variables $d$ strictly positive. Invertibility of $T_{1}$ is guaranteed by the nondegeneracy hypothesis for the PLQ, see [1, Theorem 14]. Since $B$ must be injective (see Definition 5), the full rank condition on $A$ guarantees that $T_{2}$ is invertible. The row operations yield an upper triangular matrix of form

[TABLE]

which is invertible if and only if $T_{3}$ in (17) is invertible, given the other hypotheses.

In the next section, we present synthetic examples and that show the applicability of the approach for inferring simple shape parameters.

5 Synthetic Data Experiments

We illustrate the approach using a linear regression example. Consider a regression problem over the quantile Huber family (Figure 1) for a data set contaminated by asymmetric errors and outliers. The $\tau$ parameter controls slope of the penalty, while $\kappa$ is the robustness threshold. We want to fit the regression model $x$ as well as obtain the correct parameters $\tau$ and $\kappa$.

When $\kappa>0$ in quantile Huber, $\rho(x;\theta)$ is smooth, and we can use the PALM algorithm from Section 4.1. The quantile Huber penalty is QS, so we can also apply the proposed IP method from Section 4.2. We use both and compare their performance.

The primal form for the quantile Huber penalty is

[TABLE]

We must choose a parametrization in terms of $\theta$. One option is to take $\theta=[\tau,\kappa]^{\mathsf{T}}$. But this parametrization violates assumptions of both first- and second-order approaches in Section 4. Indeed, $\nabla_{\theta}\rho(r;\theta)$ would not have a global Lipschitz constant, which is not valid for PALM. Neither could we write (13) using affine functions of $\theta$. Looking carefully at (19), we instead choose $\theta_{1}=\tau\kappa,\theta_{2}=\tau(1-\kappa)$. These new parameters must be non-negative.

The primal objective can be written as

[TABLE]

where $A\in\mathbb{R}^{m\times n}$ is the design matrix (always random Gaussian), $x\in\mathbb{R}^{n}$ is the model parameter vector, and $y\in\mathbb{R}^{m}$ is the observed data vector contaminated by outliers. From Theorem 4, $n_{c}(\theta)$ is $\mathcal{C}^{2}$ smooth. From Theorem 4, the objective in $\theta$ is the sum of a concave term $\rho(Ax-y;\theta)$ and a convex term $m\log[n_{c}(\theta)]$. The joint problem in $(x,\theta)$ is nonconvex, but both first- and second-order methods from Section 4 can be applied.

5.1 Shape-Optimized Quantile Huber vs. Least Square and LAD

We generate synthetic data with $m=1000$ and $n=50$. The measurement errors are sampled from quantile Huber distributions, to verify that the approach is able to recover ‘ground truth’ values for $(\tau,\kappa)$ parameters. We denote ground truth parameters as $x_{t}$, $\tau_{t}$, $\kappa_{t}$, while $x^{*}$, $\tau^{*}$, $\kappa^{*}$ are the solutions obtained by solving (7). We provide two reference solutions: $x_{LS}$ is the least square solution, and $x_{M}$ is the solution obtained by solving $\|Ax-b\|_{1}$. For each $\kappa$ and $\tau$ setting, we run the simulation 10 times, and show the average of the results in Table 2. Results shown are obtained by the IP method.

The maximum likelihood formulation correctly recovers the shape parameters $(\theta,\tau)$. Moreover, the solution $x^{*}$ obtained from the self-tuned regression is always better compared to reference solutions, and the improvement increases as measurement errors become more biased ($\tau$ close to [math] or to $1$).

5.2 PALM vs. IP for Penalty Optimization

We also compared the performance of PALM and IP over $m\in\{100,500,1000,2000\}$ and $n\in\{50,100,200,500\}$, for both iterations and run time. The results are shown in Figure 5 and Table 3. From Table 3, IP converges within 20 iterations independently of problem dimension. However, the total run time exceeds PALM as the problem size increases. This issue can be mitigated by implementing iterative solvers with pre-conditioners; which would be particularly valuable for nonsmooth $\rho$, where PALM cannot be applied. However, this effort is outside the scope of the paper.

5.3 Penalty Optimization vs. Grid Search

Finally, we compared the run time and shape parameter estimates for the quantile Huber family with a grid search method. We construct a $5\times 5$ grid over $\tau\in[0.1,0.2,0.5,0.8,0.9]$ and $\kappa\in[0.1,0.5,1.0,1.5,1.9]$, split the data into training and validation sets (80% and 20%) and pick the parameter that performs the best over validation data. We also provide the run time for a single instance evaluation with fixed $\theta$ (using gradient descent) as a baseline comparison.

Results are shown in Table 4. From the table we can see that both methods find good estimates of the shape parameters. However, grid search clearly suffers form the ‘curse of dimensionality’ even at dimensionality 2, as it takes approximately 20 times longer to complete.

6 Real Data Example

In this section, we consider large-scale examples, developing approaches for using ‘self-tuning’ penalties for robust principal component analysis (RPCA). RPCA has applications to alignment of occluded images [24], scene triangulation [32], model selection [13], face recognition [29] and document indexing [12]. We develop a self-tuning background separation approach. Given a sequence of images (Campus dataset)111Publicly available at http://vis-www.cs.umass.edu/~narayana/castanza/I2Rdataset/ , our goal is to separate the moving objects from the background. We pick $202$ images from the data set, convert them to grey scale and reshape them as column vectors of matrix $Y\in\mathbb{R}^{20480\times 202}$. We model the data $Y$ as the sum of low rank component $L$ and sparse noise $S$; we expect moving objects to be captured by $S$.

We take advantage of the fact that RPCA is equivalent to regularized Huber regression [15]:

[TABLE]

We introduce parameter $\sigma$ for the Huber penalty to automatically estimate the right scale of the residual. The joint $(\kappa,\sigma)$ parametrization is given by

[TABLE]

We then model the low rank component as $L=U^{\mathsf{T}}V$, where $U\in\mathbb{R}^{k\times m}$ and $V\in\mathbb{R}^{k\times n}$. The resulting self-tuning RPCA formulation is given by

[TABLE]

We solve the problem with PALM, obtaining the result in Figure 6(b). As the optimization proceeds, $\kappa$ and $\sigma$ decrease to [math] with a fixed ratio $\alpha=\kappa/\sigma$. The self-tuning Huber becomes the scaled 1-norm, recovering the original RPCA formulation [12]. The result in Figure 6(b) is an improvement over the result with initial $(\kappa,\sigma)$ values shown in Figure 6(a).

As in the Section 5, we compared the run time of solving the joint (self-tuning) problem with solving a single instance of RPCA. The results are shown in Table 5. In this case, the extended self-tuning problem twice the run time compared to a single instance. Cross-validation, grid-search and black-box optimization typically require multiple solutions of the original optimization problem to find good parameter values.

7 Discussion

In this paper we developed a simple approach that extends a regression problem using PLQ penalties to also infer unknown shape constraints. We used the statistical interpretation of the PLQ penalty to introduce an additional term that relates shape constraints to a normalization constant.

Many existing algorithms can be brought to bear on the extended problem. For smooth penalties, the PALM algorithm was quite useful in both synthetic and real examples, particularly for large-scale problems. We also developed an interior point method for the PLQ class, that uses conjugate representations of such penalties to solve an augmented saddle point system.

The approach offers several interesting avenues for future research. The nonconvex coupling between convex PLQ penalties and their shape parameters is interesting, and finding conditions when the approach is guaranteed to work in general is an open question. The maximum likelihood criterion itself is just one approach, and may have limitations as the number of shape parameters grows with respect to the data. Characterizing these limitations and trying alternatives (such as marginal likelihood) is also left to future work.

Références

• [1]

Aleksandr Y. Aravkin, James V. Burke, and Gianluigi Pillonetto.

Sparse/robust estimation and kalman smoothing with nonsmooth log-concave densities: Modeling, computation, and theory.

Journal of Machine Learning Research, 14:2689–2728, 2013.

• [2]

Aleksandr Y Aravkin, James V Burke, and Gianluigi Pillonetto.

Generalized system identification with stable spline kernels.

SIAM Journal on Scientific Computing, 40(5):B1419–B1443, 2018.

• [3]

Aleksandr Y Aravkin, Dmitriy Drusvyatskiy, and Tristan van Leeuwen.

Efficient quadratic penalization through the partial minimization technique.

IEEE Transactions on Automatic Control, 63(7):2131–2138, 2017.

• [4]

Aleksandr Y Aravkin and Tristan Van Leeuwen.

Estimating nuisance parameters in inverse problems.

Inverse Problems, 28(11):115016, 2012.

• [5]

Sylvain Arlot, Alain Celisse, et al.

A survey of cross-validation procedures for model selection.

Statistics surveys, 4:40–79, 2010.

• [6]

Hédy Attouch, Jérôme Bolte, Patrick Redont, and Antoine Soubeyran.

Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the kurdyka-łojasiewicz inequality.

Mathematics of operations research, 35(2):438–457, 2010.

• [7]

Anil K Bera, Antonio F Galvao, Gabriel V Montes-Rojas, and Sung Y Park.

Asymmetric laplace regression: Maximum likelihood, maximum entropy and quantile regression.

Journal of Econometric Methods, 5(1):79–101, 2016.

• [8]

James Bergstra and Yoshua Bengio.

Random search for hyper-parameter optimization.

Journal of Machine Learning Research, 13(Feb):281–305, 2012.

• [9]

James S Bergstra, Rémi Bardenet, Yoshua Bengio, and Balázs Kégl.

Algorithms for hyper-parameter optimization.

In Advances in Neural Information Processing Systems, pages 2546–2554, 2011.

• [10]

Jérôme Bolte, Shoham Sabach, and Marc Teboulle.

Proximal alternating linearized minimization for nonconvex and nonsmooth problems.

Mathematical Programming, 146(1-2):459–494, 2014.

• [11]

Stephen Boyd and Lieven Vandenberghe.

Convex optimization.

Cambridge university press, 2004.

• [12]

Emmanuel J Candès, Xiaodong Li, Yi Ma, and John Wright.

Robust principal component analysis?

Journal of the ACM (JACM), 58(3):11, 2011.

• [13]

Venkat Chandrasekaran, Sujay Sanghavi, Pablo A Parrilo, and Alan S Willsky.

Sparse and low-rank matrix decompositions.

IFAC Proceedings Volumes, 42(10):1493–1498, 2009.

• [14]

David L Donoho.

Compressed sensing.

IEEE Transactions on information theory, 52(4):1289–1306, 2006.

• [15]

Derek Driggs, Stephen Becker, and Aleksandr Aravkin.

Adapting regularized low-rank models for parallel architectures.

SIAM Journal on Scientific Computing, 41(1):A163–A189, 2019.

• [16]

Peter J Huber.

Robust statistics, volume 523.

John Wiley & Sons, 2004.

• [17]

Frank Hutter, Holger H Hoos, and Kevin Leyton-Brown.

Sequential model-based optimization for general algorithm configuration.

In International Conference on Learning and Intelligent Optimization, pages 507–523. Springer, 2011.

• [18]

Aaron Klein, Stefan Falkner, Simon Bartels, Philipp Hennig, and Frank Hutter.

Fast bayesian optimization of machine learning hyperparameters on large datasets.

arXiv preprint arXiv:1605.07079, 2016.

• [19]

Roger Koenker and Kevin F Hallock.

Quantile regression.

Journal of economic perspectives, 15(4):143–156, 2001.

• [20]

M Kojima, N Megiddo, T Noma, and A Yoshise.

A Unified Approach to Interior Point Algorithms for Linear Complementarity Problems, volume 538 of Lecture Notes in Computer Science.

Springer Verlag, Berlin, Germany, 1991.

• [21]

Lisha Li, Kevin Jamieson, Giulia DeSalvo, Afshin Rostamizadeh, and Ameet Talwalkar.

Hyperband: A novel bandit-based approach to hyperparameter optimization.

arXiv preprint arXiv:1603.06560, 2016.

• [22]

A Nemirovskii and Y Nesterov.

Interior-Point Polynomial Algorithms in Convex Programming, volume 13 of Studies in Applied Mathematics.

SIAM, Philadelphia, PA, USA, 1994.

• [23]

Dominique Orban and Mario Arioli.

Iterative Solution of Symmetric Quasi-Definite Linear Systems.

SIAM, 2017.

• [24]

Yigang Peng, Arvind Ganesh, John Wright, Wenli Xu, and Yi Ma.

Rasl: Robust alignment by sparse and low-rank decomposition for linearly correlated images.

IEEE Transactions on Pattern Analysis and Machine Intelligence, 34(11):2233–2246, 2012.

• [25]

R Tyrrell Rockafellar and Roger J-B Wets.

Variational analysis, volume 317.

Springer Science & Business Media, 2009.

• [26]

Jasper Snoek, Hugo Larochelle, and Ryan P Adams.

Practical bayesian optimization of machine learning algorithms.

In Advances in neural information processing systems, pages 2951–2959, 2012.

• [27]

Robert Tibshirani.

Regression shrinkage and selection via the lasso.

Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288, 1996.

• [28]

Shiyi Tu, Min Wang, and Xiaoqian Sun.

Bayesian variable selection and estimation in maximum entropy quantile regression.

Journal of Applied Statistics, 44(2):253–269, 2017.

• [29]

Matthew A Turk and Alex P Pentland.

Face recognition using eigenfaces.

In Computer Vision and Pattern Recognition, 1991. Proceedings CVPR’91., IEEE Computer Society Conference on, pages 586–591. IEEE, 1991.

• [30]

S J Wright.

Primal-Dual Interior-Point Methods.

Siam, Englewood Cliffs, N.J., USA, 1997.

• [31]

Keming Yu and Rana A Moyeed.

Bayesian quantile regression.

Statistics & Probability Letters, 54(4):437–447, 2001.

• [32]

Zhengdong Zhang, Arvind Ganesh, Xiao Liang, and Yi Ma.

Tilt: Transform invariant low-rank textures.

International journal of computer vision, 99(1):1–24, 2012.