exprBQuad.cpp

Go to the documentation of this file.

/* $Id: exprBQuad.cpp 490 2011-01-14 16:07:12Z pbelotti $ */
/*
 * Name:    exprBQuad.cpp
 * Author:  Pietro Belotti
 * Purpose: method to compute value of an expr?BQuad
 *
 * (C) Carnegie-Mellon University, 2006. 
 * This file is licensed under the Eclipse Public License (EPL)
 */

#include "CoinHelperFunctions.hpp"
#include "CouenneExprBQuad.hpp"

using namespace Couenne;

//#define DEBUG

CouNumber exprQuad::computeQBound (int sign) {

  // 1) loose: disaggregated bound
  // 2) tighter: aggregate coefficient per variable
  // 3) tightest: solve (convex) QP on available alpha-convexification

  // 1) loose bound
  //
  // w = a0 + a'x + x'Qx means that its lower bound is
  //
  // w_l = a0 + a'_+ x_l + a'_- x_u + \sum_


  // 2) tighter bound -- TODO
  //
  // Compute lower (if sign == -1) or upper (sign == +1) bound of an
  // exprQuad based on the information obtained through
  // alpha-convexification, if any, or as follows:
  //
  // w = a0 + a'x + x'Qx = 
  //   = a0 + sum{i} [(a_i + sum {j>=i} q_ij * x_j) * x_i] =
  //   = a0 + sum{i} [                          z_i * x_i]
  // 
  // Thus, some bound on z_i can be computed and a bound on the whole
  // expression should be better than what can be obtained by summing
  // all bounds separately.
  //
  // Notice that the above computation is fast and may be better than
  // the convexification after some updates in the variable bounds
  // without updating the convexification. Notice also that the
  // direction can also be vertical, not only horizontal.

  CouNumber bound = c0_, term;

  // derive linear part (obtain constant)
  for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el) {

    CouNumber 
      coe = el -> second, term = 0.,
      lb  = el -> first -> lb (),
      ub  = el -> first -> ub ();

    if (((coe < 0.) && (sign < 0)) || 
        ((coe > 0.) && (sign > 0))) 
      {if    ((term=ub) >  COUENNE_INFINITY) return (sign < 0)? -COUENNE_INFINITY : COUENNE_INFINITY;}
    else {if ((term=lb) < -COUENNE_INFINITY) return (sign < 0)? -COUENNE_INFINITY : COUENNE_INFINITY;}

    bound += coe * term;
  }

#ifdef DEBUG
  printf ("quadBound --- linear, %cb = %g\n", (sign < 0) ? 'l' : 'u', bound);
#endif

  // derive quadratic part (obtain linear part)
  for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {

    int xind = row -> first -> Index ();

      CouNumber 
        lbi = row -> first -> lb (),
        ubi = row -> first -> ub ();

    for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {

      int yind = col -> first -> Index ();

      CouNumber coe = col -> second;

      if (xind == yind) { // term of the form q_ii x_i^2

        if (((coe > 0.) && (sign < 0)) ||
            ((coe < 0.) && (sign > 0))) 
          term = (ubi < 0) ? (ubi * ubi) : (lbi > 0) ? (lbi * lbi) : 0.; //min{xi^2: xi in [lbi,ubi]
        else 
          if ((term = CoinMax (lbi*lbi, ubi*ubi)) > COUENNE_INFINITY) 
            return (sign < 0) ? -COUENNE_INFINITY : COUENNE_INFINITY;

        term *= coe;

#ifdef DEBUG
        printf ("Qii %d %g %g -> %g\n", xind, coe, term, bound + term);
#endif
      } else {

        coe *= 2;

        CouNumber
          lbj = col -> first -> lb (),
          ubj = col -> first -> ub (),
          b1 = coe * lbi * lbj, 
          b2 = coe * lbi * ubj,
          b3 = coe * ubi * lbj, 
          b4 = coe * ubi * ubj;

        if (fabs (lbi) == 0) b1 = b2 = 0;
        if (fabs (lbj) == 0) b1 = b3 = 0;
        if (fabs (ubi) == 0) b3 = b4 = 0;
        if (fabs (ubj) == 0) b2 = b4 = 0;

        if (sign < 0) {
          if ((term = CoinMin (CoinMin (b1, b2), CoinMin (b3, b4))) < -COUENNE_INFINITY) 
            return -COUENNE_INFINITY; 
        } else
          if ((term = CoinMax (CoinMax (b1, b2), CoinMax (b3, b4))) >  COUENNE_INFINITY)
            return  COUENNE_INFINITY;

#ifdef DEBUG
        printf ("Qij %d %d %g %g -> %g\n", xind, yind, coe, term, bound + term);
#endif
      }

      //      if ((i!=j) || (lbi >= 0) || (ubi <= 0))
      bound += term;
    }
  }

  return bound;
}