[ 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 <= x <= 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) |
html generated using doxygen and Python
|