[ 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> <vigra/quadprog.hxx>
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  inf = std::numeric_limits<T>::infinity();
208 
209  Matrix<T> R(n, n), r(constraintCount, 1), u(constraintCount,1);
210  T R_norm = NumericTraits<T>::one();
211 
212  // incorporate equality constraints
213  for (int i=0; i < me; ++i)
214  {
216  Matrix<T> d = J*transpose(np);
217  Matrix<T> z = transpose(J).subarray(Shape(0, i), Shape(n,n))*subVector(d, i, n);
218  linearSolveUpperTriangular(R.subarray(Shape(0, 0), Shape(i,i)),
219  subVector(d, 0, i),
220  subVector(r, 0, i));
221  // compute step in primal space so that the constraint becomes satisfied
222  T step = (squaredNorm(z) <= epsilonZ) // i.e. z == 0
223  ? 0.0
224  : (-dot(np, x) + ce(i,0)) / dot(z, np);
225 
226  x += step * z;
227  u(i,0) = step;
228  subVector(u, 0, i) -= step * subVector(r, 0, i);
229 
230  f_value += 0.5 * sq(step) * dot(z, np);
231 
232  vigra_precondition(vigra::detail::quadprogAddConstraint(R, J, d, i, R_norm),
233  "quadraticProgramming(): Equality constraints are linearly dependent.");
234  }
235  int activeConstraintCount = me;
236 
237  // determine optimum solution and corresponding active inequality constraints
238  ArrayVector<int> activeSet(mi);
239  for (int i = 0; i < mi; ++i)
240  activeSet[i] = i;
241 
242  int constraintToBeAdded = 0;
243  T ss = 0.0;
244  for (int i = activeConstraintCount-me; i < mi; ++i)
245  {
246  T s = dot(rowVector(CI, activeSet[i]), x) - ci(activeSet[i], 0);
247  if (s < ss)
248  {
249  ss = s;
250  constraintToBeAdded = i;
251  }
252  }
253 
254  int iter = 0, maxIter = 10*mi;
255  while(iter++ < maxIter)
256  {
257  if (ss >= 0.0) // all constraints are satisfied
258  return f_value; // => solved!
259 
260  // determine step direction in the primal space (through J, see the paper)
261  MultiArrayView<2, T, C5> np = rowVector(CI, activeSet[constraintToBeAdded]);
262  Matrix<T> d = J*transpose(np);
263  Matrix<T> z = transpose(J).subarray(Shape(0, activeConstraintCount), Shape(n,n))*subVector(d, activeConstraintCount, n);
264 
265  // compute negative of the step direction in the dual space
266  linearSolveUpperTriangular(R.subarray(Shape(0, 0), Shape(activeConstraintCount,activeConstraintCount)),
267  subVector(d, 0, activeConstraintCount),
268  subVector(r, 0, activeConstraintCount));
269 
270  // determine minimum step length in primal space such that activeSet[constraintToBeAdded] becomes feasible
271  T primalStep = (squaredNorm(z) <= epsilonZ) // i.e. z == 0
272  ? inf
273  : -ss / dot(z, np);
274 
275  // determine maximum step length in dual space that doesn't violate dual feasibility
276  // and the corresponding index
277  T dualStep = inf;
278  int constraintToBeRemoved = 0;
279  for (int k = me; k < activeConstraintCount; ++k)
280  {
281  if (r(k,0) > 0.0)
282  {
283  if (u(k,0) / r(k,0) < dualStep)
284  {
285  dualStep = u(k,0) / r(k,0);
286  constraintToBeRemoved = k;
287  }
288  }
289  }
290 
291  // the step is chosen as the minimum of dualStep and primalStep
292  T step = std::min(dualStep, primalStep);
293 
294  // take step and update matrices
295 
296  if (step == inf)
297  {
298  // case (i): no step in primal or dual space possible
299  return inf; // QPP is infeasible
300  }
301  if (primalStep == inf)
302  {
303  // case (ii): step in dual space
304  subVector(u, 0, activeConstraintCount) -= step * subVector(r, 0, activeConstraintCount);
305  vigra::detail::quadprogDeleteConstraint(R, J, u, activeConstraintCount, constraintToBeRemoved);
306  --activeConstraintCount;
307  std::swap(activeSet[constraintToBeRemoved-me], activeSet[activeConstraintCount-me]);
308  continue;
309  }
310 
311  // case (iii): step in primal and dual space
312  x += step * z;
313  // update the solution value
314  f_value += 0.5 * sq(step) * dot(z, np);
315  // u = [u 1]' + step * [-r 1]
316  subVector(u, 0, activeConstraintCount) -= step * subVector(r, 0, activeConstraintCount);
317  u(activeConstraintCount,0) = step;
318 
319  if (step == primalStep)
320  {
321  // add constraintToBeAdded to the active set
322  vigra::detail::quadprogAddConstraint(R, J, d, activeConstraintCount, R_norm);
323  std::swap(activeSet[constraintToBeAdded], activeSet[activeConstraintCount-me]);
324  ++activeConstraintCount;
325  }
326  else
327  {
328  // drop constraintToBeRemoved from the active set
329  vigra::detail::quadprogDeleteConstraint(R, J, u, activeConstraintCount, constraintToBeRemoved);
330  --activeConstraintCount;
331  std::swap(activeSet[constraintToBeRemoved-me], activeSet[activeConstraintCount-me]);
332  }
333 
334  // update values of inactive inequality constraints
335  ss = 0.0;
336  for (int i = activeConstraintCount-me; i < mi; ++i)
337  {
338  // compute CI*x - ci with appropriate row permutation
339  T s = dot(rowVector(CI, activeSet[i]), x) - ci(activeSet[i], 0);
340  if (s < ss)
341  {
342  ss = s;
343  constraintToBeAdded = i;
344  }
345  }
346  }
347  return inf; // too many iterations
348 }
349 
350 //@}
351 
352 } // namespace vigra
353 
354 #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 Sep 24 2013)