44#ifndef ROL_BUNDLE_U_AS_DEF_H
45#define ROL_BUNDLE_U_AS_DEF_H
49template<
typename Real>
53 const unsigned remSize)
54 :
Bundle_U<Real>(maxSize,coeff,omega,remSize), isInitialized_(false) {}
56template<
typename Real>
57void Bundle_U_AS<Real>::initialize(
const Vector<Real> &g) {
58 Bundle_U<Real>::initialize(g);
59 if ( !isInitialized_ ) {
65 isInitialized_ =
true;
69template<
typename Real>
70unsigned Bundle_U_AS<Real>::solveDual(
const Real t,
const unsigned maxit,
const Real tol) {
72 if (Bundle_U<Real>::size() == 1) {
73 iter = Bundle_U<Real>::solveDual_dim1(t,maxit,tol);
75 else if (Bundle_U<Real>::size() == 2) {
76 iter = Bundle_U<Real>::solveDual_dim2(t,maxit,tol);
79 iter = solveDual_arbitrary(t,maxit,tol);
84template<
typename Real>
85void Bundle_U_AS<Real>::initializeDualSolver(
void) {
87 Real sum(0), err(0), tmp(0), y(0);
88 for (
unsigned i = 0; i < Bundle_U<Real>::size(); ++i) {
91 y = Bundle_U<Real>::getDualVariable(i) - err;
93 err = (tmp - sum) - y;
96 for (
unsigned i = 0; i < Bundle_U<Real>::size(); ++i) {
97 tmp = Bundle_U<Real>::getDualVariable(i)/sum;
98 Bundle_U<Real>::setDualVariable(i,tmp);
100 nworkingSet_.clear();
102 for (
unsigned i = 0; i < Bundle_U<Real>::size(); ++i) {
103 if ( Bundle_U<Real>::getDualVariable(i) ==
zero ) {
104 workingSet_.insert(i);
107 nworkingSet_.insert(i);
112template<
typename Real>
113void Bundle_U_AS<Real>::computeLagMult(std::vector<Real> &lam,
const Real mu,
const std::vector<Real> &g)
const {
115 unsigned n = workingSet_.size();
118 typename std::set<unsigned>::iterator it = workingSet_.begin();
119 for (
unsigned i = 0; i < n; ++i) {
120 lam[i] = g[*it] - mu; it++;
128template<
typename Real>
129bool Bundle_U_AS<Real>::isNonnegative(
unsigned &ind,
const std::vector<Real> &x)
const {
131 unsigned n = workingSet_.size();
132 ind = Bundle_U<Real>::size();
134 Real min = ROL_OVERFLOW<Real>();
135 typename std::set<unsigned>::iterator it = workingSet_.begin();
136 for (
unsigned i = 0; i < n; ++i) {
143 flag = ((min < -ROL_EPSILON<Real>()) ? false :
true);
148template<
typename Real>
149Real Bundle_U_AS<Real>::computeStepSize(
unsigned &ind,
const std::vector<Real> &x,
const std::vector<Real> &p)
const {
151 Real alpha(1), tmp(0);
152 ind = Bundle_U<Real>::size();
153 typename std::set<unsigned>::iterator it;
154 for (it = nworkingSet_.begin(); it != nworkingSet_.end(); it++) {
155 if ( p[*it] < -ROL_EPSILON<Real>() ) {
156 tmp = -x[*it]/p[*it];
157 if ( alpha >= tmp ) {
163 return std::max(
zero,alpha);
166template<
typename Real>
167unsigned Bundle_U_AS<Real>::solveEQPsubproblem(std::vector<Real> &s, Real &mu,
168 const std::vector<Real> &g,
const Real tol)
const {
171 unsigned n = nworkingSet_.size(), cnt = 0;
173 s.assign(Bundle_U<Real>::size(),
zero);
175 std::vector<Real> gk(n,
zero);
176 typename std::set<unsigned>::iterator it = nworkingSet_.begin();
177 for (
unsigned i = 0; i < n; ++i) {
178 gk[i] = g[*it]; it++;
180 std::vector<Real> sk(n,
zero);
181 cnt = projectedCG(sk,mu,gk,tol);
182 it = nworkingSet_.begin();
183 for (
unsigned i = 0; i < n; ++i) {
184 s[*it] = sk[i]; it++;
190template<
typename Real>
191void Bundle_U_AS<Real>::applyPreconditioner(std::vector<Real> &Px,
const std::vector<Real> &x)
const {
194 std::vector<Real> tmp(Px.size(),
zero);
196 case 0: applyPreconditioner_Identity(tmp,x);
break;
197 case 1: applyPreconditioner_Jacobi(tmp,x);
break;
198 case 2: applyPreconditioner_SymGS(tmp,x);
break;
200 applyPreconditioner_Identity(Px,tmp);
203template<
typename Real>
204void Bundle_U_AS<Real>::applyG(std::vector<Real> &Gx,
const std::vector<Real> &x)
const {
207 case 0: applyG_Identity(Gx,x);
break;
208 case 1: applyG_Jacobi(Gx,x);
break;
209 case 2: applyG_SymGS(Gx,x);
break;
213template<
typename Real>
214void Bundle_U_AS<Real>::applyPreconditioner_Identity(std::vector<Real> &Px,
const std::vector<Real> &x)
const {
215 unsigned dim = nworkingSet_.size();
216 Real sum(0), err(0), tmp(0), y(0);
217 for (
unsigned i = 0; i <
dim; ++i) {
222 err = (tmp - sum) - y;
226 for (
unsigned i = 0; i <
dim; ++i) {
231template<
typename Real>
232void Bundle_U_AS<Real>::applyG_Identity(std::vector<Real> &Gx,
const std::vector<Real> &x)
const {
233 Gx.assign(x.begin(),x.end());
236template<
typename Real>
237void Bundle_U_AS<Real>::applyPreconditioner_Jacobi(std::vector<Real> &Px,
const std::vector<Real> &x)
const {
238 const Real
zero(0), one(1);
239 unsigned dim = nworkingSet_.size();
241 Real errX(0), tmpX(0), yX(0), errE(0), tmpE(0), yE(0);
242 std::vector<Real> gg(
dim,
zero);
243 typename std::set<unsigned>::iterator it = nworkingSet_.begin();
244 for (
unsigned i = 0; i <
dim; ++i) {
245 gg[i] = one/std::abs(Bundle_U<Real>::GiGj(*it,*it)); it++;
248 yX = x[i]*gg[i] - errX;
250 errX = (tmpX - sum) - yX;
256 errE = (tmpE - eHe) - yE;
260 for (
unsigned i = 0; i <
dim; ++i) {
261 Px[i] = (x[i]-sum)*gg[i];
265template<
typename Real>
266void Bundle_U_AS<Real>::applyG_Jacobi(std::vector<Real> &Gx,
const std::vector<Real> &x)
const {
267 unsigned dim = nworkingSet_.size();
268 typename std::set<unsigned>::iterator it = nworkingSet_.begin();
269 for (
unsigned i = 0; i <
dim; ++i) {
270 Gx[i] = std::abs(Bundle_U<Real>::GiGj(*it,*it))*x[i]; it++;
274template<
typename Real>
275void Bundle_U_AS<Real>::applyPreconditioner_SymGS(std::vector<Real> &Px,
const std::vector<Real> &x)
const {
276 int dim = nworkingSet_.size();
278 gx_->zero(); ge_->zero();
279 Real eHx(0), eHe(0), one(1);
281 std::vector<Real> x1(
dim,0), e1(
dim,0),gg(
dim,0);
282 typename std::set<unsigned>::iterator it, jt;
283 it = nworkingSet_.begin();
284 for (
int i = 0; i <
dim; ++i) {
285 gx_->zero(); ge_->zero(); jt = nworkingSet_.begin();
286 for (
int j = 0; j < i; ++j) {
287 Bundle_U<Real>::addGi(*jt,x1[j],*gx_);
288 Bundle_U<Real>::addGi(*jt,e1[j],*ge_);
291 gg[i] = Bundle_U<Real>::GiGj(*it,*it);
292 x1[i] = (x[i] - Bundle_U<Real>::dotGi(*it,*gx_))/gg[i];
293 e1[i] = (one - Bundle_U<Real>::dotGi(*it,*ge_))/gg[i];
297 for (
int i = 0; i <
dim; ++i) {
302 std::vector<Real> Hx(
dim,0), He(
dim,0); it = nworkingSet_.end();
303 for (
int i =
dim-1; i >= 0; --i) {
305 gx_->zero(); ge_->zero(); jt = nworkingSet_.end();
306 for (
int j =
dim-1; j >= i+1; --j) {
308 Bundle_U<Real>::addGi(*jt,Hx[j],*gx_);
309 Bundle_U<Real>::addGi(*jt,He[j],*ge_);
311 Hx[i] = (x1[i] - Bundle_U<Real>::dotGi(*it,*gx_))/gg[i];
312 He[i] = (e1[i] - Bundle_U<Real>::dotGi(*it,*ge_))/gg[i];
318 for (
int i = 0; i <
dim; ++i) {
319 Px[i] = Hx[i] - (eHx/eHe)*He[i];
323template<
typename Real>
324void Bundle_U_AS<Real>::applyG_SymGS(std::vector<Real> &Gx,
const std::vector<Real> &x)
const {
325 unsigned dim = nworkingSet_.size();
326 typename std::set<unsigned>::iterator it = nworkingSet_.begin();
327 for (
unsigned i = 0; i <
dim; ++i) {
328 Gx[i] = std::abs(Bundle_U<Real>::GiGj(*it,*it))*x[i]; it++;
332template<
typename Real>
333void Bundle_U_AS<Real>::computeResidualUpdate(std::vector<Real> &r, std::vector<Real> &g)
const {
334 unsigned n = g.size();
335 std::vector<Real> Gg(n,0);
336 Real y(0), ytmp(0), yprt(0), yerr(0);
337 applyPreconditioner(g,r);
340 for (
unsigned i = 0; i < n; ++i) {
342 yprt = (r[i] - Gg[i]) - yerr;
344 yerr = (ytmp - y) - yprt;
348 for (
unsigned i = 0; i < n; ++i) {
351 applyPreconditioner(g,r);
354template<
typename Real>
355void Bundle_U_AS<Real>::applyFullMatrix(std::vector<Real> &Hx,
const std::vector<Real> &x)
const {
357 gx_->zero(); eG_->zero();
358 for (
unsigned i = 0; i < Bundle_U<Real>::size(); ++i) {
361 yG_->set(Bundle_U<Real>::subgradient(i)); yG_->scale(x[i]); yG_->axpy(-one,*eG_);
362 tG_->set(*gx_); tG_->plus(*yG_);
363 eG_->set(*tG_); eG_->axpy(-one,*gx_); eG_->axpy(-one,*yG_);
366 for (
unsigned i = 0; i < Bundle_U<Real>::size(); ++i) {
368 Hx[i] = Bundle_U<Real>::dotGi(i,*gx_);
372template<
typename Real>
373void Bundle_U_AS<Real>::applyMatrix(std::vector<Real> &Hx,
const std::vector<Real> &x)
const {
375 gx_->zero(); eG_->zero();
376 unsigned n = nworkingSet_.size();
377 typename std::set<unsigned>::iterator it = nworkingSet_.begin();
378 for (
unsigned i = 0; i < n; ++i) {
381 yG_->set(Bundle_U<Real>::subgradient(*it)); yG_->scale(x[i]); yG_->axpy(-one,*eG_);
382 tG_->set(*gx_); tG_->plus(*yG_);
383 eG_->set(*tG_); eG_->axpy(-one,*gx_); eG_->axpy(-one,*yG_);
387 it = nworkingSet_.begin();
388 for (
unsigned i = 0; i < n; ++i) {
390 Hx[i] = Bundle_U<Real>::dotGi(*it,*gx_); it++;
394template<
typename Real>
395unsigned Bundle_U_AS<Real>::projectedCG(std::vector<Real> &x, Real &mu,
const std::vector<Real> &b,
const Real tol)
const {
396 const Real one(1),
zero(0);
397 unsigned n = nworkingSet_.size();
398 std::vector<Real> r(n,0), r0(n,0), g(n,0), d(n,0), Ad(n,0);
402 r0.assign(r.begin(),r.end());
404 computeResidualUpdate(r,g);
405 Real rg = dot(r,g), rg0(0);
408 Real alpha(0), kappa(0), beta(0), TOL(1.e-2);
409 Real CGtol = std::min(tol,TOL*rg);
411 while (rg > CGtol && cnt < 2*n+1) {
418 computeResidualUpdate(r,g);
428 Real err(0), tmp(0), y(0);
429 for (
unsigned i = 0; i < n; ++i) {
433 err = (tmp - mu) - y;
436 mu /=
static_cast<Real
>(n);
441template<
typename Real>
442Real Bundle_U_AS<Real>::dot(
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
444 Real val(0), err(0), tmp(0), y0(0);
445 unsigned n = std::min(x.size(),y.size());
446 for (
unsigned i = 0; i < n; ++i) {
448 y0 = x[i]*y[i] - err;
450 err = (tmp - val) - y0;
456template<
typename Real>
457Real Bundle_U_AS<Real>::norm(
const std::vector<Real> &x)
const {
458 return std::sqrt(dot(x,x));
461template<
typename Real>
462void Bundle_U_AS<Real>::axpy(
const Real a,
const std::vector<Real> &x, std::vector<Real> &y)
const {
463 unsigned n = std::min(y.size(),x.size());
464 for (
unsigned i = 0; i < n; ++i) {
469template<
typename Real>
470void Bundle_U_AS<Real>::scale(std::vector<Real> &x,
const Real a)
const {
471 for (
unsigned i = 0; i < x.size(); ++i) {
476template<
typename Real>
477void Bundle_U_AS<Real>::scale(std::vector<Real> &x,
const Real a,
const std::vector<Real> &y)
const {
478 unsigned n = std::min(x.size(),y.size());
479 for (
unsigned i = 0; i < n; ++i) {
484template<
typename Real>
485unsigned Bundle_U_AS<Real>::solveDual_arbitrary(
const Real t,
const unsigned maxit,
const Real tol) {
486 const Real
zero(0), one(1);
487 initializeDualSolver();
489 unsigned ind = 0, i = 0, CGiter = 0;
490 Real snorm(0), alpha(0), mu(0);
491 std::vector<Real> s(Bundle_U<Real>::size(),0), Hs(Bundle_U<Real>::size(),0);
492 std::vector<Real> g(Bundle_U<Real>::size(),0), lam(Bundle_U<Real>::size()+1,0);
493 std::vector<Real> dualVariables(Bundle_U<Real>::size(),0);
494 for (
unsigned j = 0; j < Bundle_U<Real>::size(); ++j) {
495 dualVariables[j] = Bundle_U<Real>::getDualVariable(j);
498 Bundle_U<Real>::evaluateObjective(g,dualVariables,t);
499 for (i = 0; i < maxit; ++i) {
500 CGiter += solveEQPsubproblem(s,mu,g,tol);
502 if ( snorm < ROL_EPSILON<Real>() ) {
503 computeLagMult(lam,mu,g);
504 nonneg = isNonnegative(ind,lam);
510 if ( ind < Bundle_U<Real>::size() ) {
511 workingSet_.erase(ind);
512 nworkingSet_.insert(ind);
517 alpha = computeStepSize(ind,dualVariables,s);
518 if ( alpha >
zero ) {
519 axpy(alpha,s,dualVariables);
520 applyFullMatrix(Hs,s);
523 if (ind < Bundle_U<Real>::size()) {
524 workingSet_.insert(ind);
525 nworkingSet_.erase(ind);
533 for (
unsigned j = 0; j < Bundle_U<Real>::size(); ++j) {
534 Bundle_U<Real>::setDualVariable(j,dualVariables[j]);
539template<
typename Real>
540void Bundle_U_AS<Real>::project(std::vector<Real> &x,
const std::vector<Real> &v)
const {
541 const Real
zero(0), one(1);
542 std::vector<Real> vsort(Bundle_U<Real>::size(),0);
543 vsort.assign(v.begin(),v.end());
544 std::sort(vsort.begin(),vsort.end());
545 Real sum(-1), lam(0);
546 for (
unsigned i = Bundle_U<Real>::size()-1; i > 0; i--) {
548 if ( sum >=
static_cast<Real
>(Bundle_U<Real>::size()-i)*vsort[i-1] ) {
549 lam = sum/
static_cast<Real
>(Bundle_U<Real>::size()-i);
554 lam = (sum+vsort[0])/
static_cast<Real
>(Bundle_U<Real>::size());
556 for (
unsigned i = 0; i < Bundle_U<Real>::size(); ++i) {
557 x[i] = std::max(
zero,v[i] - lam);
561template<
typename Real>
562Real Bundle_U_AS<Real>::computeCriticality(
const std::vector<Real> &g,
const std::vector<Real> &sol)
const {
563 const Real
zero(0), one(1);
564 std::vector<Real> x(Bundle_U<Real>::size(),0), Px(Bundle_U<Real>::size(),0);
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0 zero)()
Bundle_U_AS(const unsigned maxSize=10, const Real coeff=0.0, const Real omega=2.0, const unsigned remSize=2)
Provides the interface for and implements a bundle.