2 qlpsolve.cc -- implement Active_constraints, Inactive_iter
4 source file of the LilyPond music typesetter
6 (c) 1996, 1997 Han-Wen Nienhuys <hanwen@stack.nl>
9 try fixed point arithmetic, to speed up lily.
12 #include "qlpsolve.hh"
15 #include "choleski.hh"
17 const Real TOL=1e-2; // roughly 1/10 mm
20 Active_constraints::status() const
22 String s("Active|Inactive [");
23 for (int i=0; i< active.size(); i++) {
24 s += String(active[i]) + " ";
28 for (int i=0; i< inactive.size(); i++) {
29 s += String(inactive[i]) + " ";
37 Active_constraints::OK()
42 assert(active.size() +inactive.size() == opt->cons.size());
43 assert(H.dim() == opt->dim());
44 assert(active.size() == A.rows());
47 for (int i=0; i < opt->cons.size(); i++)
49 for (int i=0; i < active.size(); i++) {
53 for (int i=0; i < inactive.size(); i++) {
57 for (int i=0; i < allcons.size(); i++)
58 assert(allcons[i] == 1);
63 Active_constraints::get_lagrange(Vector gradient)
71 Active_constraints::add(int k)
77 inactive.swap(k,inactive.size()-1);
80 Vector a( opt->cons[cidx] );
84 Vector addrow(Ha.dim());
87 a != 0, so if Ha = O(EPS), then
88 Ha * aH / aHa = O(EPS^2/EPS)
90 if H*a == 0, the constraints are dependent.
92 H -= Matrix(Ha/aHa , Ha);
96 sorry, don't know how to justify this. ..
100 A -= Matrix(A*a, addrow);
101 A.insert_row(addrow,A.rows());
103 WARN << "degenerate constraints";
107 Active_constraints::drop(int k)
109 int q=active.size()-1;
112 inactive.push(active[k]);
118 if (a.norm() > EPS) {
122 Real q = a*opt->quad*a;
124 A -= A*opt->quad*Matrix(a,a/q);
126 WARN << "degenerate constraints";
128 Vector rem_row(A.row(q));
129 assert(rem_row.norm() < EPS);
136 Active_constraints::Active_constraints(Ineq_constrained_qp const *op)
141 for (int i=0; i < op->cons.size(); i++)
143 Choleski_decomposition chol(op->quad);
147 /** Find the optimum which is in the planes generated by the active
151 Active_constraints::find_active_optimum(Vector g)
156 /* *************************************************************** */
159 min_elt_index(Vector v)
161 Real m=INFTY; int idx=-1;
162 for (int i = 0; i < v.dim(); i++){
167 assert(v(i) <= INFTY);
173 /**the numerical solving. Mordecai Avriel, Nonlinear Programming: analysis and methods (1976)
178 This is a "projected gradient" algorithm. Starting from a point x
179 the next point is found in a direction determined by projecting
180 the gradient onto the active constraints. (well, not really the
181 gradient. The optimal solution obeying the active constraints is
182 tried. This is why H = Q^-1 in initialisation) )
187 Ineq_constrained_qp::solve(Vector start) const
192 Active_constraints act(this);
199 Vector gradient=quad*x+lin;
200 // Real fvalue = x*quad*x/2 + lin*x + const_term;
203 Vector last_gradient(gradient);
206 while (iterations++ < MAXITER) {
207 Vector direction= - act.find_active_optimum(gradient);
209 mtor << "gradient "<< gradient<< "\ndirection " << direction<<"\n";
211 if (direction.norm() > EPS) {
212 mtor << act.status() << '\n';
216 Inactive_iter minidx(act);
220 we know the optimum on this "hyperplane". Check if we
221 bump into the edges of the simplex
224 for (Inactive_iter ia(act); ia.ok(); ia++) {
226 if (ia.vec() * direction >= 0)
228 Real alfa= - (ia.vec()*x - ia.rhs())/
229 (ia.vec()*direction);
236 Real unbounded_alfa = 1.0;
237 Real optimal_step = min(minalf, unbounded_alfa);
239 Vector deltax=direction * optimal_step;
241 gradient += optimal_step * (quad * deltax);
243 mtor << "step = " << optimal_step<< " (|dx| = " <<
244 deltax.norm() << ")\n";
246 if (minalf < unbounded_alfa) {
247 /* bumped into an edge. try again, in smaller space. */
248 act.add(minidx.idx());
249 mtor << "adding cons "<< minidx.idx()<<'\n';
252 /*ASSERT: we are at optimal solution for this "plane"*/
257 Vector lagrange_mult=act.get_lagrange(gradient);
258 int m= min_elt_index(lagrange_mult);
260 if (m>=0 && lagrange_mult(m) > 0) {
261 break; // optimal sol.
263 assert(gradient.norm() < EPS) ;
268 mtor << "dropping cons " << m<<'\n';
271 if (iterations >= MAXITER)
272 WARN<<"didn't converge!\n";
274 mtor << ": found " << x<<" in " << iterations <<" iterations\n";