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"
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)
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);
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)
74 int newActiveConstraintCount = activeConstraintCount - 1;
76 if(constraintToBeRemoved == newActiveConstraintCount)
79 std::swap(u(constraintToBeRemoved,0), u(newActiveConstraintCount,0));
81 linalg::detail::qrGivensStepImpl(0, R.subarray(Shape(constraintToBeRemoved, constraintToBeRemoved),
82 Shape(newActiveConstraintCount,newActiveConstraintCount)),
83 J.subarray(Shape(constraintToBeRemoved, 0),
84 Shape(newActiveConstraintCount,newActiveConstraintCount)));
204 template <
class T,
class C1,
class C2,
class C3,
class C4,
class C5,
class C6,
class C7>
207 MultiArrayView<2, T, C3>
const & CE, MultiArrayView<2, T, C4>
const & ce,
208 MultiArrayView<2, T, C5>
const & CI, MultiArrayView<2, T, C6>
const & ci,
209 MultiArrayView<2, T, C7> & x)
211 using namespace linalg;
217 constraintCount = me + mi;
220 "quadraticProgramming(): Matrix shape mismatch between G and g.");
221 vigra_precondition(
rowCount(x) == n,
222 "quadraticProgramming(): Output vector x has illegal shape.");
225 "quadraticProgramming(): Matrix CE has illegal shape.");
228 "quadraticProgramming(): Matrix CI has illegal shape.");
230 Matrix<T> J = identityMatrix<T>(n);
232 Matrix<T> L(G.shape());
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)
251 MultiArrayView<2, T, C3> np =
rowVector(CE, i);
260 : (-
dot(np, x) + ce(i,0)) /
dot(z, np);
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;
274 ArrayVector<int> activeSet(mi);
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)
297 MultiArrayView<2, T, C5> np =
rowVector(CI, activeSet[constraintToBeAdded]);
299 Matrix<T> z =
transpose(J).subarray(Shape(0, activeConstraintCount), Shape(n,n))*
subVector(d, activeConstraintCount, n);
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;
MultiArrayView< 2, T, C > columnVector(MultiArrayView< 2, T, C > const &m, MultiArrayIndex d)
Definition: matrix.hxx:727
MultiArrayIndex rowCount(const MultiArrayView< 2, T, C > &x)
Definition: matrix.hxx:671
void transpose(const MultiArrayView< 2, T, C1 > &v, MultiArrayView< 2, T, C2 > &r)
Definition: matrix.hxx:965
bool linearSolveLowerTriangular(const MultiArrayView< 2, T, C1 > &l, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > x)
Definition: linear_solve.hxx:1118
void choleskySolve(MultiArrayView< 2, T, C1 > const &L, MultiArrayView< 2, T, C2 > const &b, MultiArrayView< 2, T, C3 > &x)
Definition: linear_solve.hxx:1168
linalg::TemporaryMatrix< T > sq(MultiArrayView< 2, T, C > const &v)
doxygen_overloaded_function(template<...> void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
Matrix< T, ALLLOC >::SquaredNormType squaredNorm(const Matrix< T, ALLLOC > &a)
TinyVector< MultiArrayIndex, N > type
Definition: multi_shape.hxx:272
MultiArrayView< 2, T, C > subVector(MultiArrayView< 2, T, C > const &m, int first, int end)
Definition: matrix.hxx:761
NormTraits< T >::SquaredNormType dot(const MultiArrayView< 2, T, C1 > &x, const MultiArrayView< 2, T, C2 > &y)
Definition: matrix.hxx:1342
MultiArrayView< 2, T, C > rowVector(MultiArrayView< 2, T, C > const &m, MultiArrayIndex d)
Definition: matrix.hxx:697
MultiArrayIndex columnCount(const MultiArrayView< 2, T, C > &x)
Definition: matrix.hxx:684
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
bool choleskyDecomposition(MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > &L)
Definition: linear_solve.hxx:959
bool linearSolveUpperTriangular(const MultiArrayView< 2, T, C1 > &r, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > x)
Definition: linear_solve.hxx:1068
unsigned int quadraticProgramming(...)