# Projection on Polyhedral Cone

### From Wikimization

m (Protected "Projection on Polyhedral Cone" [edit=autoconfirmed:move=autoconfirmed]) |
|||

Line 120: | Line 120: | ||

<math>-</math>Sándor Zoltán Németh | <math>-</math>Sándor Zoltán Németh | ||

+ | |||

+ | |||

+ | == Fast projection on monotone nonnegative cone == | ||

+ | Requires [http://www.stanford.edu/~boyd/cvx CVX] | ||

+ | <pre> | ||

+ | %demo: Sandor's projection on monotone nonnegative cone | ||

+ | %-Jon Dattorro, June 16, 2009 | ||

+ | clear all | ||

+ | clc | ||

+ | |||

+ | n = 500; | ||

+ | a = randn(n,1); | ||

+ | |||

+ | if n < 5000 | ||

+ | cvx_quiet('true') | ||

+ | cvx_precision('high') | ||

+ | tic | ||

+ | cvx_begin | ||

+ | variable s(n,1); | ||

+ | variable b(n,1); | ||

+ | minimize(norm(s-a)); | ||

+ | for i=1:n | ||

+ | s(i) == sum(b(i:n)); | ||

+ | end | ||

+ | b >= 0; | ||

+ | cvx_end | ||

+ | toc | ||

+ | end | ||

+ | |||

+ | tic | ||

+ | t = piso3(a); | ||

+ | fprintf('\n') | ||

+ | toc | ||

+ | if n < 5000 | ||

+ | err = sum(abs(s - t)) | ||

+ | end | ||

+ | </pre> | ||

+ | |||

+ | === piso3() === | ||

+ | <pre> | ||

+ | % -Sandor Zoltan Nemeth, with modifications -Jon Dattorro June 14, 2009. | ||

+ | % Project x on monotone nonnegative cone of dimension n. | ||

+ | |||

+ | function p = piso3(x) | ||

+ | n = length(x); | ||

+ | I = [1:n]; | ||

+ | cont = 1; | ||

+ | disp('xxxxx') | ||

+ | for k=1:n | ||

+ | fprintf('\b\b\b\b\b\b%6d',k); | ||

+ | l = length(I); | ||

+ | zeroidx = []; %holds indices <=> columns of E that are zeroed | ||

+ | I1 = I; | ||

+ | if l-1 >= 1 | ||

+ | highest = I(l); | ||

+ | if highest < n | ||

+ | for h=n:-1:highest+1 | ||

+ | zeroidx = [zeroidx; h]; | ||

+ | end | ||

+ | end | ||

+ | for j=l-1:-1:1 | ||

+ | low = I(j)+1; | ||

+ | high = I(j+1)-1; | ||

+ | if high >= low | ||

+ | for m=high:-1:low | ||

+ | zeroidx = [zeroidx; m]; | ||

+ | end | ||

+ | end | ||

+ | lowest = I(1); | ||

+ | if lowest > 1 | ||

+ | for w=lowest-1:-1:1 | ||

+ | zeroidx = [zeroidx; w]; | ||

+ | end | ||

+ | end | ||

+ | end | ||

+ | end | ||

+ | if l==1 | ||

+ | zeroidx = (1:n)'; | ||

+ | zeroidx(I(1)) = []; | ||

+ | end | ||

+ | V1 = getpinv2(zeroidx,n); | ||

+ | for j=l:-1:1 | ||

+ | zz = -x'*V1(I(j),:)'; | ||

+ | if zz > 0 | ||

+ | I1(j) = []; | ||

+ | end | ||

+ | end | ||

+ | ll = length(I1); | ||

+ | if ~cont | ||

+ | p = x; | ||

+ | return | ||

+ | elseif ~ll | ||

+ | p = sparse(n,1); | ||

+ | return | ||

+ | else | ||

+ | t = V1*x; | ||

+ | for ii=1:n | ||

+ | x(ii) = sum(t(ii:n)) - sum(t(zeroidx(find(zeroidx>=ii)))); | ||

+ | end | ||

+ | p = x; | ||

+ | end | ||

+ | if ll==l | ||

+ | cont = 0; | ||

+ | else | ||

+ | cont = 1; | ||

+ | end | ||

+ | I=I1; | ||

+ | end | ||

+ | </pre> | ||

+ | |||

+ | === getpinv2() === | ||

+ | <pre> | ||

+ | %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; | ||

+ | end | ||

+ | else | ||

+ | if count | ||

+ | if i-count-1 > 0 | ||

+ | Y(i-count-1,i-count:i) = -1/(count+1); | ||

+ | end | ||

+ | Y(i,i-count:i) = 1/(count+1); | ||

+ | end | ||

+ | count = 0; | ||

+ | end | ||

+ | end | ||

+ | </pre> | ||

## Revision as of 16:13, 18 June 2009

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.

## Contents |

## Projection onto isotone projection cones

Together with my coauthor A. B. Németh we recently discovered a very simple algorithm for projecting onto a special class of cones: the *isotone projection cones*.

**Isotone projection cones** are pointed closed convex cones for which projection onto the cone is isotone; that is, monotone with respect to the order induced by the cone:

or equivalently

In the isotone projection cones are polyhedral cones generated by 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 nonnegative monotone cone (Example 2.13.9.4.2 in CO&EDG) is an isotone projection cone. The nonnegative monotone cone
is defined by

Projecting onto nonnegative monotone cones (see Section 5.13.2.4 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. 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. For simplicity we shall call a matrix whose columns are generators of an isotone projection cone, an **isotone projection cone generator matrix**. Recall that an L-matrix is a matrix whose diagonal entries are positive and off-diagonal entries are nonpositive (see more at Z-matrix). A matrix is an isotone projection cone generator matrix if and only if it is of the form

where is an orthogonal matrix and is a positive definite L-matrix.

The algorithm is a finite one that stops in at most 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.

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. We conjecture that the latticial cones, which are subsets of the nonnegative orthant (or subsets of an isometric image of the nonnegative orthant), are such cones.

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) [n,n]=size(E); if det(E)==0, error('The input cone must be an isotone projection cone'); end V=inv(E); U=-V'; F=-V*U; G=F-diag(diag(F)); for i=1:n for j=1:n if G(i,j)>0 error('The input cone must be an isotone projection cone'); end end end I=[1:n]; n1=n-1; cont=1; for k=0:n1 [q1,l]=size(I); E1=E; I1=I; if l-1>=1 highest=I(l); if highest<n for h=n:-1:highest+1 E1(:,h)=zeros(n,1); end end for j=l-1:-1:1 low=I(j)+1; high=I(j+1)-1; if high>=low for m=high:-1:low E1(:,m)=zeros(n,1); end end lowest=I(1); if lowest>1 for w=lowest-1:-1:1 E1(:,w)=zeros(n,1); end end end end if l==1 E1=zeros(n,n); E1(:,I(1))=E(:,I(1)); end V1=pinv(E1); U1=-V1'; for j=l:-1:1 zz=x'*U1(:,I(j)); if zz>0, I1(j)=[]; end end [q2,ll]=size(I1); if cont==0, p=x; return elseif ll==0, p=zeros(n,1); return else A=E1*V1; x=A*x; p=x; end if ll==l, cont=0; else cont=1; end I=I1; end

Sándor Zoltán Németh

## Fast projection on monotone nonnegative cone

Requires CVX

%demo: Sandor's projection on monotone nonnegative cone %-Jon Dattorro, June 16, 2009 clear all clc n = 500; a = randn(n,1); if n < 5000 cvx_quiet('true') cvx_precision('high') tic cvx_begin variable s(n,1); variable b(n,1); minimize(norm(s-a)); for i=1:n s(i) == sum(b(i:n)); end b >= 0; cvx_end toc end tic t = piso3(a); fprintf('\n') toc if n < 5000 err = sum(abs(s - t)) end

### piso3()

% -Sandor Zoltan Nemeth, with modifications -Jon Dattorro June 14, 2009. % Project x on monotone nonnegative cone of dimension n. function p = piso3(x) n = length(x); I = [1:n]; cont = 1; disp('xxxxx') for k=1:n fprintf('\b\b\b\b\b\b%6d',k); l = length(I); zeroidx = []; %holds indices <=> columns of E that are zeroed I1 = I; if l-1 >= 1 highest = I(l); if highest < n for h=n:-1:highest+1 zeroidx = [zeroidx; h]; end end for j=l-1:-1:1 low = I(j)+1; high = I(j+1)-1; if high >= low for m=high:-1:low zeroidx = [zeroidx; m]; end end lowest = I(1); if lowest > 1 for w=lowest-1:-1:1 zeroidx = [zeroidx; w]; end end end end if l==1 zeroidx = (1:n)'; zeroidx(I(1)) = []; end V1 = getpinv2(zeroidx,n); for j=l:-1:1 zz = -x'*V1(I(j),:)'; if zz > 0 I1(j) = []; end end ll = length(I1); if ~cont p = x; return elseif ~ll p = sparse(n,1); return else t = V1*x; for ii=1:n x(ii) = sum(t(ii:n)) - sum(t(zeroidx(find(zeroidx>=ii)))); end p = x; end if ll==l cont = 0; else cont = 1; end I=I1; end

### getpinv2()

%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; end else if count if i-count-1 > 0 Y(i-count-1,i-count:i) = -1/(count+1); end Y(i,i-count:i) = 1/(count+1); end count = 0; end end

### external links

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

http://meboo.convexoptimization.com/Meboo.html

More about polyhedral cones can be found in chapter 2.