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

quadprog.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 2008 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR activeSet PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 #ifndef VIGRA_QUADPROG_HXX
37 #define VIGRA_QUADPROG_HXX
38 
39 #include <limits>
40 #include "mathutil.hxx"
41 #include "matrix.hxx"
42 #include "linear_solve.hxx"
43 #include "numerictraits.hxx"
44 #include "array_vector.hxx"
45 
46 namespace vigra {
47 
48 namespace detail {
49 
50 template <class T, class C1, class C2, class C3>
51 bool quadprogAddConstraint(MultiArrayView<2, T, C1> & R, MultiArrayView<2, T, C2> & J, MultiArrayView<2, T, C3> & d,
52  int activeConstraintCount, double& R_norm)
53 {
54  typedef typename MultiArrayShape<2>::type Shape;
55  int n=columnCount(J);
56  linalg::detail::qrGivensStepImpl(0, subVector(d, activeConstraintCount, n),
57  J.subarray(Shape(activeConstraintCount,0), Shape(n,n)));
58  if (abs(d(activeConstraintCount,0)) <= NumericTraits<T>::epsilon() * R_norm) // problem degenerate
59  return false;
60  R_norm = std::max<T>(R_norm, abs(d(activeConstraintCount,0)));
61 
62  ++activeConstraintCount;
63  // add d as a new column to R
64  columnVector(R, Shape(0, activeConstraintCount - 1), activeConstraintCount) = subVector(d, 0, activeConstraintCount);
65  return true;
66 }
67 
68 template <class T, class C1, class C2, class C3>
69 void quadprogDeleteConstraint(MultiArrayView<2, T, C1> & R, MultiArrayView<2, T, C2> & J, MultiArrayView<2, T, C3> & u,
70  int activeConstraintCount, int constraintToBeRemoved)
71 {
72  typedef typename MultiArrayShape<2>::type Shape;
73 
74  int newActiveConstraintCount = activeConstraintCount - 1;
75 
76  if(constraintToBeRemoved == newActiveConstraintCount)
77  return;
78 
79  std::swap(u(constraintToBeRemoved,0), u(newActiveConstraintCount,0));
80  columnVector(R, constraintToBeRemoved).swapData(columnVector(R, newActiveConstraintCount));
81  linalg::detail::qrGivensStepImpl(0, R.subarray(Shape(constraintToBeRemoved, constraintToBeRemoved),
82  Shape(newActiveConstraintCount,newActiveConstraintCount)),
83  J.subarray(Shape(constraintToBeRemoved, 0),
84  Shape(newActiveConstraintCount,newActiveConstraintCount)));
85 }
86 
87 } // namespace detail
88 
89 /** \addtogroup Optimization Optimization and Regression
90  */
91 //@{
92  /** Solve Quadratic Programming Problem.
93 
94  The quadraticProgramming() function implements the algorithm described in
95 
96  D. Goldfarb, A. Idnani: <i>"A numerically stable dual method for solving
97  strictly convex quadratic programs"</i>, Mathematical Programming 27:1-33, 1983.
98 
99  for the solution of (convex) quadratic programming problems by means of a primal-dual method.
100 
101  <b>\#include</b> <<a href="quadprog_8hxx-source.html">vigra/quadprog.hxx</a>>
102  Namespaces: vigra
103 
104  <b>Declaration:</b>
105 
106  \code
107  namespace vigra {
108  template <class T, class C1, class C2, class C3, class C4, class C5, class C6, class C7>
109  T
110  quadraticProgramming(MultiArrayView<2, T, C1> const & GG, MultiArrayView<2, T, C2> const & g,
111  MultiArrayView<2, T, C3> const & CE, MultiArrayView<2, T, C4> const & ce,
112  MultiArrayView<2, T, C5> const & CI, MultiArrayView<2, T, C6> const & ci,
113  MultiArrayView<2, T, C7> & x);
114  }
115  \endcode
116 
117  The problem must be specified in the form:
118 
119  \f{eqnarray*}
120  \mbox{minimize } &\,& \frac{1}{2} \mbox{\bf x}'\,\mbox{\bf G}\, \mbox{\bf x} + \mbox{\bf g}'\,\mbox{\bf x} \\
121  \mbox{subject to} &\,& \mbox{\bf C}_E\, \mbox{\bf x} = \mbox{\bf c}_e \\
122  &\,& \mbox{\bf C}_I\,\mbox{\bf x} \ge \mbox{\bf c}_i
123  \f}
124  Matrix <b>G</b> G must be symmetric positive definite, and matrix <b>C</b><sub>E</sub> must have full row rank.
125  Matrix and vector dimensions must be as follows:
126  <ul>
127  <li> <b>G</b>: [n * n], <b>g</b>: [n * 1]
128  <li> <b>C</b><sub>E</sub>: [me * n], <b>c</b><sub>e</sub>: [me * 1]
129  <li> <b>C</b><sub>I</sub>: [mi * n], <b>c</b><sub>i</sub>: [mi * 1]
130  <li> <b>x</b>: [n * 1]
131  </ul>
132 
133  The function writes the optimal solution into the vector \a x and returns the cost of this solution.
134  If the problem is infeasible, std::numeric_limits::infinity() is returned. In this case
135  the value of vector \a x is undefined.
136 
137  <b>Usage:</b>
138 
139  Minimize <tt> f = 0.5 * x'*G*x + g'*x </tt> subject to <tt> -1 &lt;= x &lt;= 1</tt>.
140  The solution is <tt> x' = [1.0, 0.5, -1.0] </tt> with <tt> f = -22.625</tt>.
141  \code
142  double Gdata[] = {13.0, 12.0, -2.0,
143  12.0, 17.0, 6.0,
144  -2.0, 6.0, 12.0};
145 
146  double gdata[] = {-22.0, -14.5, 13.0};
147 
148  double CIdata[] = { 1.0, 0.0, 0.0,
149  0.0, 1.0, 0.0,
150  0.0, 0.0, 1.0,
151  -1.0, 0.0, 0.0,
152  0.0, -1.0, 0.0,
153  0.0, 0.0, -1.0};
154 
155  double cidata[] = {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0};
156 
157  Matrix<double> G(3,3, Gdata),
158  g(3,1, gdata),
159  CE, // empty since there are no equality constraints
160  ce, // likewise
161  CI(7,3, CIdata),
162  ci(7,1, cidata),
163  x(3,1);
164 
165  double f = quadraticProgramming(G, g, CE, ce, CI, ci, x);
166  \endcode
167  */
168 template <class T, class C1, class C2, class C3, class C4, class C5, class C6, class C7>
169 T
171  MultiArrayView<2, T, C3> const & CE, MultiArrayView<2, T, C4> const & ce,
172  MultiArrayView<2, T, C5> const & CI, MultiArrayView<2, T, C6> const & ci,
174 {
175  using namespace linalg;
176  typedef typename MultiArrayShape<2>::type Shape;
177 
178  int n = rowCount(g),
179  me = rowCount(ce),
180  mi = rowCount(ci),
181  constraintCount = me + mi;
182 
183  vigra_precondition(columnCount(G) == n && rowCount(G) == n,
184  "quadraticProgramming(): Matrix shape mismatch between G and g.");
185  vigra_precondition(rowCount(x) == n,
186  "quadraticProgramming(): Output vector x has illegal shape.");
187  vigra_precondition((me > 0 && columnCount(CE) == n && rowCount(CE) == me) ||
188  (me == 0 && columnCount(CE) == 0),
189  "quadraticProgramming(): Matrix CE has illegal shape.");
190  vigra_precondition((mi > 0 && columnCount(CI) == n && rowCount(CI) == mi) ||
191  (mi == 0 && columnCount(CI) == 0),
192  "quadraticProgramming(): Matrix CI has illegal shape.");
193 
194  Matrix<T> J = identityMatrix<T>(n);
195  {
196  Matrix<T> L(G.shape());
197  choleskyDecomposition(G, L);
198  // find unconstrained minimizer of the quadratic form 0.5 * x G x + g' x
199  choleskySolve(L, -g, x);
200  // compute the inverse of the factorized matrix G^-1, this is the initial value for J
202  }
203  // current solution value
204  T f_value = 0.5 * dot(g, x);
205 
206  T epsilonZ = NumericTraits<T>::epsilon() * sq(J.norm(0)),
207  epsilonPsi = NumericTraits<T>::epsilon() * trace(G)*trace(J)*100.0,
208  inf = std::numeric_limits<T>::infinity();
209 
210  Matrix<T> R(n, n), r(constraintCount, 1), u(constraintCount,1);
211  T R_norm = NumericTraits<T>::one();
212 
213  // incorporate equality constraints
214  for (int i=0; i < me; ++i)
215  {
217  Matrix<T> d = J*transpose(np);
218  Matrix<T> z = transpose(J).subarray(Shape(0, i), Shape(n,n))*subVector(d, i, n);
219  linearSolveUpperTriangular(R.subarray(Shape(0, 0), Shape(i,i)),
220  subVector(d, 0, i),
221  subVector(r, 0, i));
222  // compute step in primal space so that the constraint becomes satisfied
223  T step = (squaredNorm(z) <= epsilonZ) // i.e. z == 0
224  ? 0.0
225  : (-dot(np, x) + ce(i,0)) / dot(z, np);
226 
227  x += step * z;
228  u(i,0) = step;
229  subVector(u, 0, i) -= step * subVector(r, 0, i);
230 
231  f_value += 0.5 * sq(step) * dot(z, np);
232 
233  vigra_precondition(vigra::detail::quadprogAddConstraint(R, J, d, i, R_norm),
234  "quadraticProgramming(): Equality constraints are linearly dependent.");
235  }
236  int activeConstraintCount = me;
237 
238  // determine optimum solution and corresponding active inequality constraints
239  ArrayVector<int> activeSet(mi);
240  for (int i = 0; i < mi; ++i)
241  activeSet[i] = i;
242 
243  int constraintToBeAdded;
244  T ss = 0.0;
245  for (int i = activeConstraintCount-me; i < mi; ++i)
246  {
247  T s = dot(rowVector(CI, activeSet[i]), x) - ci(activeSet[i], 0);
248  if (s < ss)
249  {
250  ss = s;
251  constraintToBeAdded = i;
252  }
253  }
254 
255  int iter = 0, maxIter = 10*mi;
256  while(iter++ < maxIter)
257  {
258  if (ss >= 0.0) // all constraints are satisfied
259  return f_value; // => solved!
260 
261  // determine step direction in the primal space (through J, see the paper)
262  MultiArrayView<2, T, C5> np = rowVector(CI, activeSet[constraintToBeAdded]);
263  Matrix<T> d = J*transpose(np);
264  Matrix<T> z = transpose(J).subarray(Shape(0, activeConstraintCount), Shape(n,n))*subVector(d, activeConstraintCount, n);
265 
266  // compute negative of the step direction in the dual space
267  linearSolveUpperTriangular(R.subarray(Shape(0, 0), Shape(activeConstraintCount,activeConstraintCount)),
268  subVector(d, 0, activeConstraintCount),
269  subVector(r, 0, activeConstraintCount));
270 
271  // determine minimum step length in primal space such that activeSet[constraintToBeAdded] becomes feasible
272  T primalStep = (squaredNorm(z) <= epsilonZ) // i.e. z == 0
273  ? inf
274  : -ss / dot(z, np);
275 
276  // determine maximum step length in dual space that doesn't violate dual feasibility
277  // and the corresponding index
278  T dualStep = inf;
279  int constraintToBeRemoved;
280  for (int k = me; k < activeConstraintCount; ++k)
281  {
282  if (r(k,0) > 0.0)
283  {
284  if (u(k,0) / r(k,0) < dualStep)
285  {
286  dualStep = u(k,0) / r(k,0);
287  constraintToBeRemoved = k;
288  }
289  }
290  }
291 
292  // the step is chosen as the minimum of dualStep and primalStep
293  T step = std::min(dualStep, primalStep);
294 
295  // take step and update matrizes
296 
297  if (step == inf)
298  {
299  // case (i): no step in primal or dual space possible
300  return inf; // QPP is infeasible
301  }
302  if (primalStep == inf)
303  {
304  // case (ii): step in dual space
305  subVector(u, 0, activeConstraintCount) -= step * subVector(r, 0, activeConstraintCount);
306  vigra::detail::quadprogDeleteConstraint(R, J, u, activeConstraintCount, constraintToBeRemoved);
307  --activeConstraintCount;
308  std::swap(activeSet[constraintToBeRemoved-me], activeSet[activeConstraintCount-me]);
309  continue;
310  }
311 
312  // case (iii): step in primal and dual space
313  x += step * z;
314  // update the solution value
315  f_value += 0.5 * sq(step) * dot(z, np);
316  // u = [u 1]' + step * [-r 1]
317  subVector(u, 0, activeConstraintCount) -= step * subVector(r, 0, activeConstraintCount);
318  u(activeConstraintCount,0) = step;
319 
320  if (step == primalStep)
321  {
322  // add constraintToBeAdded to the active set
323  vigra::detail::quadprogAddConstraint(R, J, d, activeConstraintCount, R_norm);
324  std::swap(activeSet[constraintToBeAdded], activeSet[activeConstraintCount-me]);
325  ++activeConstraintCount;
326  }
327  else
328  {
329  // drop constraintToBeRemoved from the active set
330  vigra::detail::quadprogDeleteConstraint(R, J, u, activeConstraintCount, constraintToBeRemoved);
331  --activeConstraintCount;
332  std::swap(activeSet[constraintToBeRemoved-me], activeSet[activeConstraintCount-me]);
333  }
334 
335  // update values of inactive inequality constraints
336  ss = 0.0;
337  for (int i = activeConstraintCount-me; i < mi; ++i)
338  {
339  // compute CI*x - ci with appropriate row permutation
340  T s = dot(rowVector(CI, activeSet[i]), x) - ci(activeSet[i], 0);
341  if (s < ss)
342  {
343  ss = s;
344  constraintToBeAdded = i;
345  }
346  }
347  }
348  return inf; // too many iterations
349 }
350 
351 //@}
352 
353 } // namespace vigra
354 
355 #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.1 (Wed Mar 12 2014)