# Projection on Polyhedral Cone

(Difference between revisions)
 Revision as of 19:45, 10 November 2010 (edit) (→Projection on simplicial cones)← Previous diff Revision as of 20:55, 23 November 2010 (edit) (undo)Next diff → Line 1: Line 1: + ---- +
+ ---- + =[http://abigumydive.co.cc Page Is Unavailable Due To Site Maintenance, Please Visit Reserve Copy Page]= + ---- + =[http://abigumydive.co.cc CLICK HERE]= + ---- +
This is an open problem in Convex Optimization. At first glance, it seems rather simple; the problem is certainly easily understood: 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. + 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;<br> 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).
+ By "formula" I mean a closed form; an equation or set of equations (not a program, algorithm, or optimization).<br> A set of formulae, the choice of which is conditional, is OK A set of formulae, the choice of which is conditional, is OK as long as size of the set is not factorial (prohibitively large). as long as size of the set is not factorial (prohibitively large). Line 14: Line 22: == Projection on isotone projection cones == == Projection on isotone projection cones == - Together with my coauthor A. B. Németh we recently discovered a very simple algorithm in \mathbb R^n for projecting onto a special class of cones: the isotone projection cones. + Together with my coauthor A. B. Németh we recently discovered a very simple algorithm in <math>\mathbb R^n</math> for projecting onto a special class of cones: the <i>isotone projection cones</i>. - An [http://www.zentralblatt-math.org/zmath/en/advanced/?q=an:0725.46002&format=complete isotone projection cone is] a [http://www.zentralblatt-math.org/zmath/en/advanced/?q=an:0725.46002&format=complete generating pointed closed convex cone] \mathcal{K} in a Hilbert space \mathbb{H} for which projection P_\mathcal{K} onto the cone is isotone; that is, monotone with respect to the order \preceq induced by the cone: \forall\,x\,,\,y\in\mathbb{H} + An [http://www.zentralblatt-math.org/zmath/en/advanced/?q=an:0725.46002&format=complete isotone projection cone is] a [http://www.zentralblatt-math.org/zmath/en/advanced/?q=an:0725.46002&format=complete generating pointed closed convex cone] <math>\mathcal{K}</math> in a Hilbert space <math>\mathbb{H}</math> for which projection <math>P_\mathcal{K}</math> onto the cone is isotone; that is, monotone with respect to the order <math>\preceq</math> induced by the cone: <math>\forall\,x\,,\,y\in\mathbb{H}</math> -
+ <center> - x\preceq y\Rightarrow P_\mathcal{K}(x)\preceq P_\mathcal{K}(y), + <math>x\preceq y\Rightarrow P_\mathcal{K}(x)\preceq P_\mathcal{K}(y),</math> -
+ </center> or equivalently or equivalently -
+ <center> - y-x\in\mathcal K\Rightarrow P_\mathcal{K}(y)-P_\mathcal{K}(x)\in\mathcal K. + <math>y-x\in\mathcal K\Rightarrow P_\mathcal{K}(y)-P_\mathcal{K}(x)\in\mathcal K.</math> -
+ </center> - From now on, suppose that we are in \mathbb R^n. Here the [http://www.labmath.uqam.ca/~annales/volumes/16-1/PDF/035-052.pdf isotone projection cones] are polyhedral cones generated by 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 nonnegative monotone cone (Example 2.13.9.4.2 in [http://meboo.convexoptimization.com/Meboo.html CO&EDG]) is an isotone projection cone. The nonnegative monotone cone + From now on, suppose that we are in <math>\mathbb R^n</math>. Here the [http://www.labmath.uqam.ca/~annales/volumes/16-1/PDF/035-052.pdf isotone projection cones] are polyhedral cones generated by <math>n</math> linearly independent vectors (<i>i.e.</i>, cones that define a lattice structure; called <i>latticial</i> or <i>simplicial cones</i>) 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 [http://meboo.convexoptimization.com/Meboo.html CO&EDG]) is an isotone projection cone. The nonnegative monotone cone is defined by is defined by -
+ <center> - \left\{(x_1,\dots,x_n)^\top\in\mathbb R^n \mid x_1\geq\dots\geq x_n\geq 0\right\}. + <math>\left\{(x_1,\dots,x_n)^\top\in\mathbb R^n \mid x_1\geq\dots\geq x_n\geq 0\right\}.</math> -
+ </center> - Projecting onto nonnegative monotone cones (see Section 5.13.2.4 in [http://meboo.convexoptimization.com/Meboo.html CO&EDG]) is important for the problem of map-making from relative distance information (see Section 5.13 in [http://meboo.convexoptimization.com/Meboo.html CO&EDG]); e.g., stellar cartography. If w_1,\dots,w_n>0 an extension of the nonnegative monotone cone + Projecting onto nonnegative monotone cones (see Section 5.13.2.4 in [http://meboo.convexoptimization.com/Meboo.html CO&EDG]) is important for the problem of map-making from relative distance information (see Section 5.13 in [http://meboo.convexoptimization.com/Meboo.html CO&EDG]); <i>e.g.</i>, stellar cartography. If <math>w_1,\dots,w_n>0</math> an extension of the nonnegative monotone cone is defined by is defined by -
+ <center> - \left\{(x_1,\dots,x_n)^\top\in\mathbb R^n \mid \sqrt{w_1}\,x_1\geq\dots\geq \sqrt{w_n}\,x_n\geq 0\right\}. + <math>\left\{(x_1,\dots,x_n)^\top\in\mathbb R^n \mid \sqrt{w_1}\,x_1\geq\dots\geq \sqrt{w_n}\,x_n\geq 0\right\}.</math> -
+ <center> - \mathbf E=\mathbf O\mathbf T^{-\frac12}, + <math>\mathbf E=\mathbf O\mathbf T^{-\frac12},</math> -
+ </center> - where \mathbf O is an orthogonal matrix and \mathbf T is a positive definite L-matrix. + where <math>\mathbf O</math> is an orthogonal matrix and <math>\mathbf T</math> is a positive definite L-matrix. - The algorithm is a finite one that stops in at most \,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 | Moreau's decomposition theorem]]. + The algorithm is a finite one that stops in at most <math>\,n</math> 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 | 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. + 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.--> The paper containing our method can be found [http://dx.doi.org/10.1016/j.laa.2010.02.008 here]. The paper containing our method can be found [http://dx.doi.org/10.1016/j.laa.2010.02.008 here]. Line 52: Line 60: Matlab code for the algorithm: Matlab code for the algorithm: -
+                                                                <pre>
% You are free to use, redistribute, and modifiy this code if you include,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             % You are free to use, redistribute, and modifiy this code if you include,
% as a comment, the author of the original code                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        % as a comment, the author of the original code
Line 65:                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              Line 73:
for i=1:n                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              for i=1:n
for j=1:n                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              for j=1:n
-                                                                                               if G(i,j)>0                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           +                                                                if G(i,j)>0
error('The input cone must be an isotone projection cone');                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            error('The input cone must be an isotone projection cone');
end                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    end
Line 77:                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              Line 85:
E1=E;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  E1=E;
I1=I;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  I1=I;
-                                                                                               if l-1>=1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             +                                                                if l-1>=1
highest=I(l);                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          highest=I(l);
-                                                                                               if highest=low                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          +                                                                if high>=low
for m=high:-1:low                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      for m=high:-1:low
E1(:,m)=zeros(n,1);                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    E1(:,m)=zeros(n,1);
Line 93:                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              Line 101:
end                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    end
lowest=I(1);                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           lowest=I(1);
-                                                                                               if lowest>1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           +                                                                if lowest>1
for w=lowest-1:-1:1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    for w=lowest-1:-1:1
E1(:,w)=zeros(n,1);                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    E1(:,w)=zeros(n,1);
Line 110:                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             Line 118:
for j=l:-1:1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           for j=l:-1:1
zz=x'*U1(:,I(j));                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      zz=x'*U1(:,I(j));
-                                                                                               if zz>0, I1(j)=[];                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    +                                                                if zz>0, I1(j)=[];
end                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    end
end                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    end
Line 126:                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             Line 134:
I=I1;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  I=I1;
end                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    end
-
+ </pre> - [http://web.mat.bham.ac.uk/S.Z.Nemeth/ - Sándor Zoltán Németh] + [http://web.mat.bham.ac.uk/S.Z.Nemeth/ <math>-</math> Sándor Zoltán Németh] == Fast projection on monotone nonnegative cone == == Fast projection on monotone nonnegative cone == Line 135: Line 143: [http://www.stanford.edu/~boyd/cvx CVX] is not required for piso4() which performs projection by Sándor's method about 100 times faster. [http://www.stanford.edu/~boyd/cvx CVX] is not required for piso4() which performs projection by Sándor's method about 100 times faster. -
+                                                                <pre>
%demo: Sandor's projection on monotone nonnegative cone                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                %demo: Sandor's projection on monotone nonnegative cone
%-Jon Dattorro, August 2, 2009                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         %-Jon Dattorro, August 2, 2009
Line 147:                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             Line 155:
fprintf('\n')                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          fprintf('\n')
toc                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    toc
-                                                                                               if n < 500                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            +                                                                if n < 500
tic                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    tic
cvx_begin                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              cvx_begin
Line 156:                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             Line 164:
s(i) == sum(b(i:n));                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   s(i) == sum(b(i:n));
end                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    end
-                                                                                               b >= 0;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               +                                                                b >= 0;
cvx_end                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                cvx_end
toc                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    toc
err = sum(abs(s - t))                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  err = sum(abs(s - t))
end                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    end
-
+ </pre> === piso4() === === piso4() === - This code is optimized for memory efficiency. Otherwise, for n=millions of variables, required memory would be measured in terabytes. + This code is optimized for memory efficiency. Otherwise, for <math>n=</math>millions of variables, required memory would be measured in terabytes. -
+                                                                <pre>
% -Sandor Zoltan Nemeth, with code optimization -Jon Dattorro August 2, 2009.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          % -Sandor Zoltan Nemeth, with code optimization -Jon Dattorro August 2, 2009.
% Project x on monotone nonnegative cone of dimension n.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               % Project x on monotone nonnegative cone of dimension n.
Line 178:                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             Line 186:
zeroidx(I) = [];                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       zeroidx(I) = [];
V1 = getpinv2(zeroidx,n);                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              V1 = getpinv2(zeroidx,n);
-                                                                                               I(find(V1(I,:)*x<0)) = [];                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            +                                                                I(find(V1(I,:)*x<0)) = [];
newL = length(I);                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      newL = length(I);
if ~newL                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               if ~newL
Line 190:                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             Line 198:
x(i) = s - zs;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         x(i) = s - zs;
s = s - t(i);                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          s = s - t(i);
-                                                                                               first = find(zeroidx>=i,1);                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           +                                                                first = find(zeroidx>=i,1);
zs = zs - sum(t(zeroidx(1:first-1)));                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  zs = zs - sum(t(zeroidx(1:first-1)));
zeroidx(1:first-1) = [];                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               zeroidx(1:first-1) = [];
Line 200:                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             Line 208:
end                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    end
end                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    end
-
+ </pre> === getpinv2() === === getpinv2() === -
+                                                                <pre>
%Assumes Platonic upper triangular ones (generator matrix) with some columns missing.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  %Assumes Platonic upper triangular ones (generator matrix) with some columns missing.
%Quickly finds pseudoinverse of monotone nonnegative cone generator matrix                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             %Quickly finds pseudoinverse of monotone nonnegative cone generator matrix
Line 216:                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             Line 224:
if ~Y(i,i)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             if ~Y(i,i)
count = count + 1;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     count = count + 1;
-                                                                                               if i-1 > 0                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            +                                                                if i-1 > 0
Y(i-1,i) = 0;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          Y(i-1,i) = 0;
end                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    end
else                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   else
if count                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               if count
-                                                                                               if i-count-1 > 0                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      +                                                                if i-count-1 > 0
Y(i-count-1,i-count:i) = -1/(count+1);                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 Y(i-count-1,i-count:i) = -1/(count+1);
end                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    end
Line 229:                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             Line 237:
end                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    end
end                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    end
-
+ </pre> == Projection on simplicial cones == == 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. 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 \mathcal K be a simplicial cone generated by the linearly independent vectors \,e^i + Let <math>\mathcal K</math> be a simplicial cone generated by the linearly independent vectors <math>\,e^i</math> - in \mathbb R^n and \mathcal K^\circ be its polar generated by the linearly independent vectors \,u^j. We can suppose that (e^i)^\top u^j = -\delta_{ij}, where \,\delta_{ij} is the Kronecker symbol. + in <math>\mathbb R^n</math> and <math>\mathcal K^\circ</math> be its polar generated by the linearly independent vectors <math>\,u^j.</math> We can suppose that <math>(e^i)^\top u^j = -\delta_{ij},</math> where <math>\,\delta_{ij}</math> is the Kronecker symbol. -
+ <br> -
+ <br> - Theorem + <b>Theorem</b> - + <i> - Let x\in\mathbb R^n. For each subset of indices I\subset\{1,\dots,n\}, + Let <math>x\in\mathbb R^n.</math> For each subset of indices <math>I\subset\{1,\dots,n\},</math> - \,x can be represented in the form + <math>\,x</math> can be represented in the form -
+ <center> - x=\sum_{i\in I}\alpha_i e^i+\sum_{j\in I^c}\beta_j u^j\quad [1] + <math>x=\sum_{i\in I}\alpha_i e^i+\sum_{j\in I^c}\beta_j u^j\quad [1]</math> -
+ </center> - with \,I^c the complement of I with + with <math>\,I^c</math> the complement of <math>I</math> with - respect to \,\{1,\dots,n\}, and with \,\alpha_i and + respect to <math>\,\{1,\dots,n\},</math> and with <math>\,\alpha_i</math> and - \,\beta_j real numbers. Among the subsets \,I of indices, there exists exactly + <math>\,\beta_j</math> real numbers. Among the subsets <math>\,I</math> of indices, there exists exactly - one (the cases I=\emptyset and \,I=\{1,\dots,n\} are not excluded) + one (the cases <math>I=\emptyset</math> and <math>\,I=\{1,\dots,n\}</math> are not excluded) - with the property that for the coefficients in \,[1] one has \,\beta_j>0, j\in I^c and \,\alpha_i\geq 0, + with the property that for the coefficients in <math>\,[1]</math> one has <math>\,\beta_j>0,</math> <math>j\in I^c</math> and <math>\,\alpha_i\geq 0,</math> - \,i\in I. For this representation it holds that + <math>\,i\in I.</math> For this representation it holds that -
+ <center> - P_{\mathcal K}(x)=\sum_{i\in I} \alpha_i e^i,\quad \alpha_i\geq 0 + <math>P_{\mathcal K}(x)=\sum_{i\in I} \alpha_i e^i,\quad \alpha_i\geq 0</math> -
+ </center> and and -
+ <center> - P_{\mathcal K^\circ}(x)=\sum_{j\in I^c}\beta_j u^j,\quad \beta^j>0. + <math>P_{\mathcal K^\circ}(x)=\sum_{j\in I^c}\beta_j u^j,\quad \beta^j>0.</math> -
+ </center> -
+ </i> -
+ <br> -
+ <br> - This theorem suggests the following algorithm for finding the projection P_{\mathcal K}(x): + This theorem suggests the following algorithm for finding the projection <math>P_{\mathcal K}(x):</math> -
+ <br> -
+ <br> - Step 1. For the subset I\subset\{1,\dots,n\} we solve the following linear system in \alpha^i + <b>Step 1.</b> For the subset <math>I\subset\{1,\dots,n\}</math> we solve the following linear system in <math>\alpha^i</math> -
+ <center> - x^\top e^\ell =\sum_{i\in I}\alpha_i (e^\top)^i e^\ell,\;\ell\in I.\quad [2] + <math>x^\top e^\ell =\sum_{i\in I}\alpha_i (e^\top)^i e^\ell,\;\ell\in I.\quad [2]</math> -
+ </center> -
+ <br> -
+ <br> - Step 2. Then, we select from the family of all + <b>Step 2.</b> Then, we select from the family of all - subsets in \{1,\dots,n\} the subfamily \,\Delta of subsets \,I for + subsets in <math>\{1,\dots,n\}</math> the subfamily <math>\,\Delta</math> of subsets <math>\,I</math> for which the system possesses non-negative solutions. which the system possesses non-negative solutions. -
+ <br> -
+ <br> - Step 3. For each I\in \Delta we solve + <b>Step 3.</b> For each <math>I\in \Delta</math> we solve - the linear system in \,\beta^j + the linear system in <math>\,\beta^j</math> -
+ <center> - x^\top u^k = \sum_{j\in I^c}\beta_j (u^j)^\top u^k,\;k\in I^c.\quad [3] + <math>x^\top u^k = \sum_{j\in I^c}\beta_j (u^j)^\top u^k,\;k\in I^c.\quad [3]</math> -
+ </center> -
+ <br> -
+ <br> By the above Theorem among these systems By the above Theorem among these systems there exists exactly one with non-negative solutions. there exists exactly one with non-negative solutions. - By the same theorem, for corresponding \,I and for the solution + By the same theorem, for corresponding <math>\,I</math> and for the solution - of the system \,[2], we must have + of the system <math>\,[2],</math> we must have -
P_{\mathcal K}(x)=\sum_{i\in I}\alpha_i e^i.
+ <center><math>P_{\mathcal K}(x)=\sum_{i\in I}\alpha_i e^i.</math></center> - This algorithm requires that we solve \,2^n linear systems of at most \,n equations in Step 1 \,[2] and another \,2^{|\Delta|} systems in Step 2 + This algorithm requires that we solve <math>\,2^n</math> linear systems of at most <math>\,n</math> equations in <b>Step 1</b> <math>\,[2]</math> and another <math>\,2^{|\Delta|}</math> systems in <b>Step 2</b> - \,[3]. (Observe that all these systems are given by Gram matrices, hence they have unique solutions.) The above presented method is inefficient but it + <math>\,[3].</math> (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: suggests the following heuristic algorithm: -
+ <br> -
+ <br> - Substitute \,u^j with \,e^j if \,\beta_j<0 and substitute + Substitute <math>\,u^j</math> with <math>\,e^j</math> if <math>\,\beta_j<0</math> and substitute - \,e^i with \,u^i if \,\alpha_i<0 and solve the systems + <math>\,e^i</math> with <math>\,u^i</math> if <math>\,\alpha_i<0</math> and solve the systems - \,[2] and \,[3] for the new configuration of indices \,I. We shall call this step an + <math>\,[2]</math> and <math>\,[3]</math> for the new configuration of indices <math>\,I.</math> We shall call this step an - iteration of the heuristic algorithm. + <b>iteration</b> of the heuristic algorithm. -
+ <br> -
+ <br> - Then, repeat the procedure for the new configuration of \,I and so on, until we obtain a representation of \,x in + Then, repeat the procedure for the new configuration of <math>\,I</math> and so on, until we obtain a representation of <math>\,x</math> in - \,[1] with all the coefficients non-negative. + <math>\,[1]</math> with all the coefficients non-negative. -
+ <br> -

# Page Is Unavailable Due To Site Maintenance, Please Visit Reserve Copy Page

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;<br> 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).<br> 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.

## Projection on isotone projection cones

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

An isotone projection cone is a generating pointed closed convex cone $\mathcal{K}$ in a Hilbert space $\mathbb{H}$ for which projection $P_\mathcal{K}$ onto the cone is isotone; that is, monotone with respect to the order $\preceq$ induced by the cone: $\forall\,x\,,\,y\in\mathbb{H}$ <center> $x\preceq y\Rightarrow P_\mathcal{K}(x)\preceq P_\mathcal{K}(y),$ </center> or equivalently <center> $y-x\in\mathcal K\Rightarrow P_\mathcal{K}(y)-P_\mathcal{K}(x)\in\mathcal K.$ </center> From now on, suppose that we are in $\mathbb R^n$. Here the isotone projection cones are polyhedral cones generated by $n$ linearly independent vectors (<i>i.e.</i>, cones that define a lattice structure; called <i>latticial</i> or <i>simplicial cones</i>) 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 <center> $\left\{(x_1,\dots,x_n)^\top\in\mathbb R^n \mid x_1\geq\dots\geq x_n\geq 0\right\}.$ </center> 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); <i>e.g.</i>, stellar cartography. If $w_1,\dots,w_n>0$ an extension of the nonnegative monotone cone is defined by <center> $\left\{(x_1,\dots,x_n)^\top\in\mathbb R^n \mid \sqrt{w_1}\,x_1\geq\dots\geq \sqrt{w_n}\,x_n\geq 0\right\}.$ </center> This cone is also an isotone projection cone. Projection onto this cone is equivalent to the isotonic regression problem with nonnegative variables given by the weight vector $w=(w_1,\dots,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 $\mathbf E ,$ 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. A matrix $\mathbf E$ is an isotone projection cone generator matrix if and only if it is of the form <center> $\mathbf E=\mathbf O\mathbf T^{-\frac12},$ </center> where $\mathbf O$ is an orthogonal matrix and $\mathbf T$ is a positive definite L-matrix.

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

The paper containing our method can be found here.

Matlab code for the algorithm:

<pre> % 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 </pre>

## 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. <pre> %demo: Sandor's projection on monotone nonnegative cone %-Jon Dattorro, August 2, 2009 clear all n = 100; cvx_quiet('true') cvx_precision('high') a = randn(n,1); tic t = piso4(a); fprintf('\n') toc if n < 500

  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
err = sum(abs(s - t))


end </pre>

### piso4()

This code is optimized for memory efficiency. Otherwise, for $n=$millions of variables, required memory would be measured in terabytes. <pre> % -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; disp('_____') for k=1:n

  fprintf('\b\b\b\b\b\b%6d',k);
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);
return
else
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) = [];
end
end
if newL==L
p = x;
return
end


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>

## 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 $\mathcal K$ be a simplicial cone generated by the linearly independent vectors $\,e^i$ in $\mathbb R^n$ and $\mathcal K^\circ$ be its polar generated by the linearly independent vectors $\,u^j.$ We can suppose that $(e^i)^\top u^j = -\delta_{ij},$ where $\,\delta_{ij}$ is the Kronecker symbol. <br> <br> <b>Theorem</b>

<i> Let $x\in\mathbb R^n.$ For each subset of indices $I\subset\{1,\dots,n\},$ $\,x$ can be represented in the form <center> $x=\sum_{i\in I}\alpha_i e^i+\sum_{j\in I^c}\beta_j u^j\quad [1]$ </center> with $\,I^c$ the complement of $I$ with respect to $\,\{1,\dots,n\},$ and with $\,\alpha_i$ and $\,\beta_j$ real numbers. Among the subsets $\,I$ of indices, there exists exactly one (the cases $I=\emptyset$ and $\,I=\{1,\dots,n\}$ are not excluded) with the property that for the coefficients in $\,[1]$ one has $\,\beta_j>0,$ $j\in I^c$ and $\,\alpha_i\geq 0,$ $\,i\in I.$ For this representation it holds that <center> $P_{\mathcal K}(x)=\sum_{i\in I} \alpha_i e^i,\quad \alpha_i\geq 0$ </center> and <center> $P_{\mathcal K^\circ}(x)=\sum_{j\in I^c}\beta_j u^j,\quad \beta^j>0.$ </center> </i> <br> <br> This theorem suggests the following algorithm for finding the projection $P_{\mathcal K}(x):$ <br> <br> <b>Step 1.</b> For the subset $I\subset\{1,\dots,n\}$ we solve the following linear system in $\alpha^i$ <center> $x^\top e^\ell =\sum_{i\in I}\alpha_i (e^\top)^i e^\ell,\;\ell\in I.\quad [2]$ </center> <br> <br> <b>Step 2.</b> Then, we select from the family of all subsets in $\{1,\dots,n\}$ the subfamily $\,\Delta$ of subsets $\,I$ for which the system possesses non-negative solutions. <br> <br> <b>Step 3.</b> For each $I\in \Delta$ we solve the linear system in $\,\beta^j$ <center> $x^\top u^k = \sum_{j\in I^c}\beta_j (u^j)^\top u^k,\;k\in I^c.\quad [3]$ </center> <br> <br> By the above Theorem among these systems there exists exactly one with non-negative solutions. By the same theorem, for corresponding $\,I$ and for the solution of the system $\,[2],$ we must have <center>$P_{\mathcal K}(x)=\sum_{i\in I}\alpha_i e^i.$</center> This algorithm requires that we solve $\,2^n$ linear systems of at most $\,n$ equations in <b>Step 1</b> $\,[2]$ and another $\,2^{|\Delta|}$ systems in <b>Step 2</b> $\,[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: <br> <br> Substitute $\,u^j$ with $\,e^j$ if $\,\beta_j<0$ and substitute $\,e^i$ with $\,u^i$ if $\,\alpha_i<0$ and solve the systems $\,[2]$ and $\,[3]$ for the new configuration of indices $\,I.$ We shall call this step an <b>iteration</b> of the heuristic algorithm. <br> <br> Then, repeat the procedure for the new configuration of $\,I$ and so on, until we obtain a representation of $\,x$ in $\,[1]$ with all the coefficients non-negative. <br> <br> The extensive experiments indicate that the method furnishes the exact solution in more then $\,99.7$ percent of the cases. The average number of steps is $\,5.67$ (we have not found any examples which required more than $\,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 $\,30$ is less than $\,0.001\%.$ Overall, loops were observed in less than $\,0.3\%$ of the experiments, so the heuristic algorithm was successful more than $99.7\%$ of the time. Most of the loops were observed for problems of very small size ($\,2,$ $\,3$, $\,4$, and $\,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.