ROL_MoreauYosidaObjective.hpp

Go to the documentation of this file.

// @HEADER
// ************************************************************************
//
//               Rapid Optimization Library (ROL) Package
//                 Copyright (2014) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact lead developers:
//              Drew Kouri   (dpkouri@sandia.gov) and
//              Denis Ridzal (dridzal@sandia.gov)
//
// ************************************************************************
// @HEADER

#ifndef ROL_MOREAUYOSIDAOBJECTIVE_H
#define ROL_MOREAUYOSIDAOBJECTIVE_H

#include "ROL_Objective.hpp"
#include "ROL_BoundConstraint.hpp"
#include "ROL_Vector.hpp"
#include "ROL_Types.hpp"
#include "ROL_Ptr.hpp"
#include "ROL_ScalarController.hpp"
#include "ROL_ParameterList.hpp"
#include <iostream>

namespace ROL {

template <class Real>
class MoreauYosidaObjective : public Objective<Real> {
private:
  const Ptr<Objective<Real>> obj_;
  const Ptr<BoundConstraint<Real>> bnd_;

  Ptr<Vector<Real>> l_;
  Ptr<Vector<Real>> u_;
  Ptr<Vector<Real>> l1_;
  Ptr<Vector<Real>> u1_;
  Ptr<Vector<Real>> dl1_;
  Ptr<Vector<Real>> du1_;
  Ptr<Vector<Real>> xlam_;
  Ptr<Vector<Real>> v_;
  Ptr<Vector<Real>> dv_;
  Ptr<Vector<Real>> dv2_;
  Ptr<Vector<Real>> lam_;
  Ptr<Vector<Real>> tmp_;

  Ptr<ScalarController<Real,int>> fval_;
  Ptr<VectorController<Real,int>> gradient_;

  Real mu_;
  bool isPenEvaluated_;
  int nfval_;
  int ngrad_;
  bool updateMultiplier_;
  bool updatePenalty_;

  void computePenalty(const Vector<Real> &x) {
    if ( bnd_->isActivated() ) {
      Real one = 1.0;
      if ( !isPenEvaluated_ ) {
        xlam_->set(x);
        xlam_->axpy(one/mu_,*lam_);

        if ( bnd_->isFeasible(*xlam_) ) {
          l1_->zero(); dl1_->zero();
          u1_->zero(); du1_->zero();
        }
        else {
          // Compute lower penalty component
          l1_->set(*l_);
          bnd_->pruneLowerInactive(*l1_,*xlam_);
          tmp_->set(*xlam_);
          bnd_->pruneLowerInactive(*tmp_,*xlam_);
          l1_->axpy(-one,*tmp_);

          // Compute upper penalty component
          u1_->set(*xlam_);
          bnd_->pruneUpperInactive(*u1_,*xlam_);
          tmp_->set(*u_);
          bnd_->pruneUpperInactive(*tmp_,*xlam_);
          u1_->axpy(-one,*tmp_);

          // Compute derivative of lower penalty component
          dl1_->set(l1_->dual());
          bnd_->pruneLowerInactive(*dl1_,*xlam_);

          // Compute derivative of upper penalty component
          du1_->set(u1_->dual());
          bnd_->pruneUpperInactive(*du1_,*xlam_);
        }

        isPenEvaluated_ = true;
      }
    }
  }

  void initialize(const Vector<Real> &x,
                  const Vector<Real> &g) {
    fval_     = makePtr<ScalarController<Real,int>>();
    gradient_ = makePtr<VectorController<Real,int>>();

    l_    = x.clone();
    l1_   = x.clone();
    dl1_  = g.clone();
    u_    = x.clone();
    u1_   = x.clone();
    du1_  = g.clone();
    xlam_ = x.clone();
    v_    = x.clone();
    dv_   = g.clone();
    dv2_  = g.clone();
    lam_  = x.clone();
    tmp_  = x.clone();

    l_->set(*bnd_->getLowerBound());
    u_->set(*bnd_->getUpperBound());

    lam_->zero();
    //lam_->set(*u_);
    //lam_->plus(*l_);
    //lam_->scale(0.5);
  }

public:
  MoreauYosidaObjective(const Ptr<Objective<Real>> &obj,
                      const Ptr<BoundConstraint<Real>> &bnd, 
                      const Vector<Real> &x,
                      const Vector<Real> &g,
                      const Real mu = 1e1,
                      const bool updateMultiplier = true,
                      const bool updatePenalty = true)
    : obj_(obj), bnd_(bnd), mu_(mu),
      isPenEvaluated_(false), nfval_(0), ngrad_(0),
      updateMultiplier_(updateMultiplier), updatePenalty_(updatePenalty) {
    initialize(x,g);
  }

  MoreauYosidaObjective(const Ptr<Objective<Real>> &obj,
                      const Ptr<BoundConstraint<Real>> &bnd, 
                      const Vector<Real> &x,
                      const Vector<Real> &g,
                      const Vector<Real> &lam,
                      const Real mu = 1e1,
                      const bool updateMultiplier = true,
                      const bool updatePenalty = true)
    : MoreauYosidaObjective(obj,bnd,x,g,mu,updateMultiplier,updatePenalty) {
    lam_->set(lam);
  }

  MoreauYosidaObjective(const Ptr<Objective<Real>> &obj,
                      const Ptr<BoundConstraint<Real>> &bnd,
                      const Vector<Real> &x,
                      const Vector<Real> &g,
                      ParameterList &parlist)
    : obj_(obj), bnd_(bnd),
      isPenEvaluated_(false), nfval_(0), ngrad_(0) {
    initialize(x,g);
    ParameterList &list = parlist.sublist("Step").sublist("Moreau-Yosida Penalty");
    updateMultiplier_ = list.get("Update Multiplier",true);
    updatePenalty_    = list.get("Update Penalty",true);
    mu_               = list.get("Initial Penalty Parameter",1e1);
  }

  MoreauYosidaObjective(const Ptr<Objective<Real>> &obj,
                      const Ptr<BoundConstraint<Real>> &bnd,
                      const Vector<Real> &x,
                      const Vector<Real> &g,
                      const Vector<Real> &lam,
                      ParameterList &parlist)
    : MoreauYosidaObjective(obj,bnd,x,g,parlist) {
    lam_->set(lam);
  }

  void updateMultipliers(Real mu, const Vector<Real> &x) {
    if ( bnd_->isActivated() ) {
      if ( updateMultiplier_ ) {
        const Real one(1);
        computePenalty(x);
        lam_->set(*u1_);
        lam_->axpy(-one,*l1_);
        lam_->scale(mu_);
      }
      if ( updatePenalty_ ) {
        mu_ = mu;
      }
    }
    nfval_ = 0; ngrad_ = 0;
    isPenEvaluated_ = false;
  }

  void reset(const Real mu) {
    lam_->zero();
    mu_ = mu;
    nfval_ = 0; ngrad_ = 0;
    isPenEvaluated_ = false;
  }

  Real testComplementarity(const Vector<Real> &x) {
    Real val(0);
    if (bnd_->isActivated()) {
      computePenalty(x);

      tmp_->set(x);
      tmp_->axpy(static_cast<Real>(-1), *l_);
      Real lower = mu_*std::abs(tmp_->dot(*l1_));

      tmp_->set(x);
      tmp_->axpy(static_cast<Real>(-1), *u_);
      Real upper = mu_*std::abs(tmp_->dot(*u1_));

      tmp_->set(x);
      bnd_->project(*tmp_);
      tmp_->axpy(static_cast<Real>(-1), x);
      Real xnorm = tmp_->norm();

      val = std::max(xnorm,std::max(lower,upper));
    }
    return val;
  }

  Real getObjectiveValue(const Vector<Real> &x, Real &tol) {
    int key(0);
    Real val(0);
    bool isComputed = fval_->get(val,key);
    if (!isComputed) {
      val = obj_->value(x,tol); nfval_++;
      fval_->set(val,key);
    }
    return val;
  }

  void getObjectiveGradient(Vector<Real> &g, const Vector<Real> &x, Real &tol) {
    int key(0);
    bool isComputed = gradient_->get(g,key);
    if (!isComputed) {
      obj_->gradient(g,x,tol); ngrad_++;
      gradient_->set(g,key);
    }
  }

  int getNumberFunctionEvaluations(void) {
    return nfval_;
  }

  int getNumberGradientEvaluations(void) {
    return ngrad_;
  }

  void update( const Vector<Real> &x, UpdateType type, int iter = -1 ) {
    obj_->update(x,type,iter);
    fval_->objectiveUpdate(type);
    gradient_->objectiveUpdate(type);
    // Need to do something smart here
    isPenEvaluated_  = false;
  }

  Real value( const Vector<Real> &x, Real &tol ) {
    // Compute objective function value
    Real fval = getObjectiveValue(x,tol);
    // Add value of the Moreau-Yosida penalty
    if ( bnd_->isActivated() ) {
      const Real half(0.5);
      computePenalty(x);
      fval += half*mu_*(l1_->dot(*l1_) + u1_->dot(*u1_));
    }
    return fval;
  }

  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
    // Compute gradient of objective function
    getObjectiveGradient(g,x,tol);
    // Add gradient of the Moreau-Yosida penalty
    if ( bnd_->isActivated() ) {
      computePenalty(x);
      g.axpy(-mu_,*dl1_);
      g.axpy(mu_,*du1_);
    }
  }

  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
    // Apply objective Hessian to a vector
    obj_->hessVec(hv,v,x,tol);
    // Add Hessian of the Moreau-Yosida penalty
    if ( bnd_->isActivated() ) {
      const Real one(1);
      computePenalty(x);

      v_->set(v);
      bnd_->pruneLowerActive(*v_,*xlam_);
      v_->scale(-one);
      v_->plus(v);
      dv_->set(v_->dual());
      dv2_->set(*dv_);
      bnd_->pruneLowerActive(*dv_,*xlam_);
      dv_->scale(-one);
      dv_->plus(*dv2_);
      hv.axpy(mu_,*dv_);

      v_->set(v);
      bnd_->pruneUpperActive(*v_,*xlam_);
      v_->scale(-one);
      v_->plus(v);
      dv_->set(v_->dual());
      dv2_->set(*dv_);
      bnd_->pruneUpperActive(*dv_,*xlam_);
      dv_->scale(-one);
      dv_->plus(*dv2_);
      hv.axpy(mu_,*dv_);
    }
  }
}; // class MoreauYosidaObjective

} // namespace ROL

#endif