Fundamentals of Cone Regression: Mariella Dimiccoli

Statistics Surveys
Vol. 10 (2016) 53–99

ISSN: 1935-7516
DOI: 10.1214/16-SS114
Fundamentals of cone regression

Mariella Dimiccoli
Edifici O, Campus UAB
08193, Bellaterra (Cerdanyola del Vallès), Spain
e-mail: dimiccolimariella@gmail.com
Abstract: Cone regression is a particular case of quadratic programming

that minimizes a weighted sum of squared residuals under a set of linear in-
equality constraints. Several important statistical problems such as isotonic,
concave regression or ANOVA under partial orderings, just to name a few,
can be considered as particular instances of the cone regression problem.
Given its relevance in Statistics, this paper aims to address the fundamen-
tals of cone regression from a theoretical and practical point of view. Several
formulations of the cone regression problem are considered and, focusing on
the particular case of concave regression as an example, several algorithms
are analyzed and compared both qualitatively and quantitatively through
numerical simulations. Several improvements to enhance numerical stability
and bound the computational cost are proposed. For each analyzed algo-
rithm, the pseudo-code and its corresponding code in Matlab are provided.
The results from this study demonstrate that the choice of the optimiza-
tion approach strongly impacts the numerical performances. It is also shown
that methods are not currently available to solve efficiently cone regression
problems with large dimension (more than many thousands of points). We
suggest further research to fill this gap by exploiting and adapting classical
multi-scale strategy to compute an approximate solution.
MSC 2010 subject classifications: Primary 62; secondary 90.

Keywords and phrases: cone regression, linear complementarity prob-
lems, proximal operators.
Received March 2015.
Contents
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
2 Statement of the problem, basic notations and basic facts . . . . . . . 55
2.1 Convex quadratic programming (CQP) formulation . . . . . . . . 57
2.1.1 Primal formulation . . . . . . . . . . . . . . . . . . . . . . 57
2.1.2 Dual formulation . . . . . . . . . . . . . . . . . . . . . . . 57
2.2 Linear complementarity problem (LCP) formulation . . . . . . . 59
2.3 Proximal formulation . . . . . . . . . . . . . . . . . . . . . . . . . 59
3 State of the art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.1 Algorithms with asymptotic convergence . . . . . . . . . . . . . . 61
3.1.1 Least squares in a product space (LSPS) . . . . . . . . . . 61
3.1.2 Uzawa method . . . . . . . . . . . . . . . . . . . . . . . . 62
3.1.3 Hildreth’s algorithm . . . . . . . . . . . . . . . . . . . . . 62
3.1.4 Primal-dual interior point methods . . . . . . . . . . . . . 62
3.1.5 Dykstra’s algorithm . . . . . . . . . . . . . . . . . . . . . 63
53
54 M. Dimiccoli
3.1.6 Alternating Direction Method of Multipliers (ADMM) . . 64

3.2Algorithms with finite-time convergence . . . . . . . . . . . . . . 65
3.2.1 Properties of polyhedral convex cones with m ≤ n . . . . 65
3.2.2 Early algorithms based on the properties of polyhedral
convex cones . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.2.3 Mixed primal-dual basis algorithm . . . . . . . . . . . . . 68
3.2.4 Critical index algorithm: Nearest point problem in sim-
plicial cones . . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.2.5 Meyer’s algorithm . . . . . . . . . . . . . . . . . . . . . . 72
4 Issues about effectivness for large-scale problems . . . . . . . . . . . . 73
4.1 Suitability to take advantage of available good estimates . . . . . 73
4.1.1 PAV’s inspired approximate solution . . . . . . . . . . . . 75
4.2 Suitability to be implemented in an online fashion . . . . . . . . 75
4.3 Computational issues . . . . . . . . . . . . . . . . . . . . . . . . . 75
5 Improving the active set framework for concave regression problems . 77
6 Comparative results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
A Appendix A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
A.1 Signals used for comparative evaluations . . . . . . . . . . . . . . 80
A.2 Comparative evaluations on signals of size 50 . . . . . . . . . . . 82
A.3 Comparative evaluations on signals of size 500 . . . . . . . . . . . 85
A.4 Comparative evaluations on signals of size 1000 . . . . . . . . . . 88
B Pseudocode of the algorithms reviewed . . . . . . . . . . . . . . . . . . 91
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
Supplementary Material . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
1. Introduction
Cone regression analysis is a valuable alternative to more traditional parametric-

regression models, in all cases where the functional relationships between the
response (dependent) and the explanatory (independent) variables is unknown
and nonlinear and the constraints are a set of linear inequalities. Several impor-
tant statistical problems including isotonic, concave and constrained spline re-
gression, or ANOVA under partial orderings can be seen as particular instances
of the more general cone regression problem. Cone regression admits several
formulation approaches and implementation strategies, whose choice severely
impacts numerical performances. However, due to the little exposure to the
topics of optimization theory in modern-day Statistics, many optimization and
numerical approaches are commonly ignored by statisticians. This paper is a
contribution to fill this gap in the literature, addressing the fundamentals of
cone regression from a theoretical and practical point of view. With the goal of
going in deep with comparisons and numerical issues, we focus in particular on
the concave regression problem. In spite of its theoretical simplicity, since the
number of constraints increases linearly with the data size, concave regression
Fundamentals of cone regression 55
offers a good basis to discuss the fundamentals of cone regression and related
numerical issues.
The problem of concave regression is to estimate a regression function subject
to concavity constraints represented by a set of linear inequalities. Brought to
the attention of the scientific community by micro-economists interested in es-
timating production functions [1, 2, 3], the problem of concave regression arises
not only in the field of micro-economy (indirect utility, production or cost func-
tions, Laffer curve) but also in medicine (dose response experiments) and biology
(growth curves, hazard and failure rate in survival analysis). First addressed by
Hildreth in 1954 [1], the search for efficient methods for solving large concave
regression problems is still an open issue. This may appear quite surprising con-
sidering the noticeable advances in convex optimization since then, but it can
be understood when considering that most efforts have been devoted to theo-
retical issues such as generalizations and convergence while comparatively little
attention has been paid to the issues of efficiency and numerical performances
in practice [4, 5, 6].
In this paper, we formulate the cone regression problem by different opti-
mization approaches. We highlight similarities and difference between the var-
ious algorithms passed in review, we propose several improvements to enhance
stability and to bound the computational cost. We estimate the expected per-
formance of available algorithms, establishing in particular which is the most
competitive technique for solving large instances of the problem. Finally, in the
light of this study, we give recommendations for further research.
In section 2, we state formally the problem of cone regression also introducing
some basic notations and results that will be used throughout the paper. In
section 3, we survey the state of the art distinguishing between the class of
algorithms with asymptotic convergence and the class of algorithms with finite-
time convergence. In section 6, we make a numerical comparison of performances
and finally, in section 7, we draw some concluding remarks.
2. Statement of the problem, basic notations and basic facts
The aim of a regression analysis is to produce a reasonable analysis to the

unknown response function f , which can be modeled as
y = f (z) + (2.1)
where z ∈ R is the explanatory (dependent) variable, y ∈ Rd is the response

(independent) random variable, and is an error term, which is usually assumed
to be a mean zero random variable. Typically, one has observations on y and z
for n selected values of z. For each level of input, say zi , there may be several
trials and corresponding observations of output yi . Let Ti be the number of
trials at level of input zi and let yit be the observed output for the t−trial at
this level. Than we have
yit = f (zi ) + it , i = 1, ..., n t = 1, ..., Ti (2.2)

56 M. Dimiccoli
Inference about the response function may be drawn by assuming that the
function f (z) can be approximated by some given algebraic form with several
unknown parameters to be estimated from the data. However, the difficulty with
this procedure is that the inferences often depend critically upon the algebraic
form chosen. Alternatively, one may know some properties of the relation being
studied but does not have sufficient information to put the relation into any sim-
ple parametric form. In this case, a nonparametric approach is more appropiate.
Let xi be the expected value of output at input level zi :
xi = f (zi ) i = 1, ..., n (2.3)
Estimates of xi can be derived by the method of maximum likelihood, or by
the method of least squares or other formulations. If there were no a priori
restriction on f , than the maximum likelihood estimation of xi would just be
the mean of observed output for the level of input zi , that is
1
i=1
x̃i = yit i = 1, ..., n (2.4)
Ti
Ti
Instead, since in the cone regression problem the known property of the re-
gression function f can be expressed by a set of linear inequalities, to obtain
the maximum likelihood estimates, the likelihood function should be maximized
subject to the linear inequality constraints.
Formally, given a dataset of n dependent variables represented by the vectors
w, y ∈ Rn , corresponding to the independent variable values z1 < z2 < ... < zn ,
the problem of cone regression is to estimate the closest function to the dataset
via a least squares regression subject to a set of linear inequality constraints by
solving
x̂ = argmin x − y22,w (2.5)
{x ≤0}

n
with x − y22,w = wi (yi − xi )2
i=1
Denoting by Ki those vectors that satisfy the linear inequality constraints for
a fixed i, then Ki = ∅ is a closed convex set in Rn and the feasibility set K can
be written as the nonempty intersection of a family of closed subsets Ki ⊂ Rn .
Being the intersection of closed convex sets, the set K is also a closed convex
set. More precisely, since each Ki is an half-space which contains the origin, the
feasibility set K is a convex polyhedral cone. In matrix form, K can be written as
K = {x : Ax ≤ 0}. In the case of concave regression A ∈ Rm×n with m = n − 2
is a matrix such that each row Ai represents a concave two-piece linear function
with a negative second difference at xi+1 only and the linear inequalities are as
follows
xi+2 − xi+1 xi+1 − xi
− ≤ 0, i = 1, ..., n − 2 (2.6)
zi+2 − zi+1 zi+1 − zi
In the following, we give alternative formulations of the cone regression prob-
lem that rest on optimization theory.
2.1. Convex quadratic programming (CQP) formulation
2.1.1. Primal formulation
The problem (2.5) is to find the point x̂ in the cone K that is closest to y. The
solution is found at the orthogonal projection of y onto K, written as Π(y|K)
using the metric · 2,w , represented by the symmetric positive definite matrix
W.
x̂ = Π(y|K) = argmin(y − x)T W (y − x) (2.7)
{Ax≤0}
For the problem (2.7), the matrix W is diagonal with element wi on the diagonal.
In practice, if yi is the mean value measured at zi , than wi corresponds to the
size of the sample at zi . Since K is a closed, convex and non empy set on the
Hilbert space Rn , it is a set of Chebyshev, that is the projection exists and it is
unique.
2.1.2. Dual formulation
The dual formulation of problem (2.7) rests on the Moreau decomposition theo-
rem [7], which is a generalization in convex analysis of the orthogonal projection
theorem for vectorial sub-spaces. Central to the Moreau decomposition theorem
is the definition of polar cone of a given convex cone K, which is given below.
Definition 1. The polar cone Ko to any convex cone K is given by
Ko = {x ∈ Rn : ∀k ∈ K, k , x ≤ 0} (2.8)
The Moreau decomposition theorem is as follows.

Theorem 2.1. Let K ⊆ Rn be a closed convex cone, Ko its polar cone and
y ∈ Rn a given point. Then the following assertions are equivalent:
(i) x̂ = argmin ||x − y||2 , x̂o = argmin ||x − y||2

x∈K x∈Ko
(ii) x̂ ∈ K, x̂ ∈ K , x̂, x̂ = 0, y = x̂ + x̂o
o o o
By relying on this theorem we can alternatively solve problem (2.7) by first

finding the projection on the polar cone x̂o and then computing the solution to
the primal problem as the difference x̂ = y − x̂o (see Figure 1). This alternative
is attractive since, as it will be clarified below, it implies an analytically simpler
form for the constraints.
Before stating an important Lemma about the relationship between the polar
cone and the constraint matrix A, let us introduce the definition of edges of a
polyhedral convex cone.
Definition 2. Let K be a polyhedral convex cone in Rn , then the vectors
ei∈ Rn \ {0} are the edges or generators of K if and only if K = pos({ei }) =
{ ki ei |k ≥ 0}.
58 M. Dimiccoli
Fig 1. The polar cone Ko of a given convex cone K ⊂ R2 is given by the set of all vector
whose scalar product with vectors of K is negative. The data point y can be written as the
sum of x̂, the projection onto the cone K and x̂o , the projection onto the polar cone Ko .
Intuitively speaking, the edges of a polyhedral convex cone are one-dimen-

sional rays, which always pass through a fixed point (the vertex).
Lemma 1. The mrows of the matrix A are the edges of the polar cone, that is
Ko = {x : x = i=1 ATi ai , ai ≥ 0}.
m
To see that, observe thatKo = { i=1 ai ATi , ai ≥ 0} is polar to K since
m
∀x ∈ K,∀ρ ∈ Ko : ρ, x = { i=1 ai ATi , x ≤ 0}, which is the definition of a
polar cone of K. Conversely, K is polar to Ko since: K = (Ko )o = K.
By relying on this result Khun-Tucker [8] proved the following theorem:
Theorem 2.2. The primal constrained quadratic minimization problem (2.7)
is equivalent to the dual problem
λ̂ = argmin(y − AT λ)T W (y − AT λ) (2.9)

λ≥0
Denoting by λ̂ the solution to the dual problem, the solution to the primal prob-
lem is x̂ = y − AT λ̂.
As it can be observed, in the dual formulation each element of the vector λ
must satisfy a single positivity constraint.
Goldman [9] noticed that dual problem can be also viewed as a minimum
distance problem in the same parameter space as the primal problem.
x̂ = argmin ||x||2 (2.10)

x∈C
where C = {x|x = y − AT λ, λ ≥ 0} is a rotation of the dual cone with its vertex

translated to y. x̂ also solves the re-parametrized dual problem.
2.2. Linear complementarity problem (LCP) formulation
The CQP (2.7) can be recast as a linear complementarity problem (LCP). To

see that, let us consider the Lagrangian associated with problem (2.7).
L(x, λ) = ||x − y||22,w + < λ, Ax > (2.11)
where λ ≥ 0 is the vector of dual variables associated to each of the convexity
constraints. By applying the Karush–Kuhn–Tucker (KKT) optimality condi-
tions [8] to (2.11), that is
∇L(x̂, λ̂) = 0 (2.12)
λ≥0 (2.13)
λT Ax̂ = 0 (2.14)
we obtain the equivalent LCP
w + Mλ = q (2.15)
w ≥ 0, λ ≥ 0, wT λ = 0. (2.16)
where w = −AT x, M = −AAT and q = −AT y. Note that by dropping the con-
stant term from the Lagrangian and dividing it by 2: L(x, λ) = 12 xxT − y T x +
λT Ax. Therefore: ∇L(x, λ) = xT − y T + λT A = 0. By taking the transpose and
multiplying by −A: −Ax + (−AAT )λ = −AT y. This LCP has a unique com-
plementary solution. Denoting by (ŵ, λ̂) its solution, λ̂ is the optimal solution
of the dual problem (2.9).
The condition wT λ = 0 is called the complementarity condition and the way
in which it is dealt with determines if the optimization algorithm belongs to the
class of interior point methods that will be introduced in section 3.1.4 or to the
class of active set methods that will be detailed in section 3.2.
2.3. Proximal formulation
The CQP (2.7) can be solved by using a proximity operator [10, 11]. Proximity
operators are used to solve problems of the form
argmin f1 (x) + f2 (x)... + fm (x) (2.17)
x∈Rn
where f1 , f2 , ..., fm are convex functions from Rn to ] − ∞, +∞], that are not
necessarily differentiable. Each fi is treated through its proximity operator,
which is defined as follows.
Definition 3. Let Γ0 (Rn ) be the class of lower semicontinuous convex functions
from Rn to ] − ∞, +∞] such that their domain, denoted by dom(f ) is not the
empty set. Let f be a function f ∈ Γ0 (Rn ), then the proximity operator of f is
proxf (x) : Rn → Rn such that
1
∀x ∈ Rn , proxf (y) = argmin f (x) + ||x − y||2 (2.18)
x∈Rn 2
60 M. Dimiccoli
The proximity operator is characterized by the property
∀(x, p) ∈ Rn × Rn , p = proxf (y) ⇐⇒ y − p ∈ ∂f (p), (2.19)
where ∂f : Rn → 2R is the subdifferential of f :

n

∂f = u ∈ Rn : ∀y ∈ Rn : (y − p)T u + f (p) ≤ f (y)) (2.20)
The proximity operator of a convex function is a generalization of the projection

operator onto a closed convex set C. To see that let us consider the indicator
function of C
0 if x ∈ C
ıC (x) =
+ ∞ if x ∈ C
By using the fact that minimizing J(x) over C is equivalent to minimizing J(x)+
ıC (x) over Rn
argmin J(x) = argmin{J(x) + ıC (x)}
x∈C x∈Rn
it results that proxıC = Π(·|C).

The solution to problem (2.7) can be therefore understood as the proximity
operator over K
x̂ = argmin ||x − y||2 = argmin{||x − y||2 + ıK (x)} (2.21)

x∈K x∈Rn
= proxıK y (2.22)
Alternatively, using the fact that (K1 ∩ ... ∩ Km )o = K1o + ... + Km

o
, the dual
problem (2.7) can be seen as the proximity operator over the sum of m indicator
functions of the convex sets Kio :

m
x̂o = argmin

||x − y||2 = argmin{||x − y||2 + ıKio (x)} (2.23)
x∈ m
i=1 Kio x∈Rn i=1
= prox m
i=1 ıKo
y (2.24)
i
Intuitively, the base operation of a proximal algorithm is evaluating the prox-

imal operator of a function, which itself involves solving a small convex opti-
mization problem. These subproblems often admit closed form solutions or can
be solved very quickly with standard or simple specialized methods.
3. State of the art
In this section, we review existing algorithms for solving the cone regression
problem. The algorithmic approaches, and in turn their numerical performances,
strongly depend on the choice of the problem formulation. All existing methods
are iterative and they can attain the optimal solution since K is a Chebyshev
set and therefore the optimal solution must exist in this closed set. However,
in terms of their numerical perfomances they can be classified into two broad
classes: the class of methods that never (or in very simple cases) attain the opti-
mal solution [1, 12] and those of methods that converge to the optimal solution
in a finite number of steps [13, 14, 15, 16, 17, 18, 19, 20]. As it will be clarified
in the following, methods with asymptotic convergence rest on the properties of
the sub-gradient or more in general of proximity operators and act by finding
the solution as the limit of a sequence of successive approximations. They are
typically derived from the primal, the dual or from the proximal formulation.
Methods with finite-time convergence exploit the geometric properties of poly-
hedral convex cones and find the exact solution as a non-negative linear combi-
nation of functions, forming a basis in a specified finite dimensional space. They
are typically derived from the linear complementarity problem formulation.
3.1. Algorithms with asymptotic convergence
This section includes algorithms based on the primal formulation such as Least
Squares in a Product Space (section 3.1.1), algorithms based on the dual for-
mulation such as the Uzawa’s method (section 3.1.2) and Hildreth’s method
(section 3.1.3), algorithms that solve the dual problem simultaneously with the
primal problem such as the Dykstra’s alternating projection method (section
3.1.5) and algorithms based on the proximal formulation such as Alternating
Direction Method of Multipliers (section 3.1.6).
3.1.1. Least squares in a product space (LSPS)
Since the Euclidean space Rn equipped with the dot product is an Hilbert space
H, the problem (2.7) can be recast in the m−fold Cartesian product of H, say
Hm [35]. Let Km be the Cartesian product of the sets (Ki )i∈I , i.e., the closed
convex set Km = ×i∈I Ki = {x ∈ H : ∀i ∈ I : xi ∈ Ki } and let D be the diagonal
vector subspace, i.e. D = {(x, ..., x) ∈ Hm : x ∈ H}. Then, the CQP (2.7) is
equivalent to
argmin x − ȳ22,w (3.1)
{x∈Km ∩D}
where ȳ = (y, ..., y). Using this strategy, the problem of projecting onto the
intersection of m convex sets is reduced to the problem of projecting, in a
higher dimensional space, onto only two convex sets, one of which is a simple
vector subspace. Geometrically, this is equivalent to find a point in D which is
at minimum distance from Km . This point can be obtained iteratively by
xk+1 = xk + λk (PD ◦ PKm (xk ) − xk ) (3.2)
The advantage of this strategy is that it speeds up the convergence since a bigger
(than 2, which is the upper bound for Fejér sequences [36]) relaxation interval
can be allowed.
62 M. Dimiccoli
3.1.2. Uzawa method
A classical method to solve a convex minimization problem subject to inequal-

ity constraints is the Uzawa method [37], which search directly for the sad-
dle point of the Lagrangian (2.11). In fact, if the Lagrangian L(x, λ) admits a
saddle point, say (x̂, λ̂), then the duality gap δ = minx∈K maxλ∈R+ L(x, λ) −
maxλ∈R+ minx∈K L(x, λ) is null and x̂ is a critical point of the Lagrangian. Since
the dual function H(λ) = argminx∈K L(x, λ) is differentiable, it can be mini-
mized explicitly by using the gradient descent method. Therefore the Uzawa
method alternates a minimization step over Rn with respect to x with λ fixed
and a maximization step with respect to λ onto R+ , with x fixed. The algo-
rithmic parameter ρ > 0 can be fixed to optimize convergence by relying on
theoretical considerations. Therefore the CQP (2.9) is equivalent to finding
x̂ = argmin argmax L(x, μ) = ||x − y||2 + < μ, Ax > (3.3)
x μ≥0
3.1.3. Hildreth’s algorithm
Hildreth [1] proposed to apply the Gauss Seidel algorithm [21] to the dual prob-
lem (2.9). A single cycle of the Hildreth’s algorithm consists in updating each ele-
ment of λ sequentially in an arbitrary fixed order. Therefore each cycle consists of
m steps, each of which corresponds to a projection onto the cone Ki ,i = 1, ..., m.
The algorithm gives rise to a sequence of points, each one differing from the pre-
ceding in exactly one coordinate. At the cycle k + 1, the λk+1 i is used in the
estimation of the point λk+1 i+1 so that the best available estimations are used
for computing each variable. The convergence of the Gauss Seidel algorithm is
guaranteed only if the matrix A is full row rank, so that there are not redun-
dancies among the inequality restrictions, and it is guaranteed independently of
the initial point λ0 only if A is positive definite and symmetric. The algorithm
is sensitive to the normalization as well as to the order of the projections.
3.1.4. Primal-dual interior point methods
First introduced by Karmakar in 1984 [43], primal-dual interior point methods

act by perturbing the complementarity condition wT λ = 0 in the LCP formula-
tion (2.15) and replacing with wT λ = μ. The partition of vectors w and λ into
zero and nonzero elements is gradually revealed as the algorithm progresses by
forcing a reduction of μ. All iterates satisfy the inequality constraints strictly.
The solution is approached from either the interior or exterior of the feasible
region but never lie on the boundary of this region. Let the function F (x, λ, w)
be such that the roots of this function are solutions to the first and the last
optimality conditions in (2.15).
⎛ ⎞
w − AT x
⎜ T t ⎟
Fμ (x, λ, w) = ⎝x − y + λ A⎠
T
wT λ − μe
The perturbed complementarity condition introduces a nonlinearity, there-

fore for each fixed μ > 0 a system of nonlinear equations should be solved.
The nonlinear system is typically solved by using a Newton-like algorithm [50].
Each iteration of the Newton’s method finds a search direction from the current
iterate (xk , λk , wk ) . It is computationally expensive, but can make significant
progress towards the solution. For instance in barrier methods, which are the
most efficent of the family, this is achieved by using a penalizing term, called a
barrier function, for violations of constraints whose value on a point increases
to infinity as the point approaches the boundary of the feasible region. Interior
point methods must be initialized at an interior point, or else the barrier func-
tion is undefined. The interested reader is referred to [44] for further information
about interior point methods.
3.1.5. Dykstra’s algorithm
In 1983, Dykstra [12] proposed a generalization of the Hildreth’s procedure

applicable to the case of constraints corresponding to more general convex cones
than polyhedral convex ones. The Dykstra’s algorithm is based on the idea,
before suggested by Von Neumann [22] to the case of subspaces, of computing
the projection onto the intersection of convex sets by relying on the solution of
the simpler problem of projecting onto the individual sets. In the case of concave
regression the projection onto a single convex set Ki involves only three points
and, if the constraint is not satisfied, it corresponds to the straight line fitting
the points yi ,yi+1 ,yi+2 .
Dykstra’s algorithm iterates by passing sequentially over the individual sets
and projects onto each one a deflected version of the previous iterate. More
precisely, before projecting onto the cone Ki during the (k + 1)−th cycle, the
residuum obtained when projecting onto Ki at the previous k−th cycle, say Rik
is removed and a new residuum associated to the cone Ki , say Rik+1 is computed
after the projection. In practice, each xki is the projection of y +R1k +...+Ri−1
k
+
Ri+1 + ... + Rm onto Ki , where Ri = xi − (y + R1 + Ri−1 + Ri+1 + ... + Rm ).
k k k k k k k−1 k−1
If each Ki is a subspace, then at each new cycle k +1 the residuum −Rik of the
projection onto each convex cone Ki is of course the projection of xk onto the
cone Kio . Therefore, the Dykstra’s procedure for subspaces reduces to exactly
the cyclic, iterated projections of von Neumann. In this case, for k → ∞, the
sum of the residua over the cones Kik approximates the projection of y onto
the polar cone Ko and therefore, for the Moreau decomposition theorem, xk
approximates the projection onto K.
However, if the Ki are not subspaces, Π(·|Ki ) is not a linear operator and
then the von Neumann algorithm does not necessarily converge.
The Dykstra’s algorithm can also be interpreted as a variation of the Douglas–
Rachford splitting method applied to the dual proximal formulation (2.23).
The seminal works of Hildreth and Dykstra have inspired many studies
mostly devoted to theoretical investigations about their behavior in a Hilbert
space [23, 24], about its convergence [25, 26], about its relation to other meth-
64 M. Dimiccoli
ods [27, 28] and about its interpretation in more general frameworks, such as
the proximal splitting methods [29]. Han [30], as well as Iusem and De Pierro
[25], showed that in the polyhedral case, the method of Dysktra becomes the
Hildreth’s algorithm and therefore it has the same geometric interpretation of
Gauss Seidel to the dual problem. Gaffke and Mathar (1989) [27] showed the
relation of the Dysktra algorithm to the method of component-wise cyclic min-
imization over a product space, also proposing a fully simultaneous Dykstra
algorithm. The only works devoted to give some insight about a more efficient
implementation are the ones of Ruud. Goldman and Ruud [9](1993) generalized
the method of Hildreth showing that there is no need to restrict the iterations
to one element of λ at a time. One can optimize over subsets or/and change
the order in which the elements are taken. This observation is important for
the speed of convergence since the slow speed can be understood as a symptom
of near multicollinearity among restrictions. Because the intermediate projec-
tions are so close to one another, the algorithm makes small incremental steps
towards the solution. They also remarked that Dykstra uses a parametrization
in the primal parameter space, that causes numerical round off errors in the
variable residuum. These round off errors cumulate so that the fitted value does
not satisfy the constraints of the dual problem. It would be better to use a
parametrization on the dual so that the contraints in the dual would be sat-
isfied at each iteration. Later, Ruud [31] proved that the contraction property
of the proposed generalizations rests solely on the requirement that every con-
straint appears in at least one subproblem of an iteration. As one approaches
the solution, constraints that are satisfied at the solution are eliminated. To re-
move satisfied constraints would accelerate the Hildreth procedure. The authors
propose to reduce the set of active constraints, that is the constraints satisfied
as an equation at the corresponding points, by removing as many constraints as
possible through periodic optimization over all positive elements of λ.
3.1.6. Alternating Direction Method of Multipliers (ADMM)
ADMM is an augmented Lagrangian technique [32, 33] which can be applied to

problems of the form
F ind argmin ||y − x||2 + g(Ax) (3.4)

z∈Rm ,Ax=z,z≤0
where the matrix A is assumed to be irreducible (AAT = vI, v > 0) and the
intersection of the relative interiors of the domains of the two functions is as-
sumed to be not empty (ri dom(g) ∩ ri dom(f ) = ∅). ADMM minimizes the
augmented Lagrangian L over the two variables of the problems, say x and z,
first x with z fixed, then over z with x fixed, and then applying a proximal
maximization step with respect to the Lagrange multiplier λ. The augmented
Lagrangian of index γ ∈ [0, ∞] is
1 T 1
L(x, z, y) = f (x) + g(z) + λ (Ax − z) + ||Ax − z||2 (3.5)
γ 2γ
where f (x) = ||y − x||2 . Denoting by proxA

f the proximal operator which maps
a point z ∈ R to the unique minimizer of f (x) + ||Ax − z||2 and denoting by
n
proxg = proxf ◦A the implementation detailed in the Appendix B is obtained.

The ADMM method rests on the proximal formulation 2.21. Indeed, it can
be viewed as an application of the Douglas-Rachford splitting algorithm [34].
3.2. Algorithms with finite-time convergence
All algorithms that will be reviewed in this section are active set methods resting
on the LCP formulation 2.15. Active set methods work by choosing a subset of
indices j ∈ J˜ ⊂ J = {1, ..., n} such that wj is allowed to be non-zero and forcing
the corresponding λj to be zero, while the remaining indices j ∈ J \ J˜ force wj to
be zero and allow λj to take nonzero values. In this section we will review active
set algorithm such as the mixed primal-dual basis algorithm (section 3.2.3),
the critical index algorithm (section 3.2.4) and the Meyer’s algorithm (section
3.2.5).
Before detailing the algorithms, we introduce some definitions and basic re-
sults about the geometry of polyhedral convex cones, on which are based the
algorithms presented in this section. For further details the reader is referred to
[38].
3.2.1. Properties of polyhedral convex cones with m ≤ n
Lemma 1 establishes the relationship between the constraint matrix A and the
edges of the polar cone {γ i , i = 1, ..., m}, namely AT = [γ 1 , ..., γ m ]. We would
like now to determine the edges of the constraint cone K.
Let the vectors {γ m+1 , .., γ n } be orthogonal to the vectors {γ i , i = 1, .., m}
and orthonormal to each other so that the set {γ i , i = 1, .., n} forms a basis for
Rn . By defining the dual basis of the basis {γ i , i = 1, .., n} as the set of vectors
{β i , i = 1, ..., n} that verify the relationship

−1 i = j
(β i )T γ j = (3.6)
0 i = j
the constraint cone K = {x : Ax ≤ 0} can be equivalently written as
m
n
K= x:x= bi β i + bi β i , bi ≥ 0, i = 1, ..., m (3.7)
i=1 i=m+1
To see that let B = [β 1 , ..., β n ] and C = [γ 1 , ..., γ n ]. Then Ax are the first
m coordinates of Cx. Since B T C = −In by construction, by multiplying both
members at left for B −1 and at right for x, we obtain: Cx = −B −1 x. Therefore
Cx gives the negative coordinates of x in the basis {β i , i = 1, ..., n}. Further-
n in K have their first m coordinates non-negative and can be written
more, points
as x = i=1 bi β i , where bi ≥ 0 for i = 1, .., m.
Taking into account Def. 2, (3.7) established that the vectors β i are the edges
of the constraint cone K.
66 M. Dimiccoli
Definition 4. Let K be a polyhedral convex cone in Rn , then F ⊆ K is a face

of K if and only if F is the intersection of K with a supporting hyperplane.
A polyhedral convex cone arises as the intersection of a finite number of half-
spaces whose defining hyperplanes pass through the origin. The i−th row of A
is normal to the hyperplane generating the i−th closed half-space.
The following Lemma, proved by Rockafellar in 1970 [39], establishes a re-
lationship between the first m vectors of the dual basis {β i , i = 1, ..., m} and
the faces of the cone K = K ∩ span(K), where span(K) denotes the subspace
spanned by the m edges of K.
Lemma 2. Let K = {x : Ax ≤ 0}, where AT = [γ1 , ..., γm ], be the constraint
cone and let {β i , i = 1, ..., m} be the dual basis of {γ i , i = 1, ..., m}. Denoting
by span(K) the subspace spanned by the m edges of K, let K = K ∩ span(K) =
{x ∈ R : x = j=1 bj β j , bj ≥ 0}. Then, for J ⊆ {1, ..., m} the faces of K are
n m
the sets:

x ∈ Rn : x = bj β j , bj ≥ 0
j∈J
The set of all relatively open faces

FJ = x ∈ Rn : x = bj β j , bj > 0
j∈J
forms a partition of K .
Denoting by u and v the projections of y onto span(K) and span({γ m+1 , ...,
γ }) respectively, being the two subspaces orthogonal, y can be written as y =
n
u + v. Since v can be easily computed as v = X(X T X)−1 X T y with X =

[γ m+1 , ..., γ n ], the problem (2.7) reduces to find the projection of u onto K .
The next Lemma, proved by Zarantonello [40] focuses on the set of points in
span(K) projecting on a given point x ∈ K , say Π−1 K (x). Before stating it we
need to define the support cone of a closed convex set at a given point.
Definition 5. The support cone of a closed convex set K at x denoted by LK (x)
is the smallest convex cone with vertex at the origin containing K − x.
Lemma 3. Let K , FJ , {γi , i = 1, ..., n} and {βi , i = 1, ..., n} defined as in
Lemma 2. If x is a point of K belonging to the open face FJ , then:

• Π−1
K (x) = x + L o
K (x) = x + i∈J
/ c i γ i
, c i ≥ 0 = j
j∈J bj β +

/ ci γ , bi > 0, ci ≥ 0 ,
i
i∈J

where x = j∈J bj β j .
• The sets Π−1 K (x) are disjoint closed convex cones.
• ∪x∈K Π−1K (x) = span(K)
where LoK (x) denotes the dual of the support cone of K at x.
Fig 2. The point x1 ∈ K belongs to the open face FJ = {x ∈ K : x = bβ 1 , b > 0}, with
T
J = {1}. The support cone of K at x1 is LK (x1 ) = {x : γ 2 x ≤ 0} and its dual is LoK (x1 ) =
{x : x = cγ , c ≥ 0}. The set of points that project onto x1 is given by the set Π−1
2
K (x1 ) =
{x1 + LoK (x1 )} = {x1 + cγ 2 , c ≥ 0}. The point x3 ∈ K belongs to the open face FJ = {0},
with J = ∅. The support cone of K at x3 is K and its dual is Ko , so that the set of points
projecti onto x3 is {K }. The point x4 ∈ K belongs to the open face FJ = {x ∈ K : x =
that o
i=1,2 bi β , bi > 0}, with J = {1, 2}. The support cone of K at x4 is the origin and its dual
is the origin, so that the set of points that project onto x4 is {x4 }.
Then any point in span(K) projects onto a unique point in K and belong to
a unique non-negative orthant, or sector SJ

SJ = x ∈ R n : x = bj β j + cj γ j , bj > 0, cj ≥ 0 (3.8)
j∈J j ∈J
/
Figure 2 illustrates this result for span(K) = R2 . Points in SJ project onto

the subspace spanned by the vectors {β j , j ∈ J}, that is on the face FJ =

j∈J bj β . Vectors belonging to K project onto the origin, vectors belonging
j o
to K project on themself, while each other vector of Rn project onto a unique

face of K.
Therefore, if the sector SJ containing the vector u is known, then the projec-
tion of u onto K can be easily computed as the projection of u onto the subspace
spanned by the {β j , j ∈ J}. This reduces the problem of projecting y onto K
to the problem of finding the set of indices Jˆ such that the sector SJˆ contains
u. The set complementary of Jˆ with respect to {1, ..., m} corresponds to the
indices of the constraints satisfied at equality in the optimal solution.
The algorithms described in this section propose different strategies to find
the optimal set J.ˆ
3.2.2. Early algorithms based on the properties of polyhedral convex cones
The first algorithm addressing the problem of projecting a point y ∈ Rn onto a

polyhedral convex cone K ⊂ Rn by a non-asymptotic procedure dates back to
68 M. Dimiccoli
i
work of Wilhelmsen
[13] in 1976.
Wilhelmsen assumes
that the m generators β
m
of the cone K = x ∈ R : x = i=1 bi β , bi > 0 are known and propose an
n i
algorithm which computes a sequence of nearest points xk to y in subcones Kk

of K. Each subcone Kk is chosen so that xk ∈ int(K k ) and is closer to y than
is xk−1 . This means that xk is in the near side of the supporting hyperplane of
Kk−1 with respect to y. The key step is to find xk+1 given xk and the proposed
procedure to do that is laborious.
Pshenichny and Danilin (1978) [14] proposed an algorithm similar to the
one of Wilhelmsen which also converges in a finite number of steps. In both
algorithms, m can be any integer even larger than n. A more efficient procedure,
but with the more restrictive assumption that m ≤ n has been proposed by
Fraser and Massam in 1989.
3.2.3. Mixed primal-dual basis algorithm
Fraser and Massam [16] proposed an iterative algorithm to solve the general
problem of projecting a data point y ∈ Rn onto a polyhedral convex cone
K ⊂ Rn generated by m ≤ n linear inequality constraints.
Polyhedral convex cones generated by a number of inequalities at least equal
to the dimension of the space they belong to have been the subject of section
3.2.1. As seen there, the problem of projecting a data point onto this class of
cones can be reduced to find the set of edges, or generators of the cone, indexed
by Jˆ ⊆ 1, ..., n such that the sector SJˆ contains the data point.
To this goal, the set of edges of the polar cone {γ i , i = 1, ..., m} is completed
by n − m vectors orthogonal to {γ i , i = 1, ..., m} and orthonormal to each other.
In the case of concave regression m = n − 2, the set is completed by a constant
function γ m+1 = 1/||1||, where 1 is the m-dimensional unitary vector and by
a linear
m function γ
m+2
= (x − x̄1)/||(x − x̄1)||, where x = (x1 , ..., xm ) and
x̄ = i=1 xi /m. The set of vectors {γ i , i = 1, ..., n} form a basis for Rn .
Let the vectors {βi , i = 1, ..., n} be the dual basis of the basis {γ i , i = 1, ..., n}
as defined in (3.6). Fraser and Massam called the vectors β i and γ i primal and
dual vectors respectively. A primal-dual basis for Rn , associated to the set of
indices J ⊆ {1, ..., n} ≡ L is a basis BJ = [α1 , ..., αn ] made up of a subset of
the primal basis vectors {βi }i∈J and a complementary subset of the dual basis
vector {γi }i∈L\J . For n = m the primal basis vectors, corresponding to the
edges of K, are simply the columns of −(AT )−1 . Using the above definitions,
the problem of projecting a point y ∈ Rn onto the cone K can be formulated as
follows.
is equivalent to the problem of finding
argmin ||u − x||2 (3.9)
x∈K
where u = y − v, with v = Π(y|span(γ m+1 , ..., γ n )). Denoting by xu the solution

to this problem, the solution to the primal problem (2.7) is x̂ = xu + v.
Finding the sector containing u is achieved moving along a fixed line joining
an arbitrary chosen initial point x0 inside the cone or on its boundary to the
data point u. By moving along a fixed line, many sectors are crossed. Each time
a sector is crossed the successive approximation xk is obtained by projecting the
point on the line passing through u and x0 on the face FJ k of K (see Lemma 2)
so that the distance ||u − xk || is decreasing in k.
At each iteration each basis differs from another by one vector only and
therefore the coordinates of xk onto the new basis are easy to calculate. What
changes is that one coordinate becomes equal to zero and therefore there is
no need to estimate the inverse of a matrix at each step. For this reason this
algorithm is faster that the one of Wilhelmsen.
Points belonging to the sector SJ k have non-negative coordinates in the mixed
primal-dual basis BJ k relative to the cone K. Therefore the procedure terminates
when the coordinates of the data point in the primal dual basis BJ k are all non-
negative, meaning that the point xk is on the face of the sector containing the
data point u.
The number of iterations needed is equal to the number of different sectors
that the line joining the initial point to the data point has to cross. This number
is bounded above by 2n . It is worth to remark that crossing a sector corresponds
to a pivot step in the equivalent LCP. In fact, the net effect of a pivot step is
of moving from the point xk contained in the face FJ k of K to the point xk+1
contained in the face FJ k+1 of K.
Ten years later, Meyer [17] generalized the algorithm of Fraser and Massam
to the case of more constraints than dimensions, that is when m > n.
3.2.4. Critical index algorithm: Nearest point problem in simplicial cones
Murty and Fathy (1982) [19] considered the general problem of projecting a
given vector y ∈ Rn onto a simplicial cone K ⊂ Rn . The definition of simplicial
cone and pos cone, on which the former definition is based are as follows.
Definition 6. The pos cone generated by the
mvectors in Δ = {δ i , i = 1, ..., m},
denoted by pos(Δ), is the set {x ∈ R |x = i=1 di δ , di ≥ 0}.
n i
Definition 7. A cone K ⊂ Rn is said to be simplicial if it can be expressed as

a positive linear span of n linearly independent vectors Δ = {δ i , i = 1, ..., n} in
Rn (i.e., a basis for Rn ): K = pos(Δ).
For any point x ∈ pos(Δ) the vector d = D−1 x, where D = [δ 1 , ..., δ n ] is
called combination vector corresponding to x. Therefore the projection x̂ of y
onto pos(Δ) can be expressed as a non-negative linear combination of the edges
n
of the cone: x̂ = i=1 dî δ i , where the optimal combination vector corresponding
to x̂ is dˆ = D−1 x̂.
Murty and Fathy named the set of indices Jˆ ⊆ {1, ..., n} such that dî∈Jˆ > 0
set of critical indices. Using the definitions of simplicial cone and combination
vector, the original problem of projecting the point y ∈ Rn onto the cone K can
be formulated as follows.
70 M. Dimiccoli

is equivalent to the problem
dˆ = argmin(u − Dd)T W (u − Dd) (3.10)

d≥0
where u = Π(y|span(K)) and D = [γ 1 , ..., γ n ], with {γ i = ATi , i = 1, ..., m} and

{γ i , i = m + 1, ..., n} defined as in section 3.2.3.
Denoting by dˆ the solution to this problem, the solution to the primal problem
(2.7) is x̂ = y − Ddˆ + v.
This formulation has the same structure of the dual formulation (2.9), where
the combination vector d in (3.10) plays the same role as the dual variable λ in
(2.9). The only difference is that in (2.9) the matrix A ∈ Rm×n is not squared
as the matrix D ∈ Rn×n .
As shown in section 2.2 for the dual quadratic formulation (2.9), the formu-
lation (3.10) can be recast as a LCP. This equivalency is important since the
following theorem, proved in [19] applies to the LCP formulation of (3.10).
Theorem 3.3. If a single critical index for the LCP problem of order n is
known, the problem can be reduced to a LCP of order n − 1.
The fact that finding a critical index reduces the dimension of the prob-
lem can be argued geometrically. Let l be a critical index, then denoting by
N P P [Γ; u] the subproblem (3.10), where Γ = {γ i , i = 1, ..., n} its solution is
l
(uT γ l )
also the solution to the N P P [Γ ∪ {−γ l }; u]. Defining ū = u − γ ||γ l ||2 and Γ̄ =
l T l
{γ̄ 1 , ..., γ̄ l−1 , γ̄ l+1 , ..., γ̄ n }, where γ̄ i = γ i − γ ||γ
(u γ ) i
l ||2 γ , then γ̄ , i ∈ {1, ..., n} \ l is
i
orthogonal to γ and the cone pos(Γ ∪ {−γ }) is the direct sum of the full line
l l
generated by γ l and the simplicial cone pos(Γ̄).

Solving N P P [Γ̄, ū] is an n − 1 dimensional problem. If x̂∗ is the solution
of the N P P (Γ̄, ū), then the solution x̂ to the N P P [Γ; u] is obtained as x̂ =
l
(uT γ l )
x̂∗ + γ ||γ l ||2 .
By relying on Theorem 3.3, the authors proposed an algorithm consisting of a

subroutine to identify a critical index for the problem, followed by a subroutine
which reduces the size of the problem once a critical index is found. Since the
solution dˆ is the orthogonal projection onto the linear hull of {γ i , i ∈ J}, ˆ if Jˆ
is known, the solution of the equivalent LCP (q, M ) and correspondingly the
solution to N P P [Γ; u] can be easily found.
The routine to identify a critical index operates on the N P P [Γ; u] by exploit-
ing the geometric properties of projection faces of a pos cone, whose definition
is as follows.
Definition 8. Let S ⊂ Γ. pos(S) is a face of pos(Γ) of dimension |S|. pos(S) is
said to be a projection face of pos(Γ) if Π(u|span(S)) ∈ pos(S).
In the following, we list some theorems on which the critical index algorithm
is based. Their proofs can be found in [41] (chapter 7).
Theorem 3.4. Let S ⊂ Γ, S = 0. The optimum solution of N P P [S; u] is in

the relative interior of pos(Γ)if and only if the projection of u onto the linear
span of S is in the relative interior of pos(S): Π(u|span(S)) ∈ ri(pos(S)).
Theorem 3.5. Let x̂ = Ddˆ be the optimum solution of N P P [Γ; u]. Let Jˆ the
ˆ Then pos(S) is a projection face of
set of critical indices and S = {γ j , j ∈ J}.
pos(Γ).
Theorem 3.5 tells us that the projection onto the cone Ko belongs to the
pos cone generated by the set S of vectors corresponding to critical indices and
that such pos cone is a projection face. Therefore, if the set of critical indices
is known, for Theorem 3.4 the solution can be computed as the projection onto
the linear subspace spanned by the vectors in S.
The routine maintains a non-empty subset of Γ called the current set denoted
by S, and a point called the current point denoted by x̄. At each stage of the
routine x̄ ∈ pos(S). When termination occurs the routine either finds the nearest
point in pos(Γ) to u in which case the problem is completely solved or it finds
a critical index of the problem. In the latter case, an LCP of order (n − 1) can
be constructed and the same routine can be applied to this smaller problem.
Hence, the unique solution of the original problem can be obtained after at most
n applications of the routine which finds a critical index.
A characterization useful to find a critical index or the solution to the problem
is provided by the following theorem.
Theorem 3.6. Let x̄ ∈ pos(Γ) be such that 0 ∈ T (u, x̄), where T (u, x̄) is the
tangent hyperplane at x̄ to the ball of center u and radius x̄. If there exists an
index j such that (u − x̄)T γ i ≤ 0 for all i = j and (u − x̄)T γ j then j is a critical
index of N P P (Γ, u)
A characterization of the optimal solution in terms of separating hyperplanes
is given by Robertson et al. [42].
Theorem 3.7. A point x̄ ∈ pos(Γ) is the nearest point in pos(Γ)to y if and only
if 0 ∈ T (y; x̄) and (y − x̄)T γ j ≤ 0, ∀j = 1, ..., n, where T (y, x̂) is the tangent
hyperplane to pos(Γ) in x̂.
The routine to find a critical index alternates distance reduction operations
with line-search and projection steps to find a projection face. In practice, the
routine starts by projecting on the closest edge to the data point. If the optimal-
ity condition is not satisfied, than the procedure iteratively adds vectors to S
and updates the point x̄ while consistently reducing the distance between u and
x̄. The distance reduction operation is carried out efficiently by projecting onto
the non-negative hull of two vectors in Rn , the current point x̄ and a vector γ i
satisfying one of the conditions given by the following theorem.
Theorem 3.8. Given x̄ ∈ pos(Γ), x̄ = 0 such that 0 ∈ T (y, x̄), if for some
i ∈ {1, ..., n} we have (y − x̄)T γ i > 0 and either:
• ||x̄ − y|| ≤ ||Π(y|γ i ) − y|| and {x̄, γ i } is linearly independent, or
• yT γ i ≤ 0
72 M. Dimiccoli
then the projection of y onto the linear hull of {x̄, γ i } is in the relative interior
of pos{x̄, γ i }
Once such updates are no longer possible, it employs a sequence of line-

search steps and projections in the subspace spanned by the vectors in S to find
a projection face of the corresponding pos cone. This line-search is in the same
spirit as the one proposed by Fraser and Massam since the goal is to reduce the
distance to the data point while keeping at the interior of a pos cone. In the
particular case of concave regression, for which m < n, it can be implemented
exactly in the same way.
This algorithm is much faster than the MPDB algorithm. The primary source
of its computational efficiency is that it relies mostly on distance reduction op-
erations and size reduction steps whose computational requirement is relatively
small compared to the computational effort required to find a projection face
through a line-search.
Recently, Liu and Fathy (2011) [20] generalized the work of Murty and Fathy
(1982) to polyhedral non-simplicial cones, hence allowing the set Γ to contain
more than n vectors. What allows the generalization is the equivalence between
the structure of the two problems through the concept of a polar cone.
The authors also proposed several strategies for efficient implementation
mostly based on the mathematical properties of the entities involved. We in-
corporated, where possible, these strategies, to all algorithms we evaluated.
3.2.5. Meyer’s algorithm
Meyer [18] considered the general problem of projecting a data point y ∈ Rn

onto a polyhedral convex cone K ⊂ Rn generated by a finite number m of linear
inequalities. The problem is reduced to find the set of indices Jˆ ⊂ {1, ..., M } ≡ L,
where M ≤ m is the number of linearly independent constraints, corresponding
to not saturated constraints at the solution. Meyer called these indices hinges.
When m ≤ n, the algorithm can be applied to both the primal and the dual
formulation, whereas for m > n it is applied to the dual formulation. In the
following, since for the problem of concave regression m < n, we consider how
to solve the primal problem (2.7).
The search for Jˆ is performed iteratively, starting with an initial guess by
removing or adding one index at time, until the optimal solution is obtained.
For each candidate J k two conditions are tested: the interior point condition
and optimality condition.
The interior point condition is satisfied when the current iterate belongs to
the interior of pos(S), that is when xk ∈ pos(S), where S ⊂ {β i , i ∈ J k } . By
using the following theorem
Theorem 3.9. Let S ⊂ {β i , i = 1..., m}, S = ∅. The optimum solution of

problem (2.7) is in the relative interior of pos(S) if and only if Π(y|span(S)) is
in the relative interior of pos(S).
xk can be computed as projection of y onto the linear hull spanned by the

vectors {β i , i ∈ J k }. If the feasibility condition is not satisfied, the index j ∈
L \ J k corresponding to the most negative coefficient is added to J k and the
interior point condition is checked again.
Once the feasibility condition is satisfied, the optimality condition is tested
by using the characterization given in Theorem 3.7. If it is not satisfied, the
vector β i , i ∈ J k which most violates the condition is removed. The procedure
continues until both conditions are satisfied.
Convergence is guaranteed by the fact that when the algorithm replaces just
one edge, the Sum of the Squared Errors (SSE) after is less than the SSE before
so that the algorithm never produces the same set of edges twice, which would
result in an infinite loop.
In practice, each time a hinge is added, the best solution with n + 1 hinges
where the first n hinges are already given is obtained. But this is not in general
the best fit with n + 1 hinges, so that some hinge may need to be changed.
Therefore, the optimal solution can be interpreted as the best approximation
with the biggest possible number of hinges.
4. Issues about effectivness for large-scale problems
In this section, we discuss strengths and limitations of the algorithms detailed

above in solving large-scale instances of a particular kind of cone regression,
the concave regression. In particular, we consider computational issues related
to numerical stability, computational cost and memory load as well as the suit-
ability to take advantage of a good initialization and to be implemented in an
online fashion.
4.1. Suitability to take advantage of available good estimates
One general strategy for reducing the computational cost of a large-scale opti-
mization problem is to use an initial guess that is easier to calculate and close
to the optimal solution.
Within the class of algorithms with asymptotic convergence, splitting-based
methods work by activating each of the convex constraints repetitively and by
combining them to obtain a sequence converging to a feasible point. Since the
projection point Π(y|K) is characterized by the variational inequality
x̂ = Π(y|K) ∈ K, ∀x ∈ K : y − x̂, x − x̂ ≤ 0 (4.1)
the projection operator Π(·|K) is a closed contraction. Therefore the set of fixed
points of Π(·|K) is exactly K. This prevents the use of an initialization point
belonging to the feasible set as well as the use of multiscale strategies since there
is no guarantee that the solution from a previous level does not belong to the
feasible set.
The same difficulty arises when considering interior point algorithms, since
they need to be initialized to an interior point. In [9], Goldman proved that the
74 M. Dimiccoli
Dykstra’s algorithm can potentially start from better starting values than the
given data point y. The author established the convergence to the nearest point
to the primal cone from an arbitrary point in the intersection of C and the ball
of radius ||y||, where C = {x|x = y − AT λ, λ ≥ 0}, is a rotation of π radiants
of the dual cone Ko with its vertex translated at y. A point satisfying these
conditions can be obtained efficiently by using distance reduction operations
based on Theorem 3.8. It is worthwhile to remark that this result can be easily
interpreted in the active set framework. In fact, the Dykstra’s algorithm can
be undestood as a primal active set method and its solution is a primal dual
feasible point. Therefore, any dual feasible point, that is every point belonging
to the set C, can be used as initialization.
All algorithms with time finite convergence detailed in the previous section
are primal-dual and they reduce the problem of projecting a given data point
onto a convex set to the problem of finding the set of indices corresponding to
not saturated constraints at the solution. In general, they involve the selection
of a subset from a collection of items, say Jˆ ⊆ {1, ..., m}. With this formulation,
they potentially allow us to take advantage of a good estimate of the optimal
active set. However, the adaptation is not rapid since the active set estimate is
updated using one-by-one changes preventing this class of method from being
effective general-purpose solvers for large-scale problems.
In the algorithm of Fraser and Massam, the set of successive approximations
are obtained by moving along a fixed line connecting the initial guess to the
data point. The number of iterations needed to reduce the data point reduces
the number of sectors to be crossed to join the sector containing the data point.
By contrast, in the algorithm of Meyer the proximity of the initial guess to
the data point is not a good criterion selection for the initial guess. In fact,
given an initial guess, the solution is attained by adding and/or removing one
index at time until an optimal solution is found. Taking into account that the
optimal solution can be interpreted as the best approximation with the biggest
possible number of hinges, if the optimal solution contains just a few hinges,
than using the empty set as an initial guess would be much faster than using
the full set of possible hinges. On the other hand, if just a few constraints are
satisfied at equality in the optimal solution, than the full set of indices will be
a much better initial guess than the empty set. Therefore even if the choice of
the initial guess may highly influence the performances, its choice depends on
the data and there is not a well established criterion to fix it.
The Murty and Fathy’s algorithm reduces the size of the problem each time
a critical index is found. Therefore it is not compatible with strategies that
take advantage of a good initial estimate since a good estimate does not lead to
finding a critical index faster.
To overcome the limitation of active set methods in taking advantage of a
good estimate, Curtis et al. [45] have recently proposed an heuristic framework
that allows for multiple simultaneous changes in the active-set estimate, which
often leads to a rapid identification of the optimal set. However, there is no
guarantee of the computational advantages for general problems and, further-
more, the authors recommend their approach when solving generic quadratic
programming problems with many degrees of freedom. That is not the case of
general concave regression problems.
4.1.1. PAV’s inspired approximate solution
To evaluate through numerical simulations the suitability of taking advantage

of good initial estimates, we propose an algorithm inspired to Pool Adjacent
Violators (PAV), whose computational complexity is O(n). Starting from the
original signal, violated constraints are removed one by one by projecting the
current iterate onto the convex cone corresponding to the violated constraint
until a primal feasible solution is obtained. Since the dual feasibility of each
iterate is not guaranteed, the founded solution is not optimal. However, in our
experience the solution is a very good approximation of the optimal solution.
4.2. Suitability to be implemented in an online fashion
Another strategy to deal with large projection problems would be to build and
evaluate the solution incrementally according to the order of its input, as done
in online methodologies developed for dynamic optimization [46] over a stream
of input. Even if the input is given in advance, inputs are processed sequentially
and the algorithm must respond in real-time with no knowledge of future input.
Each new input may cause rising or falling constraints and the final solution is
a sequence of feasible solutions, one for each time step, such that later solutions
build on earlier solutions incrementally.
Of course this strategy requires that the algorithm responds in real-time
for each new input. That would not possible when dealing with large matrix
inverse computations. Let x̂ ∈ Rn be the projection of y ∈ Kn onto K and let
ˆ ∈ Kn+1 be the projection of ȳ ∈ Rn+1 onto Kn+1 . When a new element is
x̄
added to the data point, a new constraint is added too, so that the constraint
cone has a new edge. If this constraint corresponds to a critical index, that is to
a constraint satisfied at equality in the optimal solution, then the projection face
will be the same so that no further computing will be needed. On the contrary,
if the new constraint does not correspond to a critical index, the projection
face will change, including the edge corresponding to the new constraint and
removing and adding some others. Therefore, the major difficulty faced by an
online strategy is the same faced in exploiting good estimates.
4.3. Computational issues
As highlighted in the previous section, despite the different strategies imple-

mented by algorithms with finite-time convergence, generally an index at time
is iteratively added or removed from the current set of indices J k until both
the feasibility condition and the optimality condition are satisfied. Checking
the optimality condition involves computing the inverse of a matrix that differs
76 M. Dimiccoli
slightly from the matrix of the previous iteration. What ’slightly’ exactly means
depends on the specific algorithm and it is detailed in the following.
The algorithm of Fraser and Massam involves the computation of an n × n
fixed size matrix inverse at each iteration. The matrix to be inverted differs from
the matrix of the previous iteration only for one column, being the change of
the form A → (A + u × v), where u is a unitary vector with only one nonzero
component corresponding to the index of the column to be changed, and v
corresponds to the difference between the elements of the dual basis vector and
the elements of the primal basis vector or vice versa. In this case the inverse can
be updated efficiently by using the Sherman-Morrison formula: (A + u × v)−1 =
−1
z×w
1+λ , where z = A u, w = (A−1 )T v, λ = v T z. Therefore only two matrix
multiplications and a scalar product are needed to compute the new inverse at
each step.
The algorithm of Meyer, as well as the one of Liu and Fathi [20] involve the
computation of an inverse matrix of variable size at each iteration. The matrix to
be inverted differs from the matrix of the previous iteration only for one column,
which has been added or removed. Since the matrix to be inverted A(J) ∈ Rr×n ,
with r ≤ m is generally rectangular, the Moore-Penrose generalized inverse or
pseudoinverse is computed:A(J)† = A(J)T (A(J)A(J)T )−1 .
In Matlab and Scilab, the computation of the pseudoinverse is based on the
Singular Value Decomposition (SVD) and singular values lower than a tolerance
are treated as zero. The advantage of this approach is that the pseudoinverse
can be computed also for a nonsingular matrix. However, the method proposed
by Liu and Fathi [20] to improve the computational efficiency of their algorithm
does not take advantage of the SVD approach since it consists in updating the
matrix (A(J)A(J)T )−1 . If the matrix A(J) is ill-conditioned, then the inverse
cannot be computed with good accuracy and for the matrix A(J)AT (J) is even
more so because this operation squares the condition number of the matrix
A(J).
A better solution would be to update directly the pseudoinverse. This can
be achieved when a column is added
byusing the method proposed in [47] Let
AT ∈ Rn×r , x ∈ Rn and B = AT x ∈ Rn×(r+1) . The pseudoinverse of B
can be computed from the pseudoinverse of AT as follows.
†
A − A† xw† wT
B† = † , where w = (I − AA† )x and w† = .
w ||w||2
Alternatively, the transformation A(J)T (A(J)A(J)T )−1 y canbe efficiently

R11
computed by a QR decomposition approach. Let AT = Q11 Q12 be the
0
QR decomposition of AT , where R11 is an m × m invertible upper triangular
matrix. Then: A(J)T (A(J)A(J)T )−1 = Q11 (R1T 1)−1 . The inverse of an upper
triangular matrix can be efficiently implemented by a left matrix division or by
more sophisticated methods as the one proposed in [48].
Courrieu [49] proposed a method based on the full rank Cholesky decomposi-
tion which has a computation time substantially shorter of the method based on
SVD decomposition. The two main operations on which his method is based are
the full rank Cholesky factorization of AT A and the inversion of LT L, where L
is a full rank matrix . On a serial processor these computations are of complexity
order O(n3 ) and O(m3 ) respectively. However, in a parallel architecture, with
as many processors as necessary, the time complexity for Cholesky factorization
of AT A could reduce to O(n), while the time complexity for the inversion of the
symmetric positive definite matrix LT L could reduce to O(log(r)). However, the
computational advantage of this method can be appreciated only when r << n,
since the inverse of a r × r matrix has to be computed, which is not in general
the case, specially for concave regression problems.
The method proposed by Zhu and Li [51] for recursive constrained least
squares problems, found that the exact solution of Linear Equality-constrained
Least Squares can be obtained by the same recursion as for the unconstrained
problem, provided that the Rescricted Least Squares procedure is appropriately
initialized. However, even this approach offers significant advantages in term of
computational cost only when the number of constraints m is small, which is
not the case for large-scale concave regression problems.
5. Improving the active set framework for concave regression

problems
An active set algorithm for solving the concave regression problem generates
a sequence of quasi stationary points. A primal (dual) feasible active set al-
gorithm generates a sequence of primal (dual) feasible quasi stationary points
with decreasing objective function values and terminates when the dual (primal)
feasibility is satisfied. An active set J induces a unique partition of the indices
{1, ..., n} into blocks. A block B of such partition is a set of consecutive integers,
say, {p, p+1, ..., q}, such that the index i of the constraint xi+1 −xi ≤ xi+2 −xi+1
is in J for each i such that p ≤ i ≤ q − 2. Conversly any such partition of indices
determines a unique active set. Denoting by λi the multiplier associated with
the ith constraint, the KTT optimality conditions can be written as:
⎧
⎪
⎪ xi+1 − xi ≤ xi+2 − xi+1 i = 1, ..., n − 2
⎪
⎪
⎪
⎪ y1 − x1 = λ1
⎪
⎪
⎪
⎪ y2 − x2 = −2λ1 + λ2
⎪
⎪
⎨ y3 − x3 = λ1 − 2λ2 + λ3
yn−1 − xn−1 = λn−4 − 2λn−3 + λn−2
⎪
⎪
⎪
⎪ yn−1 − xn−1 = λn−3 − 2λn−2
⎪
⎪
⎪
⎪ yn − xn = λn−2
⎪
⎪
⎪
⎪ λi ≥ 0 i = 1, ..., n − 2
⎩
λi (2xi+1 − xi − xi+2 ) = 0 i = 1, ..., n − 2
n n
It is easy to show that: i=1 i(yi − xi ) = 0 and i=1 (yi − xi ) = 0. Knowing
that each block can be represented by an affine function xi = α + βi, the above
78 M. Dimiccoli
system can be written as:

⎧ n1 n 1
⎪
⎪ i=1 yi = α1 i=1 1
ni + β n1 + λn1
⎪
⎪ n
⎨ i=n1 +1 yi = α i=n1 +1 (i − n1 ) + β 2 (n − n1 ) − λn1
2
n n
(in − n1 ) + β (n2− n12)
− λn1
2 2
i=n1 +1 yi = α
⎪
⎪ n
i=n1 +1
⎪
⎪ (i − ni )yi = α 2
(i − n 1 ) + β i=n1 +1n (i − n1 ) − λn1
⎩ 1 i=n1 +1 1 i=n1 +1
λ n1 + β − α2 − β 2 = 0
Therefore, for each block the unknown variables to be computed are α, β, λ.

The systems to be solved
n can bewritten as Ax = b, where x = (α1 , β 1 , α2 , β 2 ,
n1 n n
λn1 ), b = ( i=1 yi , i=n1 +1 yi , i=n1 +1 yi , i=n1 +1 (i − ni )yi , 0) and
⎛ n1 ⎞
i=1 i n1 0 0 1
⎜ 0 n−n1
⎜ 0 i n − n1 −1⎟⎟
⎜n1 2 n 1 i=1
⎟
A=⎜
⎜ i=1
i i=1 i 0 0 n1 ⎟
⎟
⎜ 0 n2 −n1 n2 −n1
⎝ 0 i=1 i2 i=1 i 0⎟⎠
n1 + 1 1 −1 −1 0
In general, if k is the number of blocks, than the system to be solved has

size 3k − 1. As observed in [52], in the case of concave regression the number
of blocks at the solution is usually much lower than n. Therefore, an active
set algorithm that starts with an empty active set, should found the solution
without the need of inverting large matrices.
6. Comparative results
In Table 1 we compare qualitatively the algorithms analyzed above in terms of

their formulation (primal, dual, or primal-dual), their suitability to take advan-
tage of a good initialization and their major limitations in dealing with large
scale problems. All analyzed methods are dual or primal-dual. Dual methods
cannot be initialized, whereas initialization in primal-dual methods is allowed
but constrained. The major limitation of asymptotic convergence methods when
dealing with large problems is the convergence, whereas finite-time convergence
methods become too slow and numerically unstable because of accumulation
errors.
In the following, we provide evidence of these limitations through numer-
ical simulations reported in the Appendix A and we compare quantitatively
their performances in term of distance from the solution for one second, one
minute and ten minutes of CPU time. Instead of using only unstructured ran-
dom data for testing, we considered signals of increasing difficulty level vary-
ing from a noisy concave signal to a concave/convex noisy signal. More pre-
cisely, we considered three signals, whose equations and plots are given in Fig-
ure 3 of the Appendix A and, for each of them, we considered three different
sizes: n ∈ {50, 500, 1000} and three increasing levels of noise (standard devia-
tions σ ∈ {0.01, 0.1, 0.5}). To evaluate robustness against initialization varia-
Table 1
Qualitative comparison of all algorithms
Algorithm Type Formulation Initialization Limitation
Hildreth dual 2.9 not possible convergence

Dykstra primal-dual 2.7 constrained convergence
ADMM dual 3.4 not possible convergence
LSPS primal 3.1 constrained convergence
Uzawa dual 3.3 not possible convergence
MPDB primal-dual 3.9 constrained slow or unstable
Meyer primal-dual 2.7 constrained and not robust slow or unstable
Critical Index dual 3.10 not possible slow or unstable
tion for finite-time active set methods, we considered three different initializa-
tions of the active set J: the empty set, the set of m indexes corresponding
to the linear inequality constraints, and the set of not saturated constraints
obtained by using the algorithm inspired to PAV described in the previous sec-
tion.
In the class of algorithms with asymptotic convergence, the most efficient
appears to be ADMM. This is evident even for very small size data, when
using difficult signals such as near-convex signals or very noisy signals (see Fig-
ure 5 and Figure 6 of the Appendix A). For noisy concave signals (Figure 4
of the Appendix A), the computational efficiency of ADMM is more evident
in the presence of a high level of noise. Already for signals of size 500 the
performance of ADMM is not very good; convergence (SSE < 10−6 ) is not
completely attained. The algorithm of Meyer is very sensitive to the initial-
ization. It gives good performances when the signal is a Gaussian white noise
and the initial active set is empty since most of constraints are expected to
be saturated at the solution (see Figure 6 and Figure 9 of the Appendix A).
Given a good initialization, the MPDB algorithm allows to compute the ex-
act solution faster than other methods. However, this algorithm is not nu-
merically stable since it requires one to compute the inverse of a matrix at
each iteration. As explained in section 4.3, this is done incrementally by using
the Shermann-Morrison formula. However, numerical round-off errors cumulate
so that, in our implementation, the exact inverse is computed each 150 iter-
ations. As it can be observed in Figure 10, Figure 11 and Figure 12 of the
Appendix A, that refer to signals of size 1000, sometimes the round-off error
dominates the calculation and the distance to the solution increases instead of
decreasing.
These results demonstrated that, although there have been theoretical and
practical advances in recent years, the use of shape-constrained nonparametric
techniques for solving large scale problems (more than many thousand of points)
is still limited by computational issues. To deal with very large scale problems,
up to a million of points, matrix inverse calculation should be avoided and more
efforts should be devoted to find a way to better initialize primal-dual methods,
by computing approximate solutions and exploiting multi-scale strategies.
80 M. Dimiccoli
7. Conclusions
In this paper we have stated, analyzed and compared qualitatively and quantita-
tively several optimization approaches for solving the cone regression problem.
We have distinguished two broad classes of methods. On the one side, methods
with asymptotic convergence that rest on the properties of proximity operators
and act by finding the solution as the limit of a sequence of successive approx-
imations. On the other side, methods with finite-time convergence that exploit
the geometric properties of polyhedral convex cones and find the exact solution
as non-negative linear combination of functions, forming a basis in a specified
finite dimensional space.
Simulations up to one thousand of points have demonstrated that the choice
of the optimization approach severely impacts algorithmic performances. In par-
ticular, it has emerged that methods with finite-time convergence are much more
efficient with respect to asymptotic-convergence methods. However, from this
study it emerged that they face a twofold difficulty to cope with large-scale
optimization. The first difficulty arises from the fact that all algorithms of this
class modify the active set estimate one-by-one, making the adaptation of the
initial active set estimation very slow; the second difficulty lies in the fact that
they involve the computation of the inverse or the pseudoinverse of a matrix
at each variation of the active set. Although there exists many methods to do
that efficiently when the matrix rank is much lower that the size of the data,
this condition cannot be assured in general. Incremental strategies to reduce the
cost of computing the inverse of a matrix when the inverse of a slightly different
matrix is known, are bounded by round-off error accumulation in an iterative
setting.
The results of this study suggest that to be able to treat very large-scale prob-
lems (up to a million of points) further research should focus on finding a way to
exploit classical multi-scale strategies and to compute an approximate solution
through a penalization or splitting method without involving any matrix inverse
calculation.
Appendix A
A.1. Signals used for comparative evaluations
⎧
⎨ 2nsin( 5n
24
z) z = 1, ..., n3
S1 (z) = α + βz z = n3 + 1, ..., 2n
⎩ 3
3
γz + δ z = 2n
3 + 1, ..., n
−2 3
where β = 1
10 , α = 2nsin( 85 ) − βn, γ = n2 3 − γ 27 .
and δ = α + β 2n 8n
S2 (x) = N (μ, σ 2 )
6
S3 (z) = sinc( z − 1)
n
Fig 3. Signals S1 (top left) and S2 (top rigth) to which has been added white noise with three
different values of standard deviation σ = 0.01, σ = 0.1, and σ = 0.5.
82 M. Dimiccoli
A.2. Comparative evaluations on signals of size 50
Fig 4. Distance to the solution (L2 norm) for a signal of type S1 of size 50. From top to
the bottom, three increasing level of noise have been added. (Left) Algorithms with finite-
time convergence: all them converge instantaneously. (Right) Algorithms with asymptotic
convergence: ADMM is more robust to noise. LSPS and Dykstra use a parametrization in the
primal parameter space which causes numerical round-off errors to cumulate so that the fitted
values do not satisfy the constraints of the dual problem and convergence is not attained.
the bottom, three increasing level of noise have been added. (Left) Algorithms with finite-
time convergence: all them converge instantaneously. (Right) Algorithms with asymptotic
convergence: both ADMM and Hildreth attain the solution when no noise is added but ADMM
is much more robust to noise.
84 M. Dimiccoli
Fig 6. Distance to the solution (L2 norm) for Gaussian white noise signals of size 50. (Left)
Algorithms with finite-time convergence: all them converge instantaneously. (Right) Algo-
rithms with asymptotic convergence: ADMM is the only algorithm that converges.
Fig 7. Distance to the solution (L2 norm) for a signal of type S1 of size 500. From top to the
bottom, three increasing levels of noise have been added. (Left) Algorithms with finite-time
convergence: Meyer’s algorithm is very sensitive to the initialization. The best performance
is achieved by MPDB with a PAV approximate initialization. (Right) Algorithms with asymp-
totic convergence: ADMM converges only when the signal is slighly noisy.
86 M. Dimiccoli
the bottom, three increasing level of noise have been added. (Left) Algorithms with finite-time
convergence: the best performance is achieved by MPDB with a PAV’s inspired initialization.
(Right) Algorithms with asymptotic convergence: no algorithm converges.
Fig 9. Distance to the solution (L2 norm) for a Gaussian white noise signal of size 500. (Left)
Algorithms with finite-time convergence: the best performance is achieved by MPDB with a
PAV’s inspired initialization. (Right) Algorithms with asymptotic convergence: no algorithm
converges.
88 M. Dimiccoli
90 M. Dimiccoli
Fig 12. Distance to the solution (L2 norm) for a signal of type S3 of size 1000. (Left)
Algorithms with finite-time convergence: the best performance is achieved by MPDB with a
PAV’s inspired initialization when the round-off error does not dominate the calculation and
by the Meyer algorithm whose initial active set is empy. This is easy to understand since, for
a pure noise signal a high degree of freedom is expected at the solution and therefore the empty
active set is close to the final active set. (Right) Algorithms with asymptotic convergence: no
algorithm converges.
Appendix B: Pseudocode of the algorithms reviewed
Basic notations:
y is the input signal
N is the iteration number (when required by the algorithm)
x̂ is the output signal
Algorithm 1: Hildreth’s algorithm

Data: y ∈ Rn , A ∈ Rm×n
Result: x̂ = argmin ||y − x||2 subject to Ax ≤ 0
begin
P ←− AAT
b = Ay
λ0i ←− 0 ∀i = 1, ..., m
for k = 1 : N do
for i = 1 : m do
i−1
λk+1
i = max(0, (bi − j=1 pij λkj − m k−1
j=i+1 pij λj )/pii )
x̂ ←− y − AT λN
end
Algorithm 2: Dykstra (proximal formulation)

Data: y ∈ Rn , A ∈ Rm×n , maximum number of iterations N
begin
x ←− y
Ii ←− 0 ∀i = 1, ..., m
for k = 1 : N do
for i = 1 : m do
z = x + Ii
Ii ←− proxAi z
x ←− z − Ii
end
AT T T
i (z Ai )
proxAi z = prox{x|Ai x≤0} z = ||Ai ||2
92 M. Dimiccoli
Algorithm 3: Alternating Directions Method of Multipliers (ADMM)

Result: x̂ = argminAx≤0 ||y − x||2 , is equivalent to x̂ = argminAx=z,z≤0 f (x) + g(z),
f (x) = ||y − x||2 , g(z) = f ◦ A = 12 ||Ax − z||2
begin
v 1 ←− 0
λ1 ←− 0
τ >0
for k = 1 : N do
x ←− proxA f (v + λ)
z ←− Ax
v ←− proxg (z − λ)
λ ←− λ + τ (v − z)
end
−1 prox (z − λ) = min
proxA z≤0 z − λ/τ
T T
f (v + λ) = RA (λ + τ v) + Ry, R = (In + τ A A) g
Algorithm 4: Least Squares in a Product Space (LSPS)

Data: y ∈ Rn , A ∈ Rm×n , maximum number of iterations N , λ > 0
Result: x̂ = argmin ||y − x||2 subject to Ax ≤ 0, which is equivalent to requires that
x ∈ S ∩ D, where D is the diagonal vector subspace, S is a closed convex set
corresponding to the cartesian product of the sets (Ai )
begin
x1 ←− y
d1i ←− y, i = 1, ..., m
for k = 1 : N do
xk+1 ←− λ(proxA1 dk1 , ..., proxAm dkm )
k+1 k+1
dk+1 ←− ( m i=1 xi , ..., mi=1 xi )
end
AT T T
i (di Ai )
proxAi di = prox{x|Ai x≤0} di = ||Ai ||2
Algorithm 5: Uzawa
Result: x̂ = argmin ||y − x||2 subject to Ax ≤ 0, which is equivalent to 3.3
begin
ρ≥0
k ←− 1
for k = 1 : N do
xk = proxg y
μk+1 = maxμ≥0 μk + ρAxk
end
g(x) =< μ, Ax >
proxg y = argminx 12 ||x − y||2 + g(μ)
Algorithm 6: Mixed Primal-Dual Bases (MPDB)

Data: (z, y) ∈ Rn × Rn , where z is the vector of independent variables, A ∈ Rm×n .
begin
1 ∈ Rn , 1 = (1, ..., 1)
γm+1←− 1/||1||
n
mi z i
z̄ = i=1
n
i=1 mi
γm+2 ←− (z − z̄1)/||z − z̄1||
V = [γm+1 , γm+2 ]
Γ = [AT T
1 , ..., Am , V ]
J ←− {1, ..., m}
v ←− V V † y
u ←− y − v
Δ ←− −(Γ T −1
m )
x0 ←− i=1 Δi
dx ←− Ax0
dy ←− Au
stop ←− F alse
while stop
= True do
if dyi ≥ 0, i = 1, ..., m then
stop ←− True
else
t ←− mini,dyi <0 dx dx i
i −dyi
x0 ←− x0 + t(u − x0 )
dxi
j ←− argmini mini,dyi <0 dxi −dyi
if j ∈ J then
J ←− J \ {j}
else
J ←− J ∪ {j}
dxj ←− 0
C = AT
C(J) = Δ(J)
dy ←− C † u
x̂ ←− Δ(J)Δ(J)† u + v
end
Δ(J) = {δj , j ∈ J}
94 M. Dimiccoli
Algorithm 7: Meyer’s Algorithm

Data: y ∈ Rn , A ∈ Rm×n , J ⊆ {1, ..., m}.
begin
stop ←− F alse
b ←− Ay
if j = argmaxi bi , bj > then
J ←− J ∪ {j}
else
stop ←− True
θi ←− 0 ∀i = 1, ..., m
while stop
= True do
a ←− A(J)† y
if j = argmini ai , aj < then
J ←− J \ {j}
stop ←− True
θi ←− 0 ∀i = 1, ..., m
else
θ ←− A(J)T a
b ←− A(y − θ)
if j = argmaxi bi , bj > , j ∈
/ J then
J ←− J ∪ {j}
stop ←− False
else
stop ←− True
x̂ ←− y − θ
end
A(J) = {Aj , j ∈ J}
Algorithm 8: Active index algorithm

Data: y ∈ Rn , A ∈ Rm×n , J ⊆ {1, ..., m}.
begin
n n
γ m+1 ←− 1/||1||, z̄ ←− i=1 mi zi / i=1 mi , γ m+2 ←− (z − z̄1)/||z − z̄1||
E ←− [γm+1 , γm+2 ], Γ ←− [AT 1 , ..., A T , E]
m
p ←− E E † y, Y ←− y − p
critical ←− F alse, solution ←− F alse, stop ←− F alse
V criticals ←− ∅, N = n
while stop == F alse do
if critical == T rue then
V ←− γ k , c ←− k
N critical ←− card(V criticals)
Jk(i) ←− T rue ∀i = 1, ..., N
if V criticals == 0 then
Jk(V criticals) ←− F alse
J ←− {i : Jk(i) == T rue}, k ←− J(k)
N critical ←− N critical + 1
V criticals(N critical) ←− k
y ←− y − v(y v)/||v||
for i = 1 : m do
γ i ←− γ i − (vv γ i ))/||v||
A ←− A \ {Ac }
critical ←− F alse
else if solution == T rue then
x ←− y − x + p, stop ←− T rue
else
for i = 1 : n do
if y T γi ≤ 0 then
Vi = 0
else if y T γı > 0 then
V i ←− γi (y T γi )/||γi ||2
l ←− argmini∈{1,...,n} {||V i − y||}, x̄ ←− V l

N (x̄) = {i ∈ {1, ..., n} : (y − x̄)T γi ) > 0}
if card(N (x̄)) == 0 then
stop ←− T rue
else if card(N (x̄)) == 1 then
î ←− N (x̄), critical ←− T rue
else
g ∈ N (x̄)
S ←− {γl , γg }
I(S) ←− {l, g}
ŷ ←− Γ(S)Γ(S)† y
x̄ ←− ŷ
N (x̄) ←− {i ∈ {1, ..., n} : (y − x̄)T γi ) > 0}
stop1 ←− F alse
while stop1 == F alse do
Ī(S) ←− {1, ..., n}\(S)
if card(N (x̄)) == 0 then
stop1 ←− T rue, stop ←− T rue
else if card(N (x̄)) == 1 then
î ←− N (x̄), stop1 ←− T rue
else if N (x̄) ∩ Ī(S) then
ŷ ←− Γ(S)Γ(S)† y
S ←− S ∪ {γg }, I(S) ←− I(S) ∪ {g}
x̄ ←− ŷ
else
y(S) ←− Γ(S)(Γ(S)T Γ(S))−1 Γ(S)T y
if y(S) ∈ pos(S) then
x̄ ←− y(S)

x̄ ←− j∈I(S) αj γ j

y(S) ←− j∈I(S) γj γj
αj
λ̄ ←− argminj:γ <0 { }
j αj −γj
x̄ ←− (1 − λ̄)x̄ + λ̄y(S)
S ←− S\{γk }
I(S) ←− I(S)\{k}
end
Γ(J) = {Γj , j ∈ J}
96 M. Dimiccoli
Acknowledgements
The author would like to thank Lionel Moisan for insighful suggestions. This
work has been carried out during a postdoctoral stage in the Laboratory of
Applied Mathematics (MAP5, CNRS UMR 8145) at Paris Descartes University.
Supplementary Material
Implementation in Scilab and Matlab of the algorithms reviewed in

“Fundamentals of cone regression”
(doi: 10.1214/16-SS114SUPP; .zip).
References
[1] C. Hildreth, Point Estimates of Ordinates of Concave Functions, Journal of

the American Statistical Association 49 (267), 598–619, (1954). MR0065093
[2] W. Dent, A Note on Least Squares Fitting of Functions Constrained to
be Either non-negative, Nondecreasing or Convex, Management Science
20 (1), 130–132, (1973). MR0337328
[3] C. A. Holloway, On the Estimation of Convex Functions, Operations Re-
search 27 (2), 401–407, (1979).
[4] C. Perkins, A Convergence Analysis of Dykstra’s Algorithm for Polyhe-
dral Sets, SIAM Journal on Numerical Analysis 40 (2), 792–804, (2003).
MR1921678
[5] N. I. Gould, How good are projection methods for convex feasibility prob-
lems?, Comput. Optim. Appl. 40 (1), 1–12, (2008). MR2386586
[6] Y. Censor, W. Chen, P. L. Combettes, R. Davidi, G. T. Herman, On the
Effectiveness of Projection Methods for Convex Feasibility Problems with
Linear Inequality Constraints, Computational Optimization and Applica-
tions, 51 (3), 1065–1088, (2012). MR2891928
[7] J. J. Moreau, Decomposition orthogonale d’un espace Hilbertien selon deux
cones mutuellement polaires, Comptes Rendus de l’Académie des Sciences
255, 238–240, (1962). MR0139919
[8] H. W. Kuhn, A. W. Tucker, Nonlinear programming, in: Proceedings of the
Second Berkeley Symposium on Mathematical Statistics and Probability,
481–492, (1951). MR0047303
[9] S. M. Goldman, P. A. Ruud, Nonparametric Multivariate Regression Sub-
ject to Constraint, Preprint. Department of Economics, UCB, (1993).
[10] J. J. Moreau, Fonctions convexes duales et points proximaux dans un espace
Hilbertien, Comptes Rendus de l’Académie des Sciences (Paris), Série A
255, 2897–2899, (1962). MR0144188
[11] J. J. Moreau, Propriétées des applications prox , Comptes Rendus
de l’Académie des Sciences (Paris), Série A 256, 1069–1071, (1963).
MR0149244
[12] R. Dykstra, An Algorithm for Restricted Least Squares Regression, Jour-

nal of the American Statistical Association 78 (384), 837–842, (1983).
MR0727568
[13] D. R. Wilhelmsen, A Nearest Point Algorithm for Convex Polyhedral Cones
and Applications to Positive Linear Approximation, Mathematics of Com-
putation 30 (133), (1976). MR0394439
[14] B. N. Pshenichny, M. Y. Danilin, Numerical methods in extremal problems,
Mir, Moscow, (1978). MR0502886
[15] C. F. Wu, Some algorithms for concave and isotonic regression, Studies in
the Management Sciences 19, 105–116, (1982). MR0723345
[16] D. A. S. Fraser, H. Massam, A Mixed Primal-Dual Bases Algorithm for Re-
gression under Inequality Constraints. Application to Concave Regression,
Scandinavian Journal of Statistics 16 (1), 65–74, (1989) . MR1003969
[17] M. C. Meyer, An extension of the mixed primal-dual bases algorithm to the
case of more constraints than dimensions, Journal of Statistical Planning
and Inference, 13–31, (1999). MR1718405
[18] M. C. Meyer, A Simple New Algorithm for Quadratic Programming with
Applications in Statistics, Communications in Statistics.
[19] K. G. Murty, Y. Fathi, A critical index algorithm for nearest point problems
on simplicial cones, Mathematical Programming 23 (1), 206–215, (1982).
MR0657080
[20] Z. Liu, Y. Fathi, An active index algorithm for the nearest point problem in
a polyhedral cone, Computational Optimimization and Applications 49 (3),
435–456, (2011). MR2803860
[21] W. Kahan, Gauss-Seidel Methods of Solving Large Systems of Linear Equa-
tions, Ph.D. thesis, University of Toronto (1958). MR2702985
[22] J. von Neumann, Functional Operators, Volume II, Princeton University
Press, USA, (1950).
[23] J. P. Boyle, R. L. Dykstra, A method for finding projections onto the in-
tersection of convex sets in Hilbert spaces, Lecture Notes in Statistics 37,
28–47, (1986). MR0875647
[24] H. R. Varian, The Nonparametric Approach to Production Analysis, Econo-
metrica 52 (3), 579–597, (1984). MR0740302
[25] A. N. Iusem, A. R. Pierro, On the convergence of Han’s method for convex
programming with quadratic objective, Mathematical Programming 52 (1),
265–284, (1991). MR1126172
[26] G. Crombez, Finding projections onto the intersection of convex sets in
Hilbert spaces, Numerical Functional Analysis and Optimization 16 (5–6),
637–652, (1995). MR1341103
[27] N. Gaffke, R. Mathar, A cyclic projection algorithm via duality, Metrika
36 (1), 29–54, (1989). MR0985010
[28] H. H. Bauschke, J. M. Borwein, A. S. Lewis, On the method of cyclic
projections for convex sets in Hilbert space, Preprint, (1994). MR1442992
[29] P. L. Combettes, J. C. Pesquet, Proximal Splitting Methods in Signal Pro-
cessing, Fixed-point algorithms for inverse problems in science and engi-
neering, Springer New York, NY, 185–212, (2011). MR2858838
98 M. Dimiccoli
[30] S. P. Han, A successive projection method, Mathematical Programming,

40 (1), 1–14, (1988). MR0923692
[31] P. A. Ruud, Restricted least squares subject to monotonicity and concav-
ity constraints, in: Advances in Economics and Econometrics: Theory and
Applications, Econometric Society Monographs, 3, 166–187, (1997).
[32] M. Hestenes, Multiplier and gradient methods, Journal of Optimization
Theory and Applications 4 (5), 303–320, (1969). MR0271809
[33] M. J. D. Powell, A method for nonlinear constraints in minimization prob-
lems, R. Fletcher Editions, Optimization, Academic Press, New York, 283–
298, (1969). MR0272403
[34] J. Eckstein, D. P. Bertsekas, On the Douglas-Rachford splitting method
and the proximal point algorithm for maximal monotone operators, Math-
ematical Programming, 55 (3), 293–318, (1992). MR1168183
[35] G. Pierra, Decomposition through formalization in a product space, Math-
ematical Programming, 28 (1), 96–115, (1984). MR0727421
[36] I. I. Eremin, Fejer mappings and problems of convex optimization, Sibirskii
Matematicheskii Zhurnal, 10, 1034–1047, (1969). MR0253750
[37] K. J. Arrow, L. Hurwicz, H. Uzawa, Studies in Linear and Non-Linear
Programming, Cambridge University Press, (1958). MR0108399
[38] M. J. Silvapulle, P. K. Sen, Constrained statistical inference: Order, in-
equality, and shape constraints, John Wiley & Sons, (2011). MR2099529
[39] R. T. Rockafellar, Convex Analysis (Princeton Landmarks in Mathematics
and Physics), Princeton University Press, (1970). MR1451876
[40] E. H. Zarantonello (Ed.), Projections on Convex Sets in Hilbert Space and
Spectral Theory, Academic Press, (1971).
[41] K. G. Murty, Linear complementarity, linear and nonlinear programming,
Berlin: Heldermann Verlag, (1988). MR0949214
[42] T. Robertson, F. Wright, R. L. Dykstra, T. Robertson, Order restricted
statistical inference, Statistical Papers, 30 (1), 316–316, (1989). MR0961262
[43] N. Karmarkar, A new polynomial-time algorithm for linear program-
ming, Proceedings of the sixteenth annual ACM symposium on Theory
of computing, STOC ’84, ACM, New York, NY, USA, 302–311, (1984).
MR0779900
[44] J. N. Singh, D. Singh, Interior-Point Methods for Linear Programming: A
Review, International Journal of Mathematical Education in Science and
Technology, 33 (3), 405–423, (2002). MR1912442
[45] F. Curtis, Z. Han, D. Robinson, A Globally Convergent Primal-Dual
Active-Set Framework for Large-Scale Convex Quadratic Optimization,
Computational Optimization and Applications, 60 (2): 311–341, (2015).
MR3316681
[46] T. Bhatia, L. T. Biegler, Dynamic Optimization in the Design and Schedul-
ing of Multiproduct Batch Plants, Industrial and Engineering Chemistry
Research. 35 (7), 2234–2246, (1996).
[47] E. Andelić, M. Schafföner, M. Katz, S. E. Krüger, A. Wendemuth, Kernel
least-squares models using updates of the pseudoinverse, Neural Computa-
tions, 18 (12), 2928–2935, (2006).
[48] R. Mahfoudhi, A fast triangular matrix inversion, Proceedings of the World

Congress on Engineering, (2012).
[49] P. Courrieu, Fast computation of Moore-Penrose inverse matrices, Neural
Information Processing - Letters and Reviews 8 (2), 25–29, (2008).
[50] B.-Israel Adi, A Newton-Raphson method for the solution of systems of
equations, Journal of Mathematical Analysis and Applications, 15 (2), 243–
252, (1966). MR0205445
[51] Y. Zhu, X. Li, Recursive least squares with linear contraints, Communica-
tions in Information and Systems 7 (3), 287–312, (2007). MR2516245
[52] T. Kuosmanen, Representation Theorem for Convex Nonparametric Least
Squares, The Econometrics Journal, 11 (2), 308–325, (2008).

Fundamentals of Cone Regression: Mariella Dimiccoli

Uploaded by

Copyright:

Available Formats

Fundamentals of Cone Regression: Mariella Dimiccoli

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Fundamentals of Cone Regression: Mariella Dimiccoli

Uploaded by

Copyright:

Available Formats

Statistics Surveys

Vol. 10 (2016) 53–99

Fundamentals of cone regression

Abstract: Cone regression is a particular case of quadratic programming

MSC 2010 subject classiﬁcations: Primary 62; secondary 90.

Received March 2015.

3.1.6 Alternating Direction Method of Multipliers (ADMM) . . 64

Cone regression analysis is a valuable alternative to more traditional parametric-

2. Statement of the problem, basic notations and basic facts

The aim of a regression analysis is to produce a reasonable analysis to the

where z ∈ R is the explanatory (dependent) variable, y ∈ Rd is the response

yit = f (zi ) + it , i = 1, ..., n t = 1, ..., Ti (2.2)

2.1. Convex quadratic programming (CQP) formulation

2.1.1. Primal formulation

2.1.2. Dual formulation

The Moreau decomposition theorem is as follows.

(i) x̂ = argmin ||x − y||2 , x̂o = argmin ||x − y||2

By relying on this theorem we can alternatively solve problem (2.7) by ﬁrst

Intuitively speaking, the edges of a polyhedral convex cone are one-dimen-

λ̂ = argmin(y − AT λ)T W (y − AT λ) (2.9)

x̂ = argmin ||x||2 (2.10)

where C = {x|x = y − AT λ, λ ≥ 0} is a rotation of the dual cone with its vertex

2.2. Linear complementarity problem (LCP) formulation

The CQP (2.7) can be recast as a linear complementarity problem (LCP). To

2.3. Proximal formulation

The proximity operator is characterized by the property

∀(x, p) ∈ Rn × Rn , p = proxf (y) ⇐⇒ y − p ∈ ∂f (p), (2.19)

where ∂f : Rn → 2R is the subdiﬀerential of f :

The proximity operator of a convex function is a generalization of the projection

it results that proxıC = Π(·|C).

x̂ = argmin ||x − y||2 = argmin{||x − y||2 + ıK (x)} (2.21)

Alternatively, using the fact that (K1 ∩ ... ∩ Km )o = K1o + ... + Km

Intuitively, the base operation of a proximal algorithm is evaluating the prox-

3. State of the art

3.1. Algorithms with asymptotic convergence

3.1.1. Least squares in a product space (LSPS)

xk+1 = xk + λk (PD ◦ PKm (xk ) − xk ) (3.2)

3.1.2. Uzawa method

A classical method to solve a convex minimization problem subject to inequal-

3.1.3. Hildreth’s algorithm

3.1.4. Primal-dual interior point methods

First introduced by Karmakar in 1984 [43], primal-dual interior point methods

The perturbed complementarity condition introduces a nonlinearity, there-

3.1.5. Dykstra’s algorithm

In 1983, Dykstra [12] proposed a generalization of the Hildreth’s procedure

3.1.6. Alternating Direction Method of Multipliers (ADMM)

ADMM is an augmented Lagrangian technique [32, 33] which can be applied to

F ind argmin ||y − x||2 + g(Ax) (3.4)

where f (x) = ||y − x||2 . Denoting by proxA

proxg = proxf ◦A the implementation detailed in the Appendix B is obtained.

3.2. Algorithms with ﬁnite-time convergence

3.2.1. Properties of polyhedral convex cones with m ≤ n

Deﬁnition 4. Let K be a polyhedral convex cone in Rn , then F ⊆ K is a face

The set of all relatively open faces

u + v. Since v can be easily computed as v = X(X T X)−1 X T y with X =

Figure 2 illustrates this result for span(K) = R2 . Points in SJ project onto

to K project on themself, while each other vector of Rn project onto a unique

3.2.2. Early algorithms based on the properties of polyhedral convex cones

Theorem 3.4. Let S ⊂ Γ, S = 0. The optimum solution of N P P [S; u] is in

Theorem 3.9. Let S ⊂ {β i , i = 1..., m}, S = ∅. The optimum solution of

x̂ = Π(y|K) ∈ K, ∀x ∈ K : y − x̂, x − x̂ ≤ 0 (4.1)

Alternatively, the transformation A(J)T (A(J)A(J)T )−1 y canbe eﬃciently