36#ifndef VIGRA_QUADPROG_HXX
37#define VIGRA_QUADPROG_HXX
40#include "mathutil.hxx"
42#include "linear_solve.hxx"
43#include "numerictraits.hxx"
44#include "array_vector.hxx"
50template <
class T,
class C1,
class C2,
class C3>
51bool quadprogAddConstraint(MultiArrayView<2, T, C1> & R, MultiArrayView<2, T, C2> & J, MultiArrayView<2, T, C3> & d,
52 int activeConstraintCount,
double& R_norm)
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)
60 R_norm = std::max<T>(R_norm,
abs(d(activeConstraintCount,0)));
62 ++activeConstraintCount;
64 columnVector(R, Shape(0, activeConstraintCount - 1), activeConstraintCount) = subVector(d, 0, activeConstraintCount);
68template <
class T,
class C1,
class C2,
class C3>
69void quadprogDeleteConstraint(MultiArrayView<2, T, C1> & R, MultiArrayView<2, T, C2> & J, MultiArrayView<2, T, C3> & u,
70 int activeConstraintCount,
int constraintToBeRemoved)
74 int newActiveConstraintCount = activeConstraintCount - 1;
76 if(constraintToBeRemoved == newActiveConstraintCount)
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)));
204template <
class T,
class C1,
class C2,
class C3,
class C4,
class C5,
class C6,
class C7>
211 using namespace linalg;
217 constraintCount = me + mi;
219 vigra_precondition(columnCount(G) == n && rowCount(G) == n,
220 "quadraticProgramming(): Matrix shape mismatch between G and g.");
221 vigra_precondition(rowCount(x) == n,
222 "quadraticProgramming(): Output vector x has illegal shape.");
223 vigra_precondition((me > 0 && columnCount(CE) == n && rowCount(CE) == me) ||
224 (me == 0 && columnCount(CE) == 0),
225 "quadraticProgramming(): Matrix CE has illegal shape.");
226 vigra_precondition((mi > 0 && columnCount(CI) == n && rowCount(CI) == mi) ||
227 (mi == 0 && columnCount(CI) == 0),
228 "quadraticProgramming(): Matrix CI has illegal shape.");
233 choleskyDecomposition(G, L);
235 choleskySolve(L, -g, x);
237 linearSolveLowerTriangular(L, J, J);
240 T f_value = 0.5 *
dot(g, x);
242 T epsilonZ = NumericTraits<T>::epsilon() *
sq(J.
norm(0)),
243 inf = std::numeric_limits<T>::infinity();
245 Matrix<T> R(n, n), r(constraintCount, 1), u(constraintCount,1);
246 T R_norm = NumericTraits<T>::one();
249 for (
int i=0; i < me; ++i)
253 Matrix<T> z = transpose(J).subarray(Shape(0, i), Shape(n,n))*subVector(d, i, n);
254 linearSolveUpperTriangular(R.subarray(Shape(0, 0), Shape(i,i)),
260 : (-
dot(np, x) + ce(i,0)) /
dot(z, np);
264 subVector(u, 0, i) -= step * subVector(r, 0, i);
266 f_value += 0.5 *
sq(step) *
dot(z, np);
268 vigra_precondition(vigra::detail::quadprogAddConstraint(R, J, d, i, R_norm),
269 "quadraticProgramming(): Equality constraints are linearly dependent.");
271 int activeConstraintCount = me;
275 for (
int i = 0; i < mi; ++i)
278 int constraintToBeAdded = 0;
280 for (
int i = activeConstraintCount-me; i < mi; ++i)
282 T s =
dot(rowVector(CI, activeSet[i]), x) - ci(activeSet[i], 0);
286 constraintToBeAdded = i;
290 int iter = 0, maxIter = 10*mi;
291 while(iter++ < maxIter)
299 Matrix<T> z = transpose(J).subarray(Shape(0, activeConstraintCount), Shape(n,n))*subVector(d, activeConstraintCount, n);
302 linearSolveUpperTriangular(R.subarray(Shape(0, 0), Shape(activeConstraintCount,activeConstraintCount)),
303 subVector(d, 0, activeConstraintCount),
304 subVector(r, 0, activeConstraintCount));
314 int constraintToBeRemoved = 0;
315 for (
int k = me; k < activeConstraintCount; ++k)
319 if (u(k,0) / r(k,0) < dualStep)
321 dualStep = u(k,0) / r(k,0);
322 constraintToBeRemoved = k;
328 T step = std::min(dualStep, primalStep);
337 if (primalStep == inf)
340 subVector(u, 0, activeConstraintCount) -= step * subVector(r, 0, activeConstraintCount);
341 vigra::detail::quadprogDeleteConstraint(R, J, u, activeConstraintCount, constraintToBeRemoved);
342 --activeConstraintCount;
343 std::swap(activeSet[constraintToBeRemoved-me], activeSet[activeConstraintCount-me]);
350 f_value += 0.5 *
sq(step) *
dot(z, np);
352 subVector(u, 0, activeConstraintCount) -= step * subVector(r, 0, activeConstraintCount);
353 u(activeConstraintCount,0) = step;
355 if (step == primalStep)
358 vigra::detail::quadprogAddConstraint(R, J, d, activeConstraintCount, R_norm);
359 std::swap(activeSet[constraintToBeAdded], activeSet[activeConstraintCount-me]);
360 ++activeConstraintCount;
365 vigra::detail::quadprogDeleteConstraint(R, J, u, activeConstraintCount, constraintToBeRemoved);
366 --activeConstraintCount;
367 std::swap(activeSet[constraintToBeRemoved-me], activeSet[activeConstraintCount-me]);
372 for (
int i = activeConstraintCount-me; i < mi; ++i)
375 T s =
dot(rowVector(CI, activeSet[i]), x) - ci(activeSet[i], 0);
379 constraintToBeAdded = i;
Definition array_vector.hxx:514
TinyVector< MultiArrayIndex, N > type
Definition multi_shape.hxx:272
Base class for, and view to, vigra::MultiArray.
Definition multi_fwd.hxx:127
const difference_type & shape() const
Definition multi_array.hxx:1650
Class for fixed size vectors.
Definition tinyvector.hxx:1008
Definition matrix.hxx:125
NormTraits< Matrix >::NormType norm() const
unsigned int quadraticProgramming(...)
NumericTraits< T >::Promote sq(T t)
The square function.
Definition mathutil.hxx:382
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition fftw3.hxx:1002
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition fftw3.hxx:1044
PromoteTraits< V1, V2 >::Promote dot(RGBValue< V1, RIDX1, GIDX1, BIDX1 > const &r1, RGBValue< V2, RIDX2, GIDX2, BIDX2 > const &r2)
dot product
Definition rgbvalue.hxx:906