vigra - vigra/quadprog.hxx Source File

[ 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> <vigra/quadprog.hxx>
    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       inf        = std::numeric_limits<T>::infinity();
    00208
    00209     Matrix<T> R(n, n), r(constraintCount, 1), u(constraintCount,1);
    00210     T R_norm = NumericTraits<T>::one();
    00211
    00212     // incorporate equality constraints
    00213     for (int i=0; i < me; ++i)
    00214     {
    00215         MultiArrayView<2, T, C3> np = rowVector(CE, i);
    00216         Matrix<T> d = J*transpose(np);
    00217         Matrix<T> z = transpose(J).subarray(Shape(0, i), Shape(n,n))*subVector(d, i, n);
    00218         linearSolveUpperTriangular(R.subarray(Shape(0, 0), Shape(i,i)),
    00219                                    subVector(d, 0, i),
    00220                                    subVector(r, 0, i));
    00221         // compute step in primal space so that the constraint becomes satisfied
    00222         T step = (squaredNorm(z) <= epsilonZ) // i.e. z == 0
    00223                      ? 0.0
    00224                      : (-dot(np, x) + ce(i,0)) / dot(z, np);
    00225
    00226         x += step * z;
    00227         u(i,0) = step;
    00228         subVector(u, 0, i) -= step * subVector(r, 0, i);
    00229
    00230         f_value += 0.5 * sq(step) * dot(z, np);
    00231
    00232         vigra_precondition(vigra::detail::quadprogAddConstraint(R, J, d, i, R_norm),
    00233             "quadraticProgramming(): Equality constraints are linearly dependent.");
    00234     }
    00235     int activeConstraintCount = me;
    00236
    00237     // determine optimum solution and corresponding active inequality constraints
    00238     ArrayVector<int> activeSet(mi);
    00239     for (int i = 0; i < mi; ++i)
    00240         activeSet[i] = i;
    00241
    00242     int constraintToBeAdded = 0;
    00243     T ss = 0.0;
    00244     for (int i = activeConstraintCount-me; i < mi; ++i)
    00245     {
    00246         T s = dot(rowVector(CI, activeSet[i]), x) - ci(activeSet[i], 0);
    00247         if (s < ss)
    00248         {
    00249             ss = s;
    00250             constraintToBeAdded = i;
    00251         }
    00252     }
    00253
    00254     int iter = 0, maxIter = 10*mi;
    00255     while(iter++ < maxIter)
    00256     {
    00257         if (ss >= 0.0)       // all constraints are satisfied
    00258             return f_value;  // => solved!
    00259
    00260         // determine step direction in the primal space (through J, see the paper)
    00261         MultiArrayView<2, T, C5> np = rowVector(CI, activeSet[constraintToBeAdded]);
    00262         Matrix<T> d = J*transpose(np);
    00263         Matrix<T> z = transpose(J).subarray(Shape(0, activeConstraintCount), Shape(n,n))*subVector(d, activeConstraintCount, n);
    00264
    00265         // compute negative of the step direction in the dual space
    00266         linearSolveUpperTriangular(R.subarray(Shape(0, 0), Shape(activeConstraintCount,activeConstraintCount)),
    00267                                    subVector(d, 0, activeConstraintCount),
    00268                                    subVector(r, 0, activeConstraintCount));
    00269
    00270         // determine minimum step length in primal space such that activeSet[constraintToBeAdded] becomes feasible
    00271         T primalStep = (squaredNorm(z) <= epsilonZ) // i.e. z == 0
    00272                           ? inf
    00273                           : -ss / dot(z, np);
    00274
    00275         // determine maximum step length in dual space that doesn't violate dual feasibility
    00276         // and the corresponding index
    00277         T dualStep = inf;
    00278         int constraintToBeRemoved = 0;
    00279         for (int k = me; k < activeConstraintCount; ++k)
    00280         {
    00281             if (r(k,0) > 0.0)
    00282             {
    00283                 if (u(k,0) / r(k,0) < dualStep)
    00284                 {
    00285                     dualStep = u(k,0) / r(k,0);
    00286                     constraintToBeRemoved = k;
    00287                 }
    00288             }
    00289         }
    00290
    00291         // the step is chosen as the minimum of dualStep and primalStep
    00292         T step = std::min(dualStep, primalStep);
    00293
    00294         // take step and update matrices
    00295
    00296         if (step == inf)
    00297         {
    00298             // case (i): no step in primal or dual space possible
    00299             return inf; // QPP is infeasible
    00300         }
    00301         if (primalStep == inf)
    00302         {
    00303             // case (ii): step in dual space
    00304             subVector(u, 0, activeConstraintCount) -= step * subVector(r, 0, activeConstraintCount);
    00305             vigra::detail::quadprogDeleteConstraint(R, J, u, activeConstraintCount, constraintToBeRemoved);
    00306             --activeConstraintCount;
    00307             std::swap(activeSet[constraintToBeRemoved-me], activeSet[activeConstraintCount-me]);
    00308             continue;
    00309         }
    00310
    00311         // case (iii): step in primal and dual space
    00312         x += step * z;
    00313         // update the solution value
    00314         f_value += 0.5 * sq(step) * dot(z, np);
    00315         // u = [u 1]' + step * [-r 1]
    00316         subVector(u, 0, activeConstraintCount) -= step * subVector(r, 0, activeConstraintCount);
    00317         u(activeConstraintCount,0) = step;
    00318
    00319         if (step == primalStep)
    00320         {
    00321             // add constraintToBeAdded to the active set
    00322             vigra::detail::quadprogAddConstraint(R, J, d, activeConstraintCount, R_norm);
    00323             std::swap(activeSet[constraintToBeAdded], activeSet[activeConstraintCount-me]);
    00324             ++activeConstraintCount;
    00325         }
    00326         else
    00327         {
    00328             // drop constraintToBeRemoved from the active set
    00329             vigra::detail::quadprogDeleteConstraint(R, J, u, activeConstraintCount, constraintToBeRemoved);
    00330             --activeConstraintCount;
    00331             std::swap(activeSet[constraintToBeRemoved-me], activeSet[activeConstraintCount-me]);
    00332         }
    00333
    00334         // update values of inactive inequality constraints
    00335         ss = 0.0;
    00336         for (int i = activeConstraintCount-me; i < mi; ++i)
    00337         {
    00338             // compute CI*x - ci with appropriate row permutation
    00339             T s = dot(rowVector(CI, activeSet[i]), x) - ci(activeSet[i], 0);
    00340             if (s < ss)
    00341             {
    00342                 ss = s;
    00343                 constraintToBeAdded = i;
    00344             }
    00345         }
    00346     }
    00347     return inf; // too many iterations
    00348 }
    00349
    00350 //@}
    00351
    00352 } // namespace vigra
    00353
    00354 #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.9.0 (Tue Nov 6 2012)