Projection on Polyhedral Cone

From Wikimization

Jump to: navigation, search

This is an open problem in Convex Optimization. At first glance, it seems rather simple; the problem is certainly easily understood:

We simply want a formula for projecting a given point in Euclidean space on a cone described by the intersection of an arbitrary number of halfspaces;
we want the closest point in the polyhedral cone.

By "formula" I mean a closed form; an equation or set of equations (not a program, algorithm, or optimization).
A set of formulae, the choice of which is conditional, is OK as long as size of the set is not factorial (prohibitively large).

This problem has many practical and theoretical applications. Its solution is certainly worth a Ph.D. thesis in any Math or Engineering Department.

You are welcome and encouraged to write your thoughts about this problem here.


Projection on isotone projection cones

Together with my coauthor A. B. Németh we recently discovered a very simple algorithm in LaTeX: \mathbb{R}^n for projecting onto a special class of cones: the isotone projection cones.

An isotone projection cone is a generating pointed closed convex cone LaTeX: \mathcal{K} in a Hilbert space LaTeX: \mathbb{H} for which projection LaTeX: P_{\mathcal{K}} onto the cone is isotone; that is, monotone with respect to the order LaTeX: \preceq induced by the cone: LaTeX: \forall\,x\,,\,y\in\mathbb{H}

LaTeX: x\preceq y\Rightarrow P_{\mathcal{K}}(x)\preceq P_{\mathcal{K}}(y),

or equivalently

LaTeX: y-x\in\mathcal{K}\Rightarrow P_{\mathcal{K}}(y)-P_{\mathcal{K}}(x)\in\mathcal{K}.

From now on, suppose that we are in LaTeX: \mathbb{R}^n. Here the isotone projection cones are polyhedral cones generated by LaTeX: n linearly independent vectors (i.e., cones that define a lattice structure; called latticial or simplicial cones) such that the generators of the polar cone form pairwise nonacute angles. For example, the monotone nonnegative cone (Example in CO&EDG) is an isotone projection cone. The monotone nonnegative cone is defined by

LaTeX: \left\{(x_1,...,x_n)^\top\in\mathbb{R}^n\,\mid\,x_1\geq...\geq x_n\geq 0\right\}.

Projecting onto monotone nonnegative cones (see Section in CO&EDG) is important for the problem of map-making from relative distance information (see Section 5.13 in CO&EDG); e.g., stellar cartography. If LaTeX: w_1,...,w_n>0 an extension of the monotone nonnegative cone is defined by

LaTeX: \left\{(x_1,...,x_n)^\top\in\mathbb{R}^n\,\mid\,\frac{x_1}{\sqrt{w_1}}\geq...\geq\frac{x_n}{\sqrt{w_n}}\geq 0\right\}.

This cone is also an isotone projection cone (see here). Projection onto this cone is equivalent to the isotonic regression problem with nonnegative variables given by the weight vector LaTeX: w=(w_1,...,w_n)^\top and a total order. The isotone projection cones were introduced by George Isac and A. B. Németh in the mid-1980's and then applied by George Isac, A. B. Németh, and S. Z. Németh to complementarity problems (see 1, 2, 3 and 4). The algorithm below opens the way for efficient implementations of the iterative methods in 1, 2 and 4. For simplicity we shall call a matrix LaTeX: \mathbf E , whose columns are generators of an isotone projection cone, an isotone projection cone generator matrix. Recall that a Stieltjes matrix is a real symmetric positive definite matrix with nonpositive off-diagonal entries. A matrix LaTeX: \mathbf E is an isotone projection cone generator matrix if and only if it is of the form

LaTeX: \mathbf E=\mathbf O\mathbf T^{-\frac12},

where LaTeX: \mathbf O is an orthogonal matrix and LaTeX: \mathbf T is Stieltjes matrix.

The algorithm is a finite one that stops in at most LaTeX: \,n steps and finds the exact projection point. Suppose that we want to project onto a latticial cone, and for each point in Euclidean space we know a proper face of the cone onto which that point projects. Then we could find the projection, recursively, by projecting onto linear subspaces of decreasing dimension. In case of isotone projection cones, the isotonicity property provides information required about such a proper face. The information is provided by geometrical structure of the polar cone. The theoretical background for the algorithm is Moreau's decomposition theorem.

If a polyhedral cone can be written as a union of isotone projection cones, reasonably small in number, then we could project a point onto the polyhedral cone by projecting onto the isotone projection cones and then taking the minimum distance of the given point from all these cones. Due to simplicity of the method for projecting onto an isotone projection cone, it is a challenging open question to find polyhedral cones that comprise a union of a small number of isotone projection cones that can be easily discerned.

The paper containing our method can be found here.

Scilab code can be downloaded here.

Matlab code for the algorithm:

% You are free to use, redistribute, and modifiy this code if you include,
% as a comment, the author of the original code
% (c) Sandor Zoltan Nemeth, 2009.
function p=piso(x,E)
if det(E)==0, error('The input cone must be an isotone projection cone'); end
for i=1:n
   for j=1:n
      if G(i,j)>0
         error('The input cone must be an isotone projection cone');
for k=0:n1
   if l-1>=1
      if highest<n
         for h=n:-1:highest+1
      for j=l-1:-1:1
         if high>=low
            for m=high:-1:low
         if lowest>1
            for w=lowest-1:-1:1
   if l==1


   for j=l:-1:1
      if zz>0, I1(j)=[];
   if cont==0, p=x; return
   elseif ll==0, p=zeros(n,1); return
   if ll==l, cont=0;
   else cont=1;

LaTeX: - Sándor Zoltán Németh

Fast projection on monotone nonnegative cone

Demo requires CVX to produce estimate of error by solving the projection problem in polynomial time.

CVX is not required for piso4() which performs projection by Sándor's method about 100 times faster.

%demo: Sandor's projection on monotone nonnegative cone
%-Jon Dattorro, August 2, 2009
clear all
n = 100;
a = randn(n,1);
t = piso4(a);
if n < 500
      variable s(n,1);
      variable b(n,1);
      for i=1:n
         s(i) == sum(b(i:n));
      b >= 0;
   err = sum(abs(s - t))


This code is optimized for memory efficiency. Otherwise, for LaTeX: n=millions of variables, required memory would be measured in terabytes.

% -Sandor Zoltan Nemeth, with code optimization -Jon Dattorro August 2, 2009.
% Project x on monotone nonnegative cone of dimension n.
function p = piso4(x)
n = length(x);
I = 1:n;
for k=1:n
   L = length(I);
   zeroidx = 1:n;
   zeroidx(I) = [];
   V1 = getpinv2(zeroidx,n);
   I(find(V1(I,:)*x<0)) = [];
   newL = length(I);
   if ~newL
      p = sparse(n,1);
      t = V1*x;
      s = sum(t);
      zs = sum(t(zeroidx));
      for i=1:n
         x(i) = s - zs;
         s = s - t(i);
         first = find(zeroidx>=i,1);
         zs = zs - sum(t(zeroidx(1:first-1)));
         zeroidx(1:first-1) = [];
   if newL==L 
      p = x;


%Assumes Platonic upper triangular ones (generator matrix) with some columns missing.
%Quickly finds pseudoinverse of monotone nonnegative cone generator matrix
%June 14, 2009  -Jon Dattorro

function Y = getpinv2(zeroidx,n);
Y = spdiags([ones(n,1), -ones(n,1)], [0 1], sparse(n,n));
Y(zeroidx,:) = 0;
%find dangling -1
count = 0;
for i=1:n
   if ~Y(i,i)
      count = count + 1;
      if i-1 > 0
         Y(i-1,i) = 0;
      if count
         if i-count-1 > 0
            Y(i-count-1,i-count:i) = -1/(count+1);
         Y(i,i-count:i) = 1/(count+1);
      count = 0;

Projection on monotone nonnegative cone using Pool Adjacent Violators type algorithms

The monotone cone is defined by

LaTeX: \left\{(x_1,...,x_n)^\top\in\mathbb{R}^n\,\mid\,x_1\geq...\geq x_n\right\}.

By using Pool Adjacent Violators (PAV) type algorithms, it is easy to project onto the monotone cone. However, the PAV type algorithms cannot be easily adapted to project onto the monotone nonnegative cone too. By using Moreau's decomposition theorem, A. B. Németh and S. Z. Németh have recently showed that one can project onto the monotone nonnegative cone by first projecting onto the monotone cone (using PAV type algorithms) and then replacing each negative coordinate of the projection by zero while keeping the nonnegative ones unchanged. The proof of this result can be found here.

Projection on simplicial cones

The following method was developed together with my coauthors A. Ekárt and A. B. Németh. Thanks are due to J. Dattorro for his continuous support and encouragement.

Let LaTeX: \mathcal{K} be a simplicial cone generated by the linearly independent vectors LaTeX: \,e^i in LaTeX: \mathbb{R}^n and LaTeX: \mathcal{K}^\circ be its polar generated by the linearly independent vectors LaTeX: \,u^j. We can suppose that LaTeX: (e^i)^\top u^j = -\delta_{ij}, where LaTeX: \,\delta_{ij} is the Kronecker symbol.


Let LaTeX: x\in\mathbb{R}^n. For each subset of indices LaTeX: I\subset\{1,...,n\}, LaTeX: \,x can be represented in the form

LaTeX: x=\sum_{i\in I}\alpha_i e^i+\sum_{j\in I^c}\beta_j u^j\quad [1]

with LaTeX: \,I^c the complement of LaTeX: I with respect to LaTeX: \,\{1,...,n\}, and with LaTeX: \,\alpha_i and LaTeX: \,\beta_j real numbers. Among the subsets LaTeX: \,I of indices, there exists exactly one (the cases LaTeX: I=\emptyset and LaTeX: \,I=\{1,...,n\} are not excluded) with the property that for the coefficients in LaTeX: \,[1] one has LaTeX: \,\beta_j>0, LaTeX: j\in I^c and LaTeX: \,\alpha_i\geq 0, LaTeX: \,i\in I. For this representation it holds that

LaTeX: P_{\mathcal{K}}(x)=\sum_{i\in I} \alpha_i e^i,\quad \alpha_i\geq 0


LaTeX: P_{\mathcal{K}^\circ}(x)=\sum_{j\in I^c}\beta_j u^j,\quad \beta^j>0.

This theorem suggests the following algorithm for finding the projection LaTeX: P_{\mathcal{K}}(x):

Step 1. For the subset LaTeX: I\subset\{1,...,n\} we solve the following linear system in LaTeX: \alpha^i

LaTeX: x^\top e^\ell =\sum_{i\in I}\alpha_i (e^i)^\top e^\ell,\;\ell\in I.\quad [2]

Step 2. Then, we select from the family of all subsets in LaTeX: \{1,...,n\} the subfamily LaTeX: \,\Delta of subsets LaTeX: \,I for which the system possesses non-negative solutions.

Step 3. For each LaTeX: I\in \Delta we solve the linear system in LaTeX: \,\beta^j

LaTeX: x^\top u^k = \sum_{j\in I^c}\beta_j (u^j)^\top u^k,\;k\in I^c.\quad [3]

By the above Theorem among these systems there exists exactly one with non-negative solutions. By the same theorem, for corresponding LaTeX: \,I and for the solution of the system LaTeX: \,[2], we must have

LaTeX: P_{\mathcal{K}}(x)=\sum_{i\in I}\alpha_i e^i.

This algorithm requires that we solve LaTeX: \,2^n linear systems of at most LaTeX: \,n equations in Step 1 LaTeX: \,[2] and another LaTeX: \,2^{|\Delta|} systems in Step 2 LaTeX: \,[3]. (Observe that all these systems are given by Gram matrices, hence they have unique solutions.) The above presented method is inefficient but it suggests the following heuristic algorithm:

Substitute LaTeX: \,u^j with LaTeX: \,e^j if LaTeX: \,\beta_j<0 and substitute LaTeX: \,e^i with LaTeX: \,u^i if LaTeX: \,\alpha_i<0 and solve the systems LaTeX: \,[2] and LaTeX: \,[3] for the new configuration of indices LaTeX: \,I. We shall call this step an iteration of the heuristic algorithm.

Then, repeat the procedure for the new configuration of LaTeX: \,I and so on, until we obtain a representation of LaTeX: \,x in LaTeX: \,[1] with all the coefficients non-negative.

The extensive experiments indicate that the method furnishes the exact solution in more then LaTeX: \,99.7 percent of the cases. The average number of steps is LaTeX: \,5.67 (we have not found any examples which required more than LaTeX: \,13 steps) and the relative number of steps with respect to the dimension decreases dramatically. Roughly speaking, for high enough dimensions the absolute number of steps is independent of the dimension. The percentage of problems where the algorithm was aborted due to going in a loop exponentially decreases as the size increases. The percentage of loops above sizes of LaTeX: \,30 is less than LaTeX: \,0.001\%. Overall, loops were observed in less than LaTeX: \,0.3\% of the experiments, so the heuristic algorithm was successful more than LaTeX: 99.7\% of the time. Most of the loops were observed for problems of very small size (LaTeX: \,2, LaTeX: \,3, LaTeX: \,4, and LaTeX: \,5). We recently observed that most of the loops for problems of very small size are false loops generated by computational errors. Moreover, we proved that for size 2 there are no loops at all. Thus, all loops for size 2 are false loops. The conclusions are based on millions of numerical experiments. More experiments are needed to see if there are any true loops or all our loops are false ones. It is also open question whether similar proofs to the 2-dimensional one can be obtained for higher dimensions or not.

You can find full details about this method here.

Scilab code can be downloaded here.

LaTeX: - Sándor Zoltán Németh

external links

More about projection on cones (and convex sets, more generally) can be found here (chapter E):

More about polyhedral cones can be found in chapter 2.

Personal tools