50 #include "Teuchos_RCP.hpp" 81 void remove(
const std::vector<unsigned> &ind) {
82 for (
unsigned j = ind.back()+1; j <
size_; j++) {
83 (subgradients_[j-1])->
set(*(subgradients_[j]));
84 linearizationErrors_[j-1] = linearizationErrors_[j];
85 distanceMeasures_[j-1] = distanceMeasures_[j];
86 dualVariables_[j-1] = dualVariables_[j];
88 (subgradients_[size_-1])->zero();
91 dualVariables_[size_-1] = 0.0;
92 for (
unsigned i = ind.size()-1; i > 0; --i) {
93 for (
unsigned j = ind[i-1]+1; j <
size_; j++) {
94 (subgradients_[j-1])->
set(*(subgradients_[j]));
95 linearizationErrors_[j-1] = linearizationErrors_[j];
96 distanceMeasures_[j-1] = distanceMeasures_[j];
97 dualVariables_[j-1] = dualVariables_[j];
104 (subgradients_[
size_])->
set(g);
105 linearizationErrors_[
size_] = le;
106 distanceMeasures_[
size_] = dm;
107 dualVariables_[
size_] = 0.0;
112 Teuchos::RCP<Vector<Real> >
tG_;
113 Teuchos::RCP<Vector<Real> >
eG_;
114 Teuchos::RCP<Vector<Real> >
yG_;
121 Bundle(
const unsigned maxSize = 10,
const Real coeff = 0.0,
const unsigned remSize = 2)
122 : size_(0), maxSize_(maxSize), remSize_(remSize), coeff_(coeff), isInitialized_(false) {
123 remSize_ = ((remSize_ < 2) ? 2 : ((remSize_ > maxSize_-1) ? maxSize_-1 : remSize_));
124 subgradients_.clear();
125 subgradients_.assign(maxSize,Teuchos::null);
126 linearizationErrors_.clear();
128 distanceMeasures_.clear();
130 dualVariables_.clear();
131 dualVariables_.assign(maxSize_,0.0);
135 if ( !isInitialized_ ) {
136 for (
unsigned i = 0; i <
maxSize_; i++) {
137 subgradients_[i] = g.
clone();
139 (subgradients_[0])->
set(g);
140 linearizationErrors_[0] = 0.0;
141 distanceMeasures_[0] = 0.0;
142 dualVariables_[0] = 1.0;
144 isInitialized_ =
true;
154 return linearizationErrors_[i];
158 return distanceMeasures_[i];
162 return *(subgradients_[i]);
168 alpha = std::max(coeff_*std::pow(dm,2.0),le);
173 const Real
alpha(
const unsigned i)
const {
174 return computeAlpha(distanceMeasures_[i],linearizationErrors_[i]);
182 aggSubGrad.
zero(); aggLinErr = 0.0; aggDistMeas = 0.0; eG_->zero();
183 Real eLE = 0.0, eDM = 0.0, yLE = 0.0, yDM = 0.0, tLE = 0.0, tDM = 0.0;
184 for (
unsigned i = 0; i <
size_; i++) {
186 tG_->set(aggSubGrad);
187 yG_->set(*subgradients_[i]); yG_->scale(dualVariables_[i]); yG_->plus(*eG_);
188 aggSubGrad.
set(*tG_); aggSubGrad.
plus(*yG_);
189 eG_->set(*tG_); eG_->axpy(-1.0,aggSubGrad); eG_->plus(*yG_);
192 yLE = dualVariables_[i]*linearizationErrors_[i] + eLE;
193 aggLinErr = tLE + yLE;
194 eLE = (tLE - aggLinErr) + yLE;
197 yDM = dualVariables_[i]*distanceMeasures_[i] + eDM;
198 aggDistMeas = tDM + yDM;
199 eDM = (tDM - aggDistMeas) + yDM;
204 if (size_ == maxSize_) {
206 unsigned loc =
size_, cnt = 0;
207 std::vector<unsigned> ind(remSize_,0);
208 for (
unsigned i = size_; i > 0; --i) {
209 if ( std::abs(linearizationErrors_[i-1]) <
ROL_EPSILON ) {
214 for (
unsigned i = 0; i <
size_; i++) {
215 if ( i < loc || i > loc ) {
219 if (cnt == remSize_) {
230 void update(
const bool flag,
const Real linErr,
const Real distMeas,
234 for (
unsigned i = 0; i <
size_; i++) {
235 linearizationErrors_[i] += linErr - subgradients_[i]->dot(s.
dual());
236 distanceMeasures_[i] += distMeas;
238 linearizationErrors_[
size_] = 0.0;
239 distanceMeasures_[
size_] = 0.0;
243 linearizationErrors_[
size_] = linErr;
244 distanceMeasures_[
size_] = distMeas;
247 (subgradients_[
size_])->
set(g);
249 dualVariables_[
size_] = 0.0;
258 return dualVariables_[i];
262 dualVariables_[i] = val;
266 dualVariables_.assign(size_,0.0);
273 Teuchos::RCP<Vector<Real> >
gx_;
274 Teuchos::RCP<Vector<Real> >
ge_;
293 Real sum = 0.0, err = 0.0, tmp = 0.0, y = 0.0;
294 for (
unsigned i = 0; i <
size_; i++) {
297 y = dualVariables_[i] + err;
299 err = (tmp - sum) + y;
301 for (
unsigned i = 0; i <
size_; i++) {
302 dualVariables_[i] /= sum;
304 nworkingSet_.clear();
306 for (
unsigned i = 0; i <
size_; i++) {
307 if ( dualVariables_[i] == 0.0 ) {
308 workingSet_.insert(i);
311 nworkingSet_.insert(i);
317 gx_->zero(); eG_->zero();
318 for (
unsigned i = 0; i <
size_; i++) {
321 yG_->set(*eG_); yG_->axpy(x[i],*(subgradients_[i]));
322 gx_->set(*tG_); gx_->plus(*yG_);
323 eG_->set(*tG_); eG_->axpy(-1.0,*gx_); eG_->plus(*yG_);
325 Real Hx = 0.0, val = 0.0, err = 0.0, tmp = 0.0, y = 0.0;
326 for (
unsigned i = 0; i <
size_; i++) {
328 Hx = gx_->dot(*(subgradients_[i]));
331 y = x[i]*(0.5*Hx +
alpha(i)/t) + err;
333 err = (tmp - val) + y;
335 g[i] = Hx +
alpha(i)/t;
341 gx_->zero(); eG_->zero();
342 for (
unsigned i = 0; i <
size_; i++) {
345 yG_->set(*eG_); yG_->axpy(x[i],*(subgradients_[i]));
346 gx_->set(*tG_); gx_->plus(*yG_);
347 eG_->set(*tG_); eG_->axpy(-1.0,*gx_); eG_->plus(*yG_);
349 for (
unsigned i = 0; i <
size_; i++) {
351 Hx[i] = subgradients_[i]->dot(*gx_);
355 void applyMatrix(std::vector<Real> &Hx,
const std::vector<Real> &x)
const {
356 gx_->zero(); eG_->zero();
357 unsigned n = nworkingSet_.size();
358 typename std::set<unsigned>::iterator it = nworkingSet_.begin();
359 for (
unsigned i = 0; i < n; i++) {
362 yG_->set(*eG_); yG_->axpy(x[i],*(subgradients_[*it]));
363 gx_->set(*tG_); gx_->plus(*yG_);
364 eG_->set(*tG_); eG_->axpy(-1.0,*gx_); eG_->plus(*yG_);
367 it = nworkingSet_.begin();
368 for (
unsigned i = 0; i < n; i++) {
370 Hx[i] = subgradients_[*it]->dot(*gx_); it++;
374 void computeLagMult(std::vector<Real> &lam,
const Real mu,
const std::vector<Real> &g)
const {
375 unsigned n = workingSet_.size();
378 typename std::set<unsigned>::iterator it = workingSet_.begin();
379 for (
unsigned i = 0; i < n; i++) {
380 lam[i] = g[*it] - mu; it++;
390 unsigned n = workingSet_.size(); ind =
size_;
393 typename std::set<unsigned>::iterator it = workingSet_.begin();
394 for (
unsigned i = 0; i < n; i++) {
406 Real
computeAlpha(
unsigned &ind,
const std::vector<Real> &x,
const std::vector<Real> &p)
const {
408 typename std::set<unsigned>::iterator it;
409 for (it = nworkingSet_.begin(); it != nworkingSet_.end(); it++) {
411 tmp = -x[*it]/p[*it];
412 if ( alpha >= tmp ) {
418 return std::max(0.0,alpha);
422 const std::vector<Real> &g,
const Real tol)
const {
424 unsigned n = nworkingSet_.size(), cnt = 0;
428 std::vector<Real> gk(n,0.0);
429 typename std::set<unsigned>::iterator it = nworkingSet_.begin();
430 for (
unsigned i = 0; i < n; i++) {
431 gk[i] = g[*it]; it++;
433 std::vector<Real> sk(n,0.0);
435 it = nworkingSet_.begin();
436 for (
unsigned i = 0; i < n; i++) {
437 s[*it] = sk[i]; it++;
445 std::vector<Real> tmp(Px.size(),0.0);
454 void applyG(std::vector<Real> &Gx,
const std::vector<Real> &x)
const {
464 unsigned dim = nworkingSet_.size();
465 Real sum = 0.0, err = 0.0, tmp = 0.0, y = 0.0;
466 for (
unsigned i = 0; i < dim; i++) {
471 err = (tmp - sum) + y;
474 for (
unsigned i = 0; i < dim; i++) {
480 Gx.assign(x.begin(),x.end());
484 unsigned dim = nworkingSet_.size();
485 Real eHe = 0.0, sum = 0.0;
486 Real errX = 0.0, tmpX = 0.0, yX = 0.0, errE = 0.0, tmpE = 0.0, yE = 0.0;
487 std::vector<Real> gg(dim,0.0);
488 typename std::set<unsigned>::iterator it = nworkingSet_.begin();
489 for (
unsigned i = 0; i < dim; i++) {
490 gg[i] = 1.0/std::abs(subgradients_[*it]->
dot(*(subgradients_[*it]))); it++;
493 yX = x[i]*gg[i] + errX;
495 errX = (tmpX - sum) + yX;
500 errE = (tmpE - eHe) + yE;
503 for (
unsigned i = 0; i < dim; i++) {
504 Px[i] = (x[i]-sum)*gg[i];
508 void applyG_Jacobi(std::vector<Real> &Gx,
const std::vector<Real> &x)
const {
509 unsigned dim = nworkingSet_.size();
510 typename std::set<unsigned>::iterator it = nworkingSet_.begin();
511 for (
unsigned i = 0; i < dim; i++) {
512 Gx[i] = std::abs(subgradients_[*it]->
dot(*(subgradients_[*it])))*x[i]; it++;
517 int dim = nworkingSet_.size();
519 gx_->zero(); ge_->zero();
520 Real eHx = 0.0, eHe = 0.0;
522 std::vector<Real> x1(dim,0.0), e1(dim,0.0),gg(dim,0.0);
523 typename std::set<unsigned>::iterator it, jt;
524 it = nworkingSet_.begin();
525 for (
int i = 0; i < dim; i++) {
526 gx_->zero(); ge_->zero(); jt = nworkingSet_.begin();
527 for (
int j = 0; j < i; j++) {
528 gx_->axpy(x1[j],*(subgradients_[*jt]));
529 ge_->axpy(e1[j],*(subgradients_[*jt]));
532 gg[i] = subgradients_[*it]->dot(*(subgradients_[*it]));
533 x1[i] = (x[i] - gx_->dot(*(subgradients_[*it])))/gg[i];
534 e1[i] = (1.0 - ge_->dot(*(subgradients_[*it])))/gg[i];
538 for (
int i = 0; i < dim; i++) {
543 std::vector<Real> Hx(dim,0.0), He(dim,0.0); it = nworkingSet_.end();
544 for (
int i = dim-1; i >= 0; --i) {
546 gx_->zero(); ge_->zero(); jt = nworkingSet_.end();
547 for (
int j = dim-1; j >= i+1; --j) {
549 gx_->axpy(Hx[j],*(subgradients_[*jt]));
550 ge_->axpy(He[j],*(subgradients_[*jt]));
552 Hx[i] = (x1[i] - gx_->dot(*(subgradients_[*it])))/gg[i];
553 He[i] = (e1[i] - ge_->dot(*(subgradients_[*it])))/gg[i];
559 for (
int i = 0; i < dim; i++) {
560 Px[i] = Hx[i] - (eHx/eHe)*He[i];
564 void applyG_SymGS(std::vector<Real> &Gx,
const std::vector<Real> &x)
const {
565 unsigned dim = nworkingSet_.size();
566 typename std::set<unsigned>::iterator it = nworkingSet_.begin();
567 for (
unsigned i = 0; i < dim; i++) {
568 Gx[i] = std::abs(subgradients_[*it]->
dot(*(subgradients_[*it])))*x[i]; it++;
573 unsigned n = g.size();
574 std::vector<Real> Gg(n,0.0);
575 Real y = 0.0, ytmp = 0.0, yprt = 0.0, yerr = 0.0;
579 for (
unsigned i = 0; i < n; i++) {
581 yprt = (r[i] - Gg[i]) + yerr;
583 yerr = (ytmp - y) + yprt;
586 for (
unsigned i = 0; i < n; i++) {
592 unsigned projectedCG(std::vector<Real> &x, Real &mu,
const std::vector<Real> &b,
const Real tol)
const {
593 unsigned n = nworkingSet_.size();
594 std::vector<Real> r(n,0.0), r0(n,0.0), g(n,0.0), d(n,0.0), Ad(n,0.0);
598 r0.assign(r.begin(),r.end());
601 Real rg =
dot(r,g), rg0 = 0.0;
604 Real
alpha = 0.0, kappa = 0.0, beta = 0.0;
605 Real CGtol = std::min(tol,1.e-2*rg);
607 while (rg > CGtol && cnt < 2*n+1) {
624 Real err = 0.0, tmp = 0.0, y = 0.0;
625 for (
unsigned i = 0; i < n; i++) {
629 err = (tmp - mu) + y;
636 Real
dot(
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
638 Real val = 0.0, err = 0.0, tmp = 0.0, y0 = 0.0;
639 unsigned n = std::min(x.size(),y.size());
640 for (
unsigned i = 0; i < n; i++) {
642 y0 = x[i]*y[i] + err;
644 err = (tmp - val) + y0;
649 Real
norm(
const std::vector<Real> &x)
const {
650 return std::sqrt(
dot(x,x));
653 void axpy(
const Real a,
const std::vector<Real> &x, std::vector<Real> &y)
const {
654 unsigned n = std::min(y.size(),x.size());
655 for (
unsigned i = 0; i < n; i++) {
660 void scale(std::vector<Real> &x,
const Real a)
const {
661 for (
unsigned i = 0; i < x.size(); i++) {
666 void scale(std::vector<Real> &x,
const Real a,
const std::vector<Real> &y)
const {
667 unsigned n = std::min(x.size(),y.size());
668 for (
unsigned i = 0; i < n; i++) {
673 unsigned solveDual_dim1(
const Real t,
const unsigned maxit = 1000,
const Real tol = 1.e-8) {
674 dualVariables_[0] = 1.0;
679 unsigned solveDual_dim2(
const Real t,
const unsigned maxit = 1000,
const Real tol = 1.e-8) {
680 gx_->set(*subgradients_[0]); gx_->axpy(-1.0,*subgradients_[1]);
681 Real diffg = gx_->dot(*gx_);
684 Real gdiffg = subgradients_[1]->dot(*gx_);
685 dualVariables_[0] = std::min(1.0,std::max(0.0,-(gdiffg+diffa)/diffg));
686 dualVariables_[1] = 1.0-dualVariables_[0];
691 dualVariables_[0] = 1.0; dualVariables_[1] = 0.0;
694 dualVariables_[0] = 0.0; dualVariables_[1] = 1.0;
698 dualVariables_[0] = 0.5; dualVariables_[1] = 0.5;
708 unsigned ind = 0, i = 0, CGiter = 0;
709 Real snorm = 0.0,
alpha = 0.0, mu = 0.0;
710 std::vector<Real> s(size_,0.0), Hs(size_,0.0), g(size_,0.0), lam(size_+1,0.0);
713 for (i = 0; i < maxit; i++) {
725 workingSet_.erase(ind);
726 nworkingSet_.insert(ind);
738 workingSet_.insert(ind);
739 nworkingSet_.erase(ind);
751 virtual unsigned solveDual(
const Real t,
const unsigned maxit = 1000,
const Real tol = 1.e-8) {
756 else if (size_ == 2) {
769 void project(std::vector<Real> &x,
const std::vector<Real> &v)
const {
770 std::vector<Real> vsort(size_,0.0);
771 vsort.assign(v.begin(),v.end());
772 std::sort(vsort.begin(),vsort.end());
773 Real sum = -1.0, lam = 0.0;
774 for (
int i = size_-1; i > 0; i--) {
776 if ( sum >= ((Real)(size_-i))*vsort[i-1] ) {
777 lam = sum/(Real)(size_-i);
782 lam = (sum+vsort[0])/(Real)
size_;
784 for (
int i = 0; i <
size_; i++) {
785 x[i] = std::max(0.0,v[i] - lam);
790 std::vector<Real> x(size_,0.0), Px(size_,0.0);
791 axpy(1.0,dualVariables_,x);
795 axpy(1.0,dualVariables_,x);
void initializeDualSolver(void)
void applyPreconditioner_Jacobi(std::vector< Real > &Px, const std::vector< Real > &x) const
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
void resetDualVariables(void)
void scale(std::vector< Real > &x, const Real a) const
std::set< unsigned > workingSet_
virtual void plus(const Vector &x)=0
Compute , where .
void update(const bool flag, const Real linErr, const Real distMeas, const Vector< Real > &g, const Vector< Real > &s)
Real norm(const std::vector< Real > &x) const
Teuchos::RCP< Vector< Real > > yG_
const Real linearizationError(const unsigned i) const
void applyFullMatrix(std::vector< Real > &Hx, const std::vector< Real > &x) const
void applyMatrix(std::vector< Real > &Hx, const std::vector< Real > &x) const
Contains definitions of custom data types in ROL.
unsigned solveEQPsubproblem(std::vector< Real > &s, Real &mu, const std::vector< Real > &g, const Real tol) const
Real computeCriticality(const std::vector< Real > &g)
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
void reset(const Vector< Real > &g, const Real le, const Real dm)
virtual void zero()
Set to zero vector.
Defines the linear algebra or vector space interface.
const Vector< Real > & subgradient(const unsigned i) const
void applyPreconditioner(std::vector< Real > &Px, const std::vector< Real > &x) const
std::set< unsigned > nworkingSet_
const Real computeAlpha(const Real dm, const Real le) const
std::vector< Real > dualVariables_
void scale(std::vector< Real > &x, const Real a, const std::vector< Real > &y) const
void applyG_SymGS(std::vector< Real > &Gx, const std::vector< Real > &x) const
Teuchos::RCP< Vector< Real > > ge_
Real dot(const std::vector< Real > &x, const std::vector< Real > &y) const
void computeLagMult(std::vector< Real > &lam, const Real mu, const std::vector< Real > &g) const
const Real alpha(const unsigned i) const
std::vector< Teuchos::RCP< Vector< Real > > > subgradients_
unsigned projectedCG(std::vector< Real > &x, Real &mu, const std::vector< Real > &b, const Real tol) const
void applyG_Identity(std::vector< Real > &Gx, const std::vector< Real > &x) const
Real computeAlpha(unsigned &ind, const std::vector< Real > &x, const std::vector< Real > &p) const
void applyG_Jacobi(std::vector< Real > &Gx, const std::vector< Real > &x) const
Real getDualVariables(const unsigned i)
std::vector< Real > distanceMeasures_
Real evaluateObjective(std::vector< Real > &g, const std::vector< Real > &x, const Real t) const
bool isNonnegative(unsigned &ind, const std::vector< Real > &x) const
void applyPreconditioner_SymGS(std::vector< Real > &Px, const std::vector< Real > &x) const
void aggregate(Vector< Real > &aggSubGrad, Real &aggLinErr, Real &aggDistMeas) const
unsigned solveDual_dim1(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
void applyPreconditioner_Identity(std::vector< Real > &Px, const std::vector< Real > &x) const
Teuchos::RCP< Vector< Real > > eG_
unsigned size(void) const
void applyG(std::vector< Real > &Gx, const std::vector< Real > &x) const
const Real distanceMeasure(const unsigned i) const
Teuchos::RCP< Vector< Real > > tG_
virtual void set(const Vector &x)
Set where .
void project(std::vector< Real > &x, const std::vector< Real > &v) const
void setDualVariables(const unsigned i, const Real val)
static const double ROL_OVERFLOW
Platform-dependent maximum double.
void axpy(const Real a, const std::vector< Real > &x, std::vector< Real > &y) const
virtual unsigned solveDual(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
std::vector< Real > linearizationErrors_
unsigned solveDual_dim2(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
void initialize(const Vector< Real > &g)
void add(const Vector< Real > &g, const Real le, const Real dm)
Teuchos::RCP< Vector< Real > > gx_
Bundle(const unsigned maxSize=10, const Real coeff=0.0, const unsigned remSize=2)
void computeResidualUpdate(std::vector< Real > &r, std::vector< Real > &g) const
Provides the interface for and implments a bundle.
static const double ROL_EPSILON
Platform-dependent machine epsilon.
unsigned solveDual_arbitrary(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)