1 | /* $Id: exprQuad.hpp 154 2009-06-16 18:52:53Z pbelotti $ */ |
---|
2 | /** |
---|
3 | * Name: exprQuad.hpp |
---|
4 | * Author: Pietro Belotti |
---|
5 | * Purpose: definition of quadratic expressions (= exprGroup + |
---|
6 | * quadratic = constant + linear + [nonlinear] + quadratic) |
---|
7 | * |
---|
8 | * (C) Carnegie-Mellon University, 2006. |
---|
9 | * This file is licensed under the Common Public License (CPL) |
---|
10 | */ |
---|
11 | |
---|
12 | #ifndef COUENNE_EXPRQUAD_H |
---|
13 | #define COUENNE_EXPRQUAD_H |
---|
14 | |
---|
15 | #include <map> |
---|
16 | #include <vector> |
---|
17 | |
---|
18 | #include "CoinPackedVector.hpp" |
---|
19 | #include "exprGroup.hpp" |
---|
20 | |
---|
21 | class quadElem; |
---|
22 | class CouenneProblem; |
---|
23 | class Domain; |
---|
24 | |
---|
25 | /** class exprQuad, with constant, linear and quadratic terms |
---|
26 | * |
---|
27 | * It represents an expression of the form \f$a_0 + \sum_{i\in I} b_i |
---|
28 | * x_i + x^T Q x + \sum_{i \in J} h_i (x)\f$, with \f$a_0 + \sum_{i\in |
---|
29 | * I} b_i x_i\f$ an affine term, \f$x^T Q x\f$ a quadratic term, and |
---|
30 | * a nonlinear sum \f$\sum_{i \in J} h_i (x)\f$. Standardization |
---|
31 | * checks possible quadratic or linear terms in the latter and |
---|
32 | * includes them in the former parts. |
---|
33 | * |
---|
34 | * If \f$h_i(x)\f$ is a product of two nonlinear, nonquadratic |
---|
35 | * functions \f$h'(x)h''(x)\f$, two auxiliary variables |
---|
36 | * \f$w'=f'(x)\f$ and \f$w''=h''(x)\f$ are created and the product |
---|
37 | * \f$w'w''\f$ is included in the quadratic part of the exprQuad. If |
---|
38 | * \f$h(x)\f$ nonquadratic, nonlinear function, an auxiliary variable |
---|
39 | * \f$w=h(x)\f$ is created and included in the linear part. |
---|
40 | */ |
---|
41 | |
---|
42 | class exprQuad: public exprGroup { |
---|
43 | |
---|
44 | public: |
---|
45 | |
---|
46 | /// matrix |
---|
47 | typedef std::vector <std::pair <exprVar *, CouNumber> > sparseQcol; |
---|
48 | typedef std::vector <std::pair <exprVar *, sparseQcol> > sparseQ; |
---|
49 | |
---|
50 | protected: |
---|
51 | |
---|
52 | /** \name Q matrix storage |
---|
53 | * Sparse implementation: given expression of the form \f$\sum_{i \in N, |
---|
54 | * j \in N} q_{ij} x_i x_j\f$, qindexI_ and qindexJ_ contain |
---|
55 | * respectively entries \f$i\f$ and \f$j\f$ for which \f$q_{ij}\f$ is nonzero |
---|
56 | * in \f$q_{ij} x_i x_j\f$: */ |
---|
57 | /** @{ */ |
---|
58 | |
---|
59 | mutable sparseQ matrix_; |
---|
60 | |
---|
61 | /** \name Convexification data structures |
---|
62 | * |
---|
63 | * These are filled by alphaConvexify, which implements the |
---|
64 | * alpha-convexification method described in the LaGO paper by |
---|
65 | * Nowak and Vigerske -- see also Adjiman and Floudas. |
---|
66 | */ |
---|
67 | |
---|
68 | /// eigenvalues and eigenvectors |
---|
69 | mutable std::vector <std::pair <CouNumber, |
---|
70 | std::vector <std::pair <exprVar *, |
---|
71 | CouNumber> > > > eigen_; |
---|
72 | |
---|
73 | /// current bounds (checked before re-computing eigenvalues/vectors) |
---|
74 | std::map <exprVar *, std::pair <CouNumber, CouNumber> > bounds_; |
---|
75 | |
---|
76 | /// number of non-zeroes in Q |
---|
77 | int nqterms_; |
---|
78 | |
---|
79 | public: |
---|
80 | |
---|
81 | /// Constructor |
---|
82 | exprQuad (CouNumber c0, |
---|
83 | std::vector <std::pair <exprVar *, CouNumber> > &lcoeff, |
---|
84 | std::vector <quadElem> &qcoeff, |
---|
85 | expression **al = NULL, |
---|
86 | int n = 0); |
---|
87 | |
---|
88 | /// Copy constructor |
---|
89 | exprQuad (const exprQuad &src, Domain *d = NULL); |
---|
90 | |
---|
91 | // get indices and coefficients vectors of the quadratic part |
---|
92 | sparseQ &getQ () const |
---|
93 | {return matrix_;} |
---|
94 | |
---|
95 | int getnQTerms () ///< Get number of quadratic terms |
---|
96 | {return nqterms_;} |
---|
97 | |
---|
98 | /// cloning method |
---|
99 | virtual expression *clone (Domain *d = NULL) const |
---|
100 | {return new exprQuad (*this, d);} |
---|
101 | |
---|
102 | /// Print expression to an iostream |
---|
103 | virtual void print (std::ostream & = std::cout, bool = false) const; |
---|
104 | |
---|
105 | /// Function for the evaluation of the expression |
---|
106 | virtual CouNumber operator () (); |
---|
107 | |
---|
108 | /// return l-2 norm of gradient at given point |
---|
109 | CouNumber gradientNorm (const double *x); |
---|
110 | |
---|
111 | /// Compute derivative of this expression with respect to variable |
---|
112 | /// whose index is passed as argument |
---|
113 | virtual expression *differentiate (int index); |
---|
114 | |
---|
115 | /// Simplify expression |
---|
116 | virtual expression *simplify (); |
---|
117 | |
---|
118 | /// Get a measure of "how linear" the expression is |
---|
119 | virtual int Linearity () { |
---|
120 | int |
---|
121 | lin = exprSum::Linearity (), // >= NONLINEAR, |
---|
122 | lin2 = (matrix_ .size () > 0) ? QUADRATIC : |
---|
123 | (lcoeff_ .size () > 0) ? LINEAR : |
---|
124 | (fabs (c0_) > COUENNE_EPS) ? CONSTANT : ZERO; |
---|
125 | |
---|
126 | return ((lin > lin2) ? lin : lin2); |
---|
127 | } |
---|
128 | |
---|
129 | /// Get lower and upper bound of an expression (if any) |
---|
130 | virtual void getBounds (expression *&, expression *&); |
---|
131 | |
---|
132 | /// Get lower and upper bound of an expression (if any) |
---|
133 | virtual void getBounds (CouNumber &, CouNumber &); |
---|
134 | |
---|
135 | /// Generate cuts for the quadratic expression, which are supporting |
---|
136 | /// hyperplanes of the concave upper envelope and the convex lower |
---|
137 | /// envelope. |
---|
138 | virtual void generateCuts (expression *w, const OsiSolverInterface &si, |
---|
139 | OsiCuts &cs, const CouenneCutGenerator *cg, |
---|
140 | t_chg_bounds * = NULL, int = -1, |
---|
141 | CouNumber = -COUENNE_INFINITY, |
---|
142 | CouNumber = COUENNE_INFINITY); |
---|
143 | |
---|
144 | /// Compute data for \f$\alpha\f$-convexification of a quadratic form |
---|
145 | /// (fills in dCoeff_ and dIndex_ for the convex underestimator) |
---|
146 | virtual bool alphaConvexify (const CouenneProblem *, const OsiSolverInterface &); |
---|
147 | |
---|
148 | /** method exprQuad::quadCuts |
---|
149 | * |
---|
150 | * \brief Based on the information (dIndex_, dCoeffLo_, dCoeffUp_) |
---|
151 | * created/modified by alphaConvexify(), create convexification cuts |
---|
152 | * for this expression. |
---|
153 | * |
---|
154 | * The original constraint is : |
---|
155 | * \f[ |
---|
156 | * \eta = a_0 + a^T x + x^T Q x |
---|
157 | * \f] |
---|
158 | * where \f$ \eta \f$ is the auxiliary corresponding to this |
---|
159 | * expression and \f$ w_j \f$ are the auxiliaries corresponding to |
---|
160 | * the other non-linear terms contained in the expression. |
---|
161 | * |
---|
162 | * The under-estimator of \f$ x^T Q x\f$ is given by \f[ x^T Q x + |
---|
163 | * \sum \lambda_{\min,i} (x_i - l_i ) (u_i - x_i ) \f] and its |
---|
164 | * over-estimator is given by |
---|
165 | * |
---|
166 | * \f[ x^T Q x + \sum \lambda_{\max, i} (x_i - l_i ) (u_i - x_i ) |
---|
167 | * \f] (where \f$ \lambda_{\min, i} = \frac{\lambda_{\min}}{w_i^2} |
---|
168 | * \f$, \f$ \lambda_{\max, i} = \frac{\lambda_{\max}}{w_i^2} \f$, |
---|
169 | * and \f$w_i = u_i - l_i\f$), where \f$\lambda_{\max}\f$ |
---|
170 | * (\f$\lambda_{\max}\f$) is the minimum (maximum) eigenvalue of the |
---|
171 | * matrix \f$A={\rm Diag}({\bf u} - {\bf l}) Q {\rm Diag}({\bf u} - |
---|
172 | * {\bf l})\f$, obtained by pre- and post-multiplying \f$ Q \f$ by |
---|
173 | * the diagonal matrix whose \f$i\f$-th element is \f$u_i - l_i\f$. |
---|
174 | * |
---|
175 | * Let \f$ \tilde a_0(\lambda)\f$, \f$ \tilde a(\lambda) \f$ and |
---|
176 | * \f$ \tilde Q (\lambda) \f$ be |
---|
177 | * |
---|
178 | * \f[ \tilde a_0(\lambda) = a_0 - \sum_{i = 1}^n \lambda_i l_i u_i \f] |
---|
179 | * |
---|
180 | * \f[ \tilde a(\lambda) = a + \left[ \begin{array}{c} \lambda_1 |
---|
181 | * (u_1 + l_1) \\ \vdots \\ \lambda_n (u_n + l_n) \end{array} |
---|
182 | * \right], \f] |
---|
183 | * |
---|
184 | * \f[ \tilde Q(\lambda) = Q - \left( \begin{array}{ccc} |
---|
185 | * {\lambda_1} & & 0 \\ & \ddots & \\ 0 & & \lambda_n \end{array} |
---|
186 | * \right). \f] |
---|
187 | * |
---|
188 | * The convex relaxation of the initial constraint is then given by |
---|
189 | * the two constraints |
---|
190 | * |
---|
191 | * \f[ \eta \geq \tilde a_0(\lambda_{\min}) + \tilde |
---|
192 | * a(\lambda_{\min})^T x + x^T \tilde Q(\lambda_{\min}) x \f] |
---|
193 | * |
---|
194 | * \f[ \eta \leq \tilde a_0(\lambda_{\max}) + \tilde |
---|
195 | * a(\lambda_{\max})^T x + x^T \tilde Q(\lambda_{\max}) x \f] |
---|
196 | * |
---|
197 | * The cut is computed as follow. Let \f$ (x^*, \eta^*) \f$ be |
---|
198 | * the solution at hand. The two outer-approximation cuts are: |
---|
199 | * |
---|
200 | * \f[ \eta \geq \tilde a_0(\lambda_{\min}) + \tilde |
---|
201 | * a(\lambda_{\min})^T x + {x^*}^T \tilde Q(\lambda_{\min}) (2x - |
---|
202 | * x^*) \f] |
---|
203 | * |
---|
204 | * and |
---|
205 | * |
---|
206 | * \f[ \eta \leq \tilde a_0(\lambda_{\max}) + \tilde |
---|
207 | * a(\lambda_{\max})^T x + {x^*}^T \tilde Q(\lambda_{\max}) (2x - |
---|
208 | * x^*); \f] |
---|
209 | * |
---|
210 | * grouping coefficients, we get: |
---|
211 | * |
---|
212 | * \f[ {x^*}^T \tilde Q(\lambda_{\min}) x^* - \tilde |
---|
213 | * a_0(\lambda_{\min}) \geq (\tilde a(\lambda_{\min}) + 2 |
---|
214 | * \tilde Q(\lambda_{\min} ) x^*)^T x - \eta \f] |
---|
215 | * |
---|
216 | * and |
---|
217 | * |
---|
218 | * \f[ {x^*}^T \tilde Q(\lambda_{\max}) x^* - \tilde |
---|
219 | * a_0(\lambda_{\max}) \leq (\tilde a(\lambda_{\max}) + 2 |
---|
220 | * \tilde Q (\lambda_{\max}) x^* )^T x - \eta \f] |
---|
221 | */ |
---|
222 | |
---|
223 | void quadCuts (expression *w, OsiCuts & cs, const CouenneCutGenerator * cg); |
---|
224 | |
---|
225 | /// Compare two exprQuad |
---|
226 | virtual int compare (exprQuad &); |
---|
227 | |
---|
228 | /// Code for comparisons |
---|
229 | virtual enum expr_type code () |
---|
230 | {return COU_EXPRQUAD;} |
---|
231 | |
---|
232 | /// Used in rank-based branching variable choice |
---|
233 | virtual int rank (); |
---|
234 | |
---|
235 | /// is this expression integer? |
---|
236 | virtual bool isInteger (); |
---|
237 | |
---|
238 | /// fill in the set with all indices of variables appearing in the |
---|
239 | /// expression |
---|
240 | virtual int DepList (std::set <int> &deplist, |
---|
241 | enum dig_type type = ORIG_ONLY); |
---|
242 | |
---|
243 | /// Set up branching object by evaluating many branching points for |
---|
244 | /// each expression's arguments |
---|
245 | virtual CouNumber selectBranch (const CouenneObject *obj, |
---|
246 | const OsiBranchingInformation *info, |
---|
247 | expression * &var, |
---|
248 | double * &brpts, |
---|
249 | double * &brDist, // distance of current LP |
---|
250 | // point to new convexifications |
---|
251 | int &way); |
---|
252 | |
---|
253 | /// Fill dependence set of the expression associated with this |
---|
254 | /// auxiliary variable |
---|
255 | virtual void fillDepSet (std::set <DepNode *, compNode> *dep, DepGraph *g); |
---|
256 | |
---|
257 | /// replace variable x with new (aux) w |
---|
258 | virtual void replace (exprVar *x, exprVar *w); |
---|
259 | |
---|
260 | /// replace variable x with new (aux) w |
---|
261 | virtual void realign (const CouenneProblem *p); |
---|
262 | |
---|
263 | /// implied bound processing |
---|
264 | virtual bool impliedBound (int, CouNumber *, CouNumber *, t_chg_bounds *); |
---|
265 | |
---|
266 | /// method to compute the bound based on sign: -1 for lower, +1 for |
---|
267 | /// upper |
---|
268 | CouNumber computeQBound (int sign); |
---|
269 | |
---|
270 | /// compute $y^{lv}$ and $y^{uv}$ for Violation Transfer algorithm |
---|
271 | virtual void closestFeasible (expression *varind, |
---|
272 | expression *vardep, |
---|
273 | CouNumber &left, |
---|
274 | CouNumber &right) const; |
---|
275 | protected: |
---|
276 | |
---|
277 | /// return lower and upper bound of quadratic expression |
---|
278 | void computeQuadFiniteBound (CouNumber &qMin, CouNumber &qMax, |
---|
279 | CouNumber *l, CouNumber *u, |
---|
280 | int &indInfLo, int &indInfUp); |
---|
281 | |
---|
282 | /// can this expression be further linearized or are we on its |
---|
283 | /// concave ("bad") side |
---|
284 | virtual bool isCuttable (CouenneProblem *problem, int index) const |
---|
285 | {return false;} // !!! TODO: specialize |
---|
286 | }; |
---|
287 | |
---|
288 | |
---|
289 | /// Compute sum of linear and nonlinear terms |
---|
290 | |
---|
291 | inline CouNumber exprQuad::operator () () { |
---|
292 | |
---|
293 | // compute non-quadratic part (linear+constant) |
---|
294 | CouNumber ret = exprGroup::operator () (); |
---|
295 | |
---|
296 | for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) { |
---|
297 | |
---|
298 | int xind = row -> first -> Index (); |
---|
299 | CouNumber x = (*(row -> first)) (); |
---|
300 | |
---|
301 | for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) { |
---|
302 | |
---|
303 | CouNumber term = x * (*(col -> first)) () * col -> second; |
---|
304 | ret += (col -> first -> Index () == xind) ? term : 2. * term; |
---|
305 | } |
---|
306 | } |
---|
307 | |
---|
308 | return (CouNumber) ret; |
---|
309 | } |
---|
310 | |
---|
311 | #endif |
---|