source: stable/0.2/Couenne/src/expression/operators/exprQuad.hpp @ 159

Last change on this file since 159 was 159, checked in by pbelotti, 11 years ago

created new stable branch 0.2 from trunk (rev. 157)

File size: 10.2 KB
Line 
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
21class quadElem;
22class CouenneProblem;
23class 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
42class exprQuad: public exprGroup {
43
44public:
45
46  /// matrix
47  typedef std::vector <std::pair <exprVar *, CouNumber> >  sparseQcol;
48  typedef std::vector <std::pair <exprVar *, sparseQcol> > sparseQ;
49
50protected:
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
79public:
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;
275protected:
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
291inline 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
Note: See TracBrowser for help on using the repository browser.