diis.h
Go to the documentation of this file.
1 /*
2  * This file is a part of TiledArray.
3  * Copyright (C) 2013 Virginia Tech
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Eduard Valeyev
19  * Department of Chemistry, Virginia Tech
20  *
21  * diis.h
22  * May 20, 2013
23  *
24  */
25 
26 #ifndef TILEDARRAY_MATH_LINALG_DIIS_H__INCLUDED
27 #define TILEDARRAY_MATH_LINALG_DIIS_H__INCLUDED
28 
30 #include "TiledArray/dist_array.h"
32 
33 #include <Eigen/QR>
34 #include <deque>
35 
36 namespace TiledArray::math::linalg {
37 
39 
82 template <typename D>
83 class DIIS {
84  public:
85  typedef typename D::element_type value_type;
87  typedef Eigen::Matrix<value_type, Eigen::Dynamic, Eigen::Dynamic,
88  Eigen::RowMajor>
90  typedef Eigen::Matrix<value_type, Eigen::Dynamic, 1> Vector;
91 
93 
116  DIIS(unsigned int strt = 1, unsigned int ndi = 5, scalar_type dmp = 0,
117  unsigned int ngr = 1, unsigned int ngrdiis = 1, scalar_type mf = 0,
118  scalar_type adt = 0)
119  : error_(0),
120  errorset_(false),
121  start(strt),
122  ndiis(ndi),
123  iter(0),
124  ngroup(ngr),
125  ngroupdiis(ngrdiis),
126  damping_factor(dmp),
127  mixing_fraction(mf),
128  attenuated_damping_threshold(adt) {
129  init();
130  }
131  ~DIIS() {
132  x_.clear();
133  errors_.clear();
134  x_extrap_.clear();
135  }
136 
144  void extrapolate(D& x, D& error, bool extrapolate_error = false) {
145  iter++;
146 
147  // compute extrapolation coefficients C_ and number of skipped vectors
148  // nskip_
150 
151  // extrapolate x using above computed parameters (C_ and nskip_)
152  extrapolate(x, C_, nskip_);
153 
154  const unsigned int nvec = errors_.size();
155 
156  // sizes of the x set and the error set should equal, otherwise throw
157  TA_ASSERT(x_.size() == errors_.size() &&
158  "DIIS: numbers of guess and error vectors do not match, "
159  "likely due to a programming error");
160 
161  // extrapolate the error if needed
162  if (extrapolate_error && (mixing_fraction == 0.0 || x_extrap_.empty())) {
163  for (unsigned int k = nskip_, kk = 1; k < nvec; ++k, ++kk) {
164  axpy(error, C_[kk], errors_[k]);
165  }
166  }
167  }
168 
177  void extrapolate(D& x, const Vector& c, unsigned int nskip = 0,
178  bool increase_iter = false) {
179  if (increase_iter) {
180  iter++;
181  }
182 
183  const bool do_mixing = (mixing_fraction != 0.0);
184 
185  // if have ndiis vectors
186  if (x_.size() ==
187  ndiis) { // holding max # of vectors already? drop the least recent x
188  x_.pop_front();
189  if (not x_extrap_.empty()) x_extrap_.pop_front();
190  }
191 
192  // push x to the set
193  x_.push_back(x);
194 
195  if (iter == 1) { // the first iteration
196  if (not x_extrap_.empty() && do_mixing) {
197  zero(x);
198  axpy(x, (1.0 - mixing_fraction), x_[0]);
199  axpy(x, mixing_fraction, x_extrap_[0]);
200  }
201  } else if (iter > start && (((iter - start) % ngroup) <
202  ngroupdiis)) { // not the first iteration and
203  // need to extrapolate?
204 
205  const unsigned int nvec = x_.size();
206  const unsigned int rank = nvec - nskip + 1; // size of coefficients
207 
208  TA_ASSERT(c.size() == rank &&
209  "DIIS: numbers of coefficients and x's do not match");
210  zero(x);
211  for (unsigned int k = nskip, kk = 1; k < nvec; ++k, ++kk) {
212  if (not do_mixing || x_extrap_.empty()) {
213  // std::cout << "contrib " << k << " c=" << c[kk] << ":" << std::endl
214  // << x_[k] << std::endl;
215  axpy(x, c[kk], x_[k]);
216  } else {
217  axpy(x, c[kk] * (1.0 - mixing_fraction), x_[k]);
218  axpy(x, c[kk] * mixing_fraction, x_extrap_[k]);
219  }
220  }
221 
222  } // do DIIS
223 
224  // only need to keep extrapolated x if doing mixing
225  if (do_mixing) x_extrap_.push_back(x);
226  }
227 
233  void compute_extrapolation_parameters(const D& error,
234  bool increase_iter = false) {
235  if (increase_iter) {
236  iter++;
237  }
238 
239  // if have ndiis vectors
240  if (errors_.size() == ndiis) { // holding max # of vectors already? drop
241  // the least recent error
242  errors_.pop_front();
243  Matrix Bcrop = B_.bottomRightCorner(ndiis - 1, ndiis - 1);
244  Bcrop.conservativeResize(ndiis, ndiis);
245  B_ = Bcrop;
246  }
247 
248  // push error to the set
249  errors_.push_back(error);
250  const unsigned int nvec = errors_.size();
251 
252  // and compute the most recent elements of B, B(i,j) = <ei|ej>
253  for (unsigned int i = 0; i < nvec - 1; i++)
254  B_(i, nvec - 1) = B_(nvec - 1, i) =
255  inner_product(errors_[i], errors_[nvec - 1]);
256  B_(nvec - 1, nvec - 1) =
257  inner_product(errors_[nvec - 1], errors_[nvec - 1]);
258  using std::abs;
259  using std::sqrt;
260  const auto current_error_2norm = sqrt(abs(B_(nvec - 1, nvec - 1)));
261 
262  const scalar_type zero_determinant = 1.0e-15;
263  const scalar_type zero_norm = 1.0e-10;
264  const auto current_damping_factor =
265  attenuated_damping_threshold > 0 &&
266  current_error_2norm < attenuated_damping_threshold
267  ? damping_factor *
268  (current_error_2norm / attenuated_damping_threshold)
269  : damping_factor;
270  const scalar_type scale = 1.0 + current_damping_factor;
271 
272  // compute extrapolation coefficients C_ and number of skipped vectors
273  // nskip_
274  if (iter > start &&
275  (((iter - start) % ngroup) <
276  ngroupdiis)) { // not the first iteration and need to extrapolate?
277 
278  scalar_type absdetA;
279  nskip_ = 0; // how many oldest vectors to skip for the sake of
280  // conditioning? try zero
281  do {
282  const unsigned int rank = nvec - nskip_ + 1; // size of matrix A
283 
284  // set up the DIIS linear system: A c = rhs
285  Matrix A(rank, rank);
286  C_.resize(rank);
287 
288  A.col(0).setConstant(-1.0);
289  A.row(0).setConstant(-1.0);
290  A(0, 0) = 0.0;
291  Vector rhs = Vector::Zero(rank);
292  rhs[0] = -1.0;
293 
294  scalar_type norm = 1.0;
295  if (std::abs(B_(nskip_, nskip_)) > zero_norm)
296  norm = 1.0 / std::abs(B_(nskip_, nskip_));
297 
298  A.block(1, 1, rank - 1, rank - 1) =
299  B_.block(nskip_, nskip_, rank - 1, rank - 1) * norm;
300  A.diagonal() *= scale;
301  // for (unsigned int i=1; i < rank ; i++) {
302  // for (unsigned int j=1; j <= i ; j++) {
303  // A(i, j) = A(j, i) = B_(i+nskip-1, j+nskip-1) * norm;
304  // if (i==j) A(i, j) *= scale;
305  // }
306  //}
307 
308 #if 0
309  std::cout << "DIIS: iter=" << iter << " nskip=" << nskip << " nvec=" << nvec << std::endl;
310  std::cout << "DIIS: B=" << B_ << std::endl;
311  std::cout << "DIIS: A=" << A << std::endl;
312  std::cout << "DIIS: rhs=" << rhs << std::endl;
313 #endif
314 
315  // finally, solve the DIIS linear system
316  Eigen::ColPivHouseholderQR<Matrix> A_QR = A.colPivHouseholderQr();
317  C_ = A_QR.solve(rhs);
318  absdetA = A_QR.absDeterminant();
319 
320  // std::cout << "DIIS: |A|=" << absdetA << " sol=" << c << std::endl;
321 
322  ++nskip_;
323 
324  } while (absdetA < zero_determinant &&
325  nskip_ < nvec); // while (system is poorly conditioned)
326 
327  // failed?
328  if (absdetA < zero_determinant) {
329  std::ostringstream oss;
330  oss << "DIIS::extrapolate: poorly-conditioned system, |A| = "
331  << absdetA;
332  throw std::domain_error(oss.str());
333  }
334  --nskip_; // undo the last ++ :-(
335 
336  parameters_computed_ = true;
337  }
338  }
339 
344  if (start > iter) start = iter + 1;
345  }
346 
347  void reinitialize(const D* data = 0) {
348  iter = 0;
349  if (data) {
350  const bool do_mixing = (mixing_fraction != 0.0);
351  if (do_mixing) x_extrap_.push_front(*data);
352  }
353  }
354 
356  const Vector& get_coeffs() {
357  TA_ASSERT(parameters_computed_ && C_.size() > 0 &&
358  "DIIS: empty coefficients, because they have not been computed");
359  return C_;
360  }
361 
363  unsigned int get_nskip() { return nskip_; }
364 
367  bool parameters_computed() { return parameters_computed_; }
368 
369  private:
370  scalar_type error_;
371  bool errorset_;
372 
373  unsigned int start;
374  unsigned int ndiis;
375  unsigned int iter;
376  unsigned int ngroup;
377  unsigned int ngroupdiis;
378  scalar_type damping_factor;
379  scalar_type mixing_fraction;
380  scalar_type attenuated_damping_threshold;
381 
384  Matrix B_;
385  Vector C_;
386  bool parameters_computed_;
387  unsigned int nskip_;
389 
390  std::deque<D>
391  x_;
392  std::deque<D> errors_;
393  std::deque<D> x_extrap_;
394 
395  void set_error(scalar_type e) {
396  error_ = e;
397  errorset_ = true;
398  }
399  scalar_type error() { return error_; }
400 
401  void init() {
402  iter = 0;
403 
404  B_ = Matrix::Zero(ndiis, ndiis);
405  C_.resize(0);
406  parameters_computed_ = false;
407  nskip_ = 0;
408 
409  x_.clear();
410  errors_.clear();
411  x_extrap_.clear();
412  // x_.resize(ndiis);
413  // errors_.resize(ndiis);
414  // x_extrap_ is bigger than the other because
415  // it must hold data associated with the next iteration
416  // x_extrap_.resize(diis+1);
417  }
418 
419 }; // class DIIS
420 
421 } // namespace TiledArray::math::linalg
422 
423 namespace TiledArray {
425 }
426 
427 #endif // TILEDARRAY_MATH_LINALG_DIIS_H__INCLUDED
void reinitialize(const D *data=0)
Definition: diis.h:347
void compute_extrapolation_parameters(const D &error, bool increase_iter=false)
Definition: diis.h:233
auto inner_product(const DistArray< Tile, Policy > &a, const DistArray< Tile, Policy > &b)
Definition: dist_array.h:1647
Eigen::Matrix< value_type, Eigen::Dynamic, 1 > Vector
Definition: diis.h:90
DIIS (`‘direct inversion of iterative subspace’') extrapolation.
Definition: diis.h:83
decltype(auto) norm(const Tile< Arg > &arg)
Vector 2-norm of a tile.
Definition: tile.h:1527
void scale(DistArray< Tile, Policy > &a, S scaling_factor)
Definition: basic.h:44
auto rank(const DistArray< Tile, Policy > &a)
Definition: dist_array.h:1617
void extrapolate(D &x, const Vector &c, unsigned int nskip=0, bool increase_iter=false)
Definition: diis.h:177
Eigen::Matrix< value_type, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > Matrix
Definition: diis.h:89
#define TA_ASSERT(EXPR,...)
Definition: error.h:39
typename TiledArray::detail::scalar_type< T >::type scalar_t
scalar_t<T> is an alias for scalar_type<T>::type
Definition: type_traits.h:760
void extrapolate(D &x, D &error, bool extrapolate_error=false)
Definition: diis.h:144
TiledArray::detail::scalar_t< value_type > scalar_type
Definition: diis.h:86
auto abs(const ComplexConjugate< T > &a)
Definition: complex.h:270
void zero(DistArray< Tile, Policy > &a)
Definition: basic.h:51
void axpy(DistArray< Tile, Policy > &y, S alpha, const DistArray< Tile, Policy > &x)
Definition: basic.h:56
D::element_type value_type
Definition: diis.h:85
const Vector & get_coeffs()
calling this function returns extrapolation coefficients
Definition: diis.h:356
::Eigen::Matrix< T, ::Eigen::Dynamic, ::Eigen::Dynamic, Options > Matrix
Definition: blas.h:63
DIIS(unsigned int strt=1, unsigned int ndi=5, scalar_type dmp=0, unsigned int ngr=1, unsigned int ngrdiis=1, scalar_type mf=0, scalar_type adt=0)
Constructor.
Definition: diis.h:116
unsigned int get_nskip()
calling this function returns number of skipped vectors in extrapolation
Definition: diis.h:363