[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

vigra/quadprog.hxx
00001 /************************************************************************/
00002 /*                                                                      */
00003 /*                  Copyright 2008 by Ullrich Koethe                    */
00004 /*                                                                      */
00005 /*    This file is part of the VIGRA computer vision library.           */
00006 /*    The VIGRA Website is                                              */
00007 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
00008 /*    Please direct questions, bug reports, and contributions to        */
00009 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00010 /*        vigra@informatik.uni-hamburg.de                               */
00011 /*                                                                      */
00012 /*    Permission is hereby granted, free of charge, to any person       */
00013 /*    obtaining a copy of this software and associated documentation    */
00014 /*    files (the "Software"), to deal in the Software without           */
00015 /*    restriction, including without limitation the rights to use,      */
00016 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00017 /*    sell copies of the Software, and to permit persons to whom the    */
00018 /*    Software is furnished to do so, subject to the following          */
00019 /*    conditions:                                                       */
00020 /*                                                                      */
00021 /*    The above copyright notice and this permission notice shall be    */
00022 /*    included in all copies or substantial portions of the             */
00023 /*    Software.                                                         */
00024 /*                                                                      */
00025 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00026 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00027 /*    OF MERCHANTABILITY, FITNESS FOR activeSet PARTICULAR PURPOSE AND  */
00028 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00029 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00030 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00031 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00032 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
00033 /*                                                                      */
00034 /************************************************************************/
00035 
00036 #ifndef VIGRA_QUADPROG_HXX
00037 #define VIGRA_QUADPROG_HXX
00038 
00039 #include <limits>
00040 #include "mathutil.hxx"
00041 #include "matrix.hxx"
00042 #include "linear_solve.hxx"
00043 #include "numerictraits.hxx"
00044 #include "array_vector.hxx"
00045 
00046 namespace vigra {
00047 
00048 namespace detail {
00049 
00050 template <class T, class C1, class C2, class C3>
00051 bool quadprogAddConstraint(MultiArrayView<2, T, C1> & R, MultiArrayView<2, T, C2> & J, MultiArrayView<2, T, C3> & d, 
00052                            int activeConstraintCount, double& R_norm)
00053 {
00054     typedef typename MultiArrayShape<2>::type Shape;
00055     int n=columnCount(J);
00056     linalg::detail::qrGivensStepImpl(0, subVector(d, activeConstraintCount, n),
00057                                      J.subarray(Shape(activeConstraintCount,0), Shape(n,n)));
00058     if (abs(d(activeConstraintCount,0)) <= NumericTraits<T>::epsilon() * R_norm) // problem degenerate
00059         return false;
00060     R_norm = std::max<T>(R_norm, abs(d(activeConstraintCount,0)));
00061 
00062     ++activeConstraintCount;   
00063     // add d as a new column to R
00064     columnVector(R, Shape(0, activeConstraintCount - 1), activeConstraintCount) = subVector(d, 0, activeConstraintCount);  
00065     return true;
00066 }
00067 
00068 template <class T, class C1, class C2, class C3>
00069 void quadprogDeleteConstraint(MultiArrayView<2, T, C1> & R, MultiArrayView<2, T, C2> & J, MultiArrayView<2, T, C3> & u, 
00070                               int activeConstraintCount,  int constraintToBeRemoved)
00071 {
00072     typedef typename MultiArrayShape<2>::type Shape;
00073     
00074     int newActiveConstraintCount = activeConstraintCount - 1;
00075 
00076     if(constraintToBeRemoved == newActiveConstraintCount)
00077         return;
00078 
00079     std::swap(u(constraintToBeRemoved,0), u(newActiveConstraintCount,0));
00080     columnVector(R, constraintToBeRemoved).swapData(columnVector(R, newActiveConstraintCount));
00081     linalg::detail::qrGivensStepImpl(0, R.subarray(Shape(constraintToBeRemoved, constraintToBeRemoved), 
00082                                                    Shape(newActiveConstraintCount,newActiveConstraintCount)),
00083                                         J.subarray(Shape(constraintToBeRemoved, 0), 
00084                                                    Shape(newActiveConstraintCount,newActiveConstraintCount)));
00085 }
00086 
00087 } // namespace detail
00088 
00089 /** \addtogroup Optimization Optimization and Regression
00090  */
00091 //@{
00092    /** Solve Quadratic Programming Problem.
00093 
00094      The quadraticProgramming() function implements the algorithm described in
00095      
00096      D. Goldfarb, A. Idnani: <i>"A numerically stable dual method for solving
00097                  strictly convex quadratic programs"</i>, Mathematical Programming 27:1-33, 1983. 
00098      
00099      for the solution of (convex) quadratic programming problems by means of a primal-dual method.
00100          
00101      <b>\#include</b> <<a href="quadprog_8hxx-source.html">vigra/quadprog.hxx</a>>
00102          Namespaces: vigra
00103 
00104      <b>Declaration:</b>
00105 
00106      \code
00107      namespace vigra { 
00108          template <class T, class C1, class C2, class C3, class C4, class C5, class C6, class C7>
00109          T 
00110          quadraticProgramming(MultiArrayView<2, T, C1> const & GG, MultiArrayView<2, T, C2> const & g,  
00111                               MultiArrayView<2, T, C3> const & CE, MultiArrayView<2, T, C4> const & ce,  
00112                               MultiArrayView<2, T, C5> const & CI, MultiArrayView<2, T, C6> const & ci, 
00113                               MultiArrayView<2, T, C7> & x);
00114      }
00115      \endcode
00116 
00117      The problem must be specified in the form:
00118 
00119      \f{eqnarray*}
00120         \mbox{minimize } &\,& \frac{1}{2} \mbox{\bf x}'\,\mbox{\bf G}\, \mbox{\bf x} + \mbox{\bf g}'\,\mbox{\bf x} \\
00121         \mbox{subject to} &\,& \mbox{\bf C}_E\, \mbox{\bf x} = \mbox{\bf c}_e \\
00122          &\,& \mbox{\bf C}_I\,\mbox{\bf x} \ge \mbox{\bf c}_i
00123      \f}            
00124      Matrix <b>G</b> G must be symmetric positive definite, and matrix <b>C</b><sub>E</sub> must have full row rank. 
00125      Matrix and vector dimensions must be as follows:
00126      <ul>
00127      <li> <b>G</b>: [n * n], <b>g</b>: [n * 1]
00128      <li> <b>C</b><sub>E</sub>: [me * n], <b>c</b><sub>e</sub>: [me * 1]
00129      <li> <b>C</b><sub>I</sub>: [mi * n], <b>c</b><sub>i</sub>: [mi * 1]
00130      <li> <b>x</b>: [n * 1]
00131      </ul>
00132      
00133      The function writes the optimal solution into the vector \a x and returns the cost of this solution. 
00134      If the problem is infeasible, std::numeric_limits::infinity() is returned. In this case
00135      the value of vector \a x is undefined.
00136      
00137      <b>Usage:</b>
00138      
00139      Minimize <tt> f = 0.5 * x'*G*x + g'*x </tt> subject to <tt> -1 &lt;= x &lt;= 1</tt>. 
00140      The solution is <tt> x' = [1.0, 0.5, -1.0] </tt> with <tt> f = -22.625</tt>.
00141      \code
00142       double Gdata[] = {13.0, 12.0, -2.0,
00143                         12.0, 17.0,  6.0,
00144                         -2.0,  6.0, 12.0};
00145 
00146       double gdata[] = {-22.0, -14.5, 13.0};
00147 
00148       double CIdata[] = { 1.0,  0.0,  0.0,
00149                           0.0,  1.0,  0.0,
00150                           0.0,  0.0,  1.0,
00151                          -1.0,  0.0,  0.0,
00152                           0.0, -1.0,  0.0,
00153                           0.0,  0.0, -1.0};
00154                         
00155       double cidata[] = {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0};
00156 
00157       Matrix<double> G(3,3, Gdata), 
00158                      g(3,1, gdata), 
00159                      CE,             // empty since there are no equality constraints
00160                      ce,             // likewise
00161                      CI(7,3, CIdata), 
00162                      ci(7,1, cidata), 
00163                      x(3,1);
00164                    
00165       double f = quadraticProgramming(G, g, CE, ce, CI, ci, x);
00166      \endcode
00167    */
00168 template <class T, class C1, class C2, class C3, class C4, class C5, class C6, class C7>
00169 T 
00170 quadraticProgramming(MultiArrayView<2, T, C1> const & G, MultiArrayView<2, T, C2> const & g,  
00171                MultiArrayView<2, T, C3> const & CE, MultiArrayView<2, T, C4> const & ce,  
00172                MultiArrayView<2, T, C5> const & CI, MultiArrayView<2, T, C6> const & ci, 
00173                MultiArrayView<2, T, C7> & x)
00174 {
00175     using namespace linalg;
00176     typedef typename MultiArrayShape<2>::type Shape;
00177     
00178     int n  = rowCount(g),
00179         me = rowCount(ce),
00180         mi = rowCount(ci),
00181         constraintCount = me + mi;
00182         
00183     vigra_precondition(columnCount(G) == n && rowCount(G) == n,
00184         "quadraticProgramming(): Matrix shape mismatch between G and g.");
00185     vigra_precondition(rowCount(x) == n,
00186         "quadraticProgramming(): Output vector x has illegal shape.");
00187     vigra_precondition((me > 0 && columnCount(CE) == n && rowCount(CE) == me) || 
00188                        (me == 0 && columnCount(CE) == 0),
00189         "quadraticProgramming(): Matrix CE has illegal shape.");
00190     vigra_precondition((mi > 0 && columnCount(CI) == n && rowCount(CI) == mi) || 
00191                        (mi == 0 && columnCount(CI) == 0),
00192         "quadraticProgramming(): Matrix CI has illegal shape.");
00193 
00194     Matrix<T> J = identityMatrix<T>(n);
00195     {
00196         Matrix<T> L(G.shape());
00197         choleskyDecomposition(G, L);
00198         // find unconstrained minimizer of the quadratic form  0.5 * x G x + g' x
00199         choleskySolve(L, -g, x);
00200         // compute the inverse of the factorized matrix G^-1, this is the initial value for J
00201         linearSolveLowerTriangular(L, J, J);
00202     }
00203     // current solution value
00204     T f_value = 0.5 * dot(g, x);
00205     
00206     T epsilonZ   = NumericTraits<T>::epsilon() * sq(J.norm(0)),
00207       epsilonPsi = NumericTraits<T>::epsilon() * trace(G)*trace(J)*100.0,
00208       inf        = std::numeric_limits<T>::infinity();
00209     
00210     Matrix<T> R(n, n), r(constraintCount, 1), u(constraintCount,1);
00211     T R_norm = NumericTraits<T>::one();
00212     
00213     // incorporate equality constraints
00214     for (int i=0; i < me; ++i)
00215     {
00216         MultiArrayView<2, T, C3> np = rowVector(CE, i);
00217         Matrix<T> d = J*transpose(np);
00218         Matrix<T> z = transpose(J).subarray(Shape(0, i), Shape(n,n))*subVector(d, i, n);
00219         linearSolveUpperTriangular(R.subarray(Shape(0, 0), Shape(i,i)), 
00220                                    subVector(d, 0, i), 
00221                                    subVector(r, 0, i));
00222         // compute step in primal space so that the constraint becomes satisfied
00223         T step = (squaredNorm(z) <= epsilonZ) // i.e. z == 0
00224                      ? 0.0 
00225                      : (-dot(np, x) + ce(i,0)) / dot(z, np);
00226     
00227         x += step * z;    
00228         u(i,0) = step;
00229         subVector(u, 0, i) -= step * subVector(r, 0, i);
00230         
00231         f_value += 0.5 * sq(step) * dot(z, np);
00232     
00233         vigra_precondition(vigra::detail::quadprogAddConstraint(R, J, d, i, R_norm),
00234             "quadraticProgramming(): Equality constraints are linearly dependent.");
00235     }
00236     int activeConstraintCount = me;
00237   
00238     // determine optimum solution and corresponding active inequality constraints
00239     ArrayVector<int> activeSet(mi);
00240     for (int i = 0; i < mi; ++i)
00241         activeSet[i] = i;
00242 
00243     int constraintToBeAdded;
00244     T ss = 0.0;
00245     for (int i = activeConstraintCount-me; i < mi; ++i)
00246     {
00247         T s = dot(rowVector(CI, activeSet[i]), x) - ci(activeSet[i], 0);
00248         if (s < ss)
00249         {
00250             ss = s;
00251             constraintToBeAdded = i;
00252         }
00253     }
00254 
00255     int iter = 0, maxIter = 10*mi;    
00256     while(iter++ < maxIter)
00257     {        
00258         if (ss >= 0.0)       // all constraints are satisfied
00259             return f_value;  // => solved!
00260 
00261         // determine step direction in the primal space (through J, see the paper)
00262         MultiArrayView<2, T, C5> np = rowVector(CI, activeSet[constraintToBeAdded]);
00263         Matrix<T> d = J*transpose(np);
00264         Matrix<T> z = transpose(J).subarray(Shape(0, activeConstraintCount), Shape(n,n))*subVector(d, activeConstraintCount, n);
00265         
00266         // compute negative of the step direction in the dual space
00267         linearSolveUpperTriangular(R.subarray(Shape(0, 0), Shape(activeConstraintCount,activeConstraintCount)), 
00268                                    subVector(d, 0, activeConstraintCount), 
00269                                    subVector(r, 0, activeConstraintCount));
00270 
00271         // determine minimum step length in primal space such that activeSet[constraintToBeAdded] becomes feasible
00272         T primalStep = (squaredNorm(z) <= epsilonZ) // i.e. z == 0
00273                           ? inf
00274                           : -ss / dot(z, np);
00275       
00276         // determine maximum step length in dual space that doesn't violate dual feasibility
00277         // and the corresponding index
00278         T dualStep = inf; 
00279         int constraintToBeRemoved;
00280         for (int k = me; k < activeConstraintCount; ++k)
00281         {
00282             if (r(k,0) > 0.0)
00283             {
00284                 if (u(k,0) / r(k,0) < dualStep)
00285                 {
00286                     dualStep = u(k,0) / r(k,0);
00287                     constraintToBeRemoved = k;
00288                 }
00289             }
00290         }
00291         
00292         // the step is chosen as the minimum of dualStep and primalStep
00293         T step = std::min(dualStep, primalStep);
00294       
00295         // take step and update matrizes
00296       
00297         if (step == inf)
00298         {
00299             // case (i): no step in primal or dual space possible
00300             return inf; // QPP is infeasible 
00301         }
00302         if (primalStep == inf)
00303         {
00304             // case (ii): step in dual space
00305             subVector(u, 0, activeConstraintCount) -= step * subVector(r, 0, activeConstraintCount);
00306             vigra::detail::quadprogDeleteConstraint(R, J, u, activeConstraintCount, constraintToBeRemoved);
00307             --activeConstraintCount;
00308             std::swap(activeSet[constraintToBeRemoved-me], activeSet[activeConstraintCount-me]);
00309             continue;
00310         }
00311       
00312         // case (iii): step in primal and dual space      
00313         x += step * z;
00314         // update the solution value
00315         f_value += 0.5 * sq(step) * dot(z, np);
00316         // u = [u 1]' + step * [-r 1]
00317         subVector(u, 0, activeConstraintCount) -= step * subVector(r, 0, activeConstraintCount);
00318         u(activeConstraintCount,0) = step;
00319       
00320         if (step == primalStep)
00321         {
00322             // add constraintToBeAdded to the active set
00323             vigra::detail::quadprogAddConstraint(R, J, d, activeConstraintCount, R_norm);
00324             std::swap(activeSet[constraintToBeAdded], activeSet[activeConstraintCount-me]);
00325             ++activeConstraintCount;
00326         }
00327         else
00328         {
00329             // drop constraintToBeRemoved from the active set
00330             vigra::detail::quadprogDeleteConstraint(R, J, u, activeConstraintCount, constraintToBeRemoved);
00331             --activeConstraintCount;
00332             std::swap(activeSet[constraintToBeRemoved-me], activeSet[activeConstraintCount-me]);
00333         }
00334         
00335         // update values of inactive inequality constraints
00336         ss = 0.0;
00337         for (int i = activeConstraintCount-me; i < mi; ++i)
00338         {
00339             // compute CI*x - ci with appropriate row permutation
00340             T s = dot(rowVector(CI, activeSet[i]), x) - ci(activeSet[i], 0);
00341             if (s < ss)
00342             {
00343                 ss = s;
00344                 constraintToBeAdded = i;
00345             }
00346         }
00347     }
00348     return inf; // too many iterations
00349 }
00350 
00351 //@}
00352 
00353 } // namespace vigra
00354 
00355 #endif

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.7.0 (Thu Aug 25 2011)