1 /* ========================================================================== */
2 /* === Source/Mongoose_QPNapsack.cpp ======================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * Mongoose Graph Partitioning Library Copyright (C) 2017-2018,
7 * Scott P. Kolodziej, Nuri S. Yeralan, Timothy A. Davis, William W. Hager
8 * Mongoose is licensed under Version 3 of the GNU General Public License.
9 * Mongoose is also available under other licenses; contact authors for details.
10 * -------------------------------------------------------------------------- */
11
12 /* ========================================================================== */
13 /* === QPNapsack ============================================================ */
14 /* ========================================================================== */
15 /* Find x that minimizes ||x-y|| while satisfying 0 <= x <= 1,
16 a'x = b, lo <= b <= hi. It is assumed that the column vector a is strictly
17 positive since, in our application, the vector a is the vertex weights,
18 which are >= 1. If a is NULL, then it is assumed that a is identically 1.
19 The approach is to solve the dual problem obtained by introducing
20 a multiplier lambda for the constraint a'x = b. The dual function is
21
22 L (lambda) = min { ||x-y||^2 + lambda (a'x - b): 0 <= x <= 1, lo <= b <= hi}
23
24 The dual function is concave. It is continuously differentiable
25 except at lambda = 0. If mu denotes the maximizer of the dual function,
26 then the solution of the primal problem is
27
28 x = proj (y - mu*a) ,
29
30 where proj (z) is the projection of z onto the set { x : 0 <= x <= 1}.
31 Thus we have
32
33 proj (z)_i = 1 if z_i >= 1,
34 0 if z_i <= 0,
35 z_i otherwise .
36
37 Note that for any lambda, the minimizing x in the dual function is
38
39 x (lambda) = proj (y - lambda*a).
40
41 The slope of the dual function is
42
43 L'(lambda) = a'proj (x(lambda)) - hi (if lambda > 0)
44 a'proj (x(lambda)) - lo (if lambda < 0)
45
46 The minimizing b in the dual function is b = hi if lambda > 0 and b = lo
47 if b <= 0. When L' (lamdbda) = 0 with lambda != 0, either x'a = hi or
48 x'a = lo. The minimum is attained at lambda = 0 if and only if the
49 slope of L is negative at lambda = 0+ and positive at lambda = 0-.
50 This is equivalent to the inequalities
51
52 lo <= a' proj (y) <= hi .
53
54 The solution technique is to start with an initial guess lambda for
55 mu and search for a zero of L'. We have the following cases:
56
57 1. lambda >= 0, L'(lambda+) >= 0: mu >= lambda. If L' = 0, then done.
58 Otherwise, increase lambda using napup until
59 slope vanishes
60
61 2. lambda <= 0, L'(lambda-) <= 0: mu <= lambda. If L' = 0, then done.
62 Otherwise, decrease lambda using napdown
63 until slope vanishes
64
65 3. lambda >= 0, L'(lambda+) < 0: If L' (0-) < 0, then mu < 0. Call napdown
66 with lambda = 0 as starting guess. If
67 L' (0+) > 0, then 0 < mu < lambda. Call
68 napdown with given starting guess lambda.
69 Otherwise, if L' (0+) <= 0, then mu = 0.
70
71 4. lambda <= 0, L'(lambda-) > 0: If L' (0+) > 0, then mu > 0. Call napup
72 with lambda = 0 as starting guess. If
73 L' (0-) < 0, then lambda < mu < 0. Call
74 napup with given starting guess lambda.
75 Otherwise, if L' (0-) >= 0, then mu = 0.
76
77 By the "free set" we mean those i for which 0 < x_i (lambda) < 1. The
78 total time taken by napsack is O (n + h log n), where n is the size of y,
79 h is the number of times an element of x (lambda) moves off the boundary
80 into the free set (entries between zero and one) plus the number of times
81 elements move from the free set to the opposite boundary. A heap is used
82 to hold the entries in the boundary and in the free set. If the slope
83 vanishes at either the starting lambda or at lambda = 0, then no heap is
84 constructed, and the time is just O (n).
85
86 If we have a guess for which components of x will be free at the optimal
87 solution, then we can obtain a good guess for the starting lambda by
88 setting the slope of the dual function to zero and solving for lambda. If
89 FreeSet_status is not NULL, then the FreeSet_status array is used to
90 compute a starting guess for lambda based on the estimated free indices.
91 Note that FreeSet_status is an INPUT array, it is not modified by this
92 routine.
93 ========================================================================== */
94
95 #include "Mongoose_Debug.hpp"
96 #include "Mongoose_Internal.hpp"
97 #include "Mongoose_Logger.hpp"
98 #include "Mongoose_QPNapDown.hpp"
99 #include "Mongoose_QPNapUp.hpp"
100
101 #include <cfloat>
102
103 namespace Mongoose
104 {
105
106 #ifndef NDEBUG
checkatx(double * x,double * a,Int n,double lo,double hi,double tol)107 void checkatx(double *x, double *a, Int n, double lo, double hi, double tol)
108 {
109 double atx = 0.;
110 int ok = 1;
111 for (Int k = 0; k < n; k++)
112 {
113 if (x[k] < 0.)
114 {
115 ok = 0;
116 PR(("x [%ld] = %g < 0!\n", k, x[k]));
117 }
118 if (x[k] > 1.)
119 {
120 ok = 0;
121 PR(("x [%ld] = %g > 1!\n", k, x[k]));
122 }
123 if (a != NULL)
124 {
125 double ak = (a) ? a[k] : 1;
126 PR(("a'x = %g * %g = %g\n", ak, x[k], ak * x[k]));
127 atx += ak * x[k];
128 }
129 else
130 {
131 PR(("a'x = %g * %g = %g\n", 1.0, x[k], x[k]));
132 atx += x[k];
133 }
134 }
135 if (atx < lo - tol)
136 {
137 ok = 0;
138 }
139 if (atx > hi + tol)
140 {
141 ok = 0;
142 }
143 if (!ok)
144 {
145 PR(("tol = %g\n", tol));
146 PR(("napsack error! lo %g a'x %g hi %g\n", lo, atx, hi));
147 FFLUSH;
148 ASSERT(0);
149 }
150 }
151 #endif
152
QPNapsack(double * x,Int n,double lo,double hi,double * Gw,double Lambda,const Int * FreeSet_status,double * w,Int * heap1,Int * heap2,double tol)153 double QPNapsack /* return the final lambda */
154 (double *x, /* holds y on input, and the solution x on output */
155 Int n, /* size of x, constraint lo <= a'x <= hi */
156 double lo, /* partition lower bound */
157 double hi, /* partition upper bound */
158 double *Gw, /* vector of nodal weights */
159 double Lambda, /* initial guess for lambda */
160 const Int *FreeSet_status,
161 /* FreeSet_status [i] = +1,-1, or 0 on input,
162 for 3 cases: x_i =1,0, or 0< x_i< 1. Not modified. */
163 double *w, /* work array of size n */
164 Int *heap1, /* work array of size n+1 */
165 Int *heap2, /* work array of size n+1 */
166 double tol /* Gradient projection tolerance */
167 )
168 {
169 (void)tol; // unused variable except during debug
170 double lambda = Lambda;
171 PR(("QPNapsack start [\n"));
172
173 /* ---------------------------------------------------------------------- */
174 /* compute starting guess if FreeSet_status is provided and lambda != 0 */
175 /* ---------------------------------------------------------------------- */
176
177 if ((FreeSet_status != NULL) && (lambda != 0))
178 {
179 double asum = (lambda > 0 ? -hi : -lo);
180 double a2sum = 0.;
181
182 for (Int k = 0; k < n; k++)
183 {
184 if (FreeSet_status[k] == 1)
185 {
186 asum += (Gw) ? Gw[k] : 1;
187 }
188 else if (FreeSet_status[k] == 0)
189 {
190 double ai = (Gw) ? Gw[k] : 1;
191 asum += x[k] * ai;
192 a2sum += ai * ai;
193 }
194 }
195
196 if (a2sum != 0.)
197 lambda = asum / a2sum;
198 }
199
200 /* ---------------------------------------------------------------------- */
201 /* compute the initial slope */
202 /* ---------------------------------------------------------------------- */
203
204 double slope = 0;
205 for (Int k = 0; k < n; k++)
206 {
207 double xi = x[k] - ((Gw) ? Gw[k] : 1) * lambda;
208 if (xi >= 1.)
209 {
210 slope += ((Gw) ? Gw[k] : 1);
211 }
212 else if (xi > 0.)
213 {
214 slope += ((Gw) ? Gw[k] : 1) * xi;
215 }
216 }
217 PR(("slope %g lo %g hi %g\n", slope, lo, hi));
218
219 /* remember: must still adjust slope by "-hi" or "-lo" for its final value
220 */
221
222 if ((lambda >= 0.) && (slope >= hi)) /* case 1 */
223 {
224 if (slope > hi)
225 {
226 PR(("napsack case 1 up\n"));
227 lambda = QPNapUp(x, n, lambda, Gw, hi, w, heap1, heap2);
228 lambda = std::max(0., lambda);
229 }
230 else
231 {
232 PR(("napsack case 1 nothing\n"));
233 }
234 }
235 else if ((lambda <= 0.) && (slope <= lo)) /* case 2 */
236 {
237 if (slope < lo)
238 {
239 PR(("napsack case 2 down\n"));
240 lambda = QPNapDown(x, n, lambda, Gw, lo, w, heap1, heap2);
241 lambda = std::min(lambda, 0.);
242 }
243 else
244 {
245 PR(("napsack case 2 nothing\n"));
246 }
247 }
248 else /* case 3 or 4 */
249 {
250 if (lambda != 0.)
251 {
252 double slope0 = 0.;
253 for (Int k = 0; k < n; k++)
254 {
255 double xi = x[k];
256 if (xi >= 1.)
257 {
258 slope0 += ((Gw) ? Gw[k] : 1);
259 }
260 else if (xi > 0.)
261 {
262 slope0 += ((Gw) ? Gw[k] : 1) * xi;
263 }
264 }
265
266 if ((lambda >= 0) && (slope < hi)) /* case 3 */
267 {
268 if (slope0 < lo)
269 {
270 PR(("napsack case 3a down\n"));
271 lambda = 0.;
272 lambda = QPNapDown(x, n, lambda, Gw, lo, w, heap1, heap2);
273 if (lambda > 0.)
274 {
275 lambda = 0.;
276 }
277 }
278 else if (slope0 > hi)
279 {
280 PR(("napsack case 3b down\n"));
281 lambda = QPNapDown(x, n, lambda, Gw, hi, w, heap1, heap2);
282 if (lambda < 0.)
283 lambda = 0.;
284 }
285 else
286 {
287 PR(("napsack case 3c nothing\n"));
288 lambda = 0.;
289 }
290 }
291 else /* ( (lambda <= 0) && (slope > lo) ) case 4 */
292 {
293 if (slope0 > hi)
294 {
295 PR(("napsack case 4a up\n"));
296 lambda = 0.;
297 lambda = QPNapUp(x, n, lambda, Gw, hi, w, heap1, heap2);
298 lambda = std::max(lambda, 0.);
299 }
300 else if (slope0 < lo)
301 {
302 PR(("napsack case 4b up\n"));
303 lambda = QPNapUp(x, n, lambda, Gw, lo, w, heap1, heap2);
304 lambda = std::min(0., lambda);
305 }
306 else
307 {
308 PR(("napsack case 4c nothing\n"));
309 lambda = 0.;
310 }
311 }
312 }
313 else /* lambda == 0 */
314 {
315 if (slope < hi) /* case 3 */
316 {
317 if (slope < lo)
318 {
319 PR(("napsack case 3d down\n"));
320 lambda = QPNapDown(x, n, lambda, Gw, lo, w, heap1, heap2);
321 lambda = std::min(0., lambda);
322 }
323 else
324 {
325 PR(("napsack case 3e nothing\n"));
326 }
327 }
328 else /* ( slope > lo ) case 4 */
329 {
330 if (slope > hi)
331 {
332 PR(("napsack case 4d up\n"));
333 lambda = QPNapUp(x, n, lambda, Gw, hi, w, heap1, heap2);
334 lambda = std::max(lambda, 0.);
335 }
336 else
337 {
338 PR(("napsack case 4e nothing\n"));
339 }
340 }
341 }
342 }
343
344 /* ---------------------------------------------------------------------- */
345 /* replace x by x (lambda) */
346 /* ---------------------------------------------------------------------- */
347
348 PR(("lambda %g\n", lambda));
349 double atx = 0;
350 Int last_move = 0;
351 for (Int k = 0; k < n; k++)
352 {
353 double xi = x[k] - ((Gw) ? Gw[k] : 1) * lambda;
354
355 if (xi < 0)
356 {
357 x[k] = 0;
358 }
359 else if (xi > 1)
360 {
361 x[k] = 1;
362 }
363 else
364 {
365 x[k] = xi;
366 last_move = k;
367 }
368
369 double newatx = atx + ((Gw) ? Gw[k] : 1) * x[k];
370
371 // Correction step if we go too far
372 if (newatx > hi)
373 {
374 double diff = hi - atx - FLT_MIN;
375 // Need diff = Gw[k] * x[k], so...
376 x[k] = diff / ((Gw) ? Gw[k] : 1);
377 newatx = atx + ((Gw) ? Gw[k] : 1) * x[k];
378 }
379 atx = newatx;
380 }
381
382 // Correction step if we didn't go far enough
383 for (Int kk = 0; kk < n && atx < lo; kk++)
384 {
385 Int k = last_move;
386 atx -= ((Gw) ? Gw[k] : 1) * x[k];
387 double diff = lo - atx;
388 // Need diff = Gw[k] * x[k], so...
389 x[k] = std::min(1., diff / ((Gw) ? Gw[k] : 1));
390 atx += ((Gw) ? Gw[k] : 1) * x[k];
391 last_move = (k + 1) % n;
392 }
393
394 #ifndef NDEBUG
395 // Define check tolerance by lambda values
396 double atx_tol = log10(std::max(fabs(lambda), fabs(Lambda))
397 / (1e-9 + std::min(fabs(lambda), fabs(Lambda))));
398 atx_tol = std::max(atx_tol, tol);
399
400 checkatx(x, Gw, n, lo, hi, atx_tol);
401 #endif
402
403 PR(("QPNapsack done ]\n"));
404
405 return lambda;
406 }
407
408 } // end namespace Mongoose
409