diis.h

Go to the documentation of this file.

/*
 *  This file is a part of TiledArray.
 *  Copyright (C) 2013  Virginia Tech
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 *  Eduard Valeyev
 *  Department of Chemistry, Virginia Tech
 *
 *  diis.h
 *  May 20, 2013
 *
 */

#ifndef TILEDARRAY_ALGEBRA_DIIS_H__INCLUDED
#define TILEDARRAY_ALGEBRA_DIIS_H__INCLUDED

#include <deque>
#include <TiledArray/math/eigen.h>
#include <TiledArray/algebra/utils.h>
#include "../dist_array.h"

namespace TiledArray {


  template <typename D>
  class DIIS {
    public:
      typedef typename D::element_type value_type;
      typedef typename detail::scalar_t<value_type> scalar_type;
      typedef Eigen::Matrix<value_type, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> EigenMatrixX;
      typedef Eigen::Matrix<value_type, Eigen::Dynamic, 1> EigenVectorX;


      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) :
             error_(0), errorset_(false),
             start(strt), ndiis(ndi),
             iter(0), ngroup(ngr),
             ngroupdiis(ngrdiis),
             damping_factor(dmp),
             mixing_fraction(mf)
           {
            init();
           }
      ~DIIS() {
        x_.clear();
        errors_.clear();
        x_extrap_.clear();
      }

      void extrapolate(D& x,
                       D& error,
                       bool extrapolate_error = false)
      {
        iter++;

        // compute extrapolation coefficients C_ and number of skipped vectors nskip_
        compute_extrapolation_parameters(error);

        // extrapolate x using above computed parameters (C_ and nskip_)
        extrapolate(x, C_, nskip_);

        const unsigned int nvec = errors_.size();

        // sizes of the x set and the error set should equal, otherwise throw
        TA_USER_ASSERT(x_.size() == errors_.size(),
                       "DIIS: numbers of guess and error vectors do not match, likely due to a programming error");

        // extrapolate the error if needed
        if (extrapolate_error && (mixing_fraction == 0.0 || x_extrap_.empty())) {
          for (unsigned int k=nskip_, kk=1; k < nvec; ++k, ++kk) {
            axpy(error, C_[kk], errors_[k]);
          }
        }
      }

      void extrapolate(D& x,
                       const EigenVectorX &c,
                       unsigned int nskip = 0,
                       bool increase_iter = false) {
        if (increase_iter) {
          iter++;
        }

        const bool do_mixing = (mixing_fraction != 0.0);

        // if have ndiis vectors
        if (x_.size() == ndiis) { // holding max # of vectors already? drop the least recent x
          x_.pop_front();
          if (not x_extrap_.empty()) x_extrap_.pop_front();
        }

        // push x to the set
        x_.push_back(x);

        if (iter == 1) { // the first iteration
          if (not x_extrap_.empty() && do_mixing) {
            zero(x);
            axpy(x, (1.0-mixing_fraction), x_[0]);
            axpy(x, mixing_fraction, x_extrap_[0]);
          }
        }
        else if (iter > start && (((iter - start) % ngroup) < ngroupdiis)) { // not the first iteration and need to extrapolate?

          const unsigned int nvec = x_.size();
          const unsigned int rank = nvec - nskip + 1; // size of coefficients

          TA_USER_ASSERT(c.size() == rank,
                         "DIIS: numbers of coefficients and x's do not match");
          zero(x);
          for (unsigned int k=nskip, kk=1; k < nvec; ++k, ++kk) {
            if (not do_mixing || x_extrap_.empty()) {
              //std::cout << "contrib " << k << " c=" << c[kk] << ":" << std::endl << x_[k] << std::endl;
              axpy(x, c[kk], x_[k]);
            } else {
              axpy(x, c[kk] * (1.0 - mixing_fraction), x_[k]);
              axpy(x, c[kk] * mixing_fraction, x_extrap_[k]);
            }
          }

        } // do DIIS

        // only need to keep extrapolated x if doing mixing
        if (do_mixing) x_extrap_.push_back(x);
      }

      void compute_extrapolation_parameters(const D &error,
                                            bool increase_iter = false) {
        if (increase_iter) {
          iter++;
        }

        const scalar_type zero_determinant = 1.0e-15;
        const scalar_type zero_norm = 1.0e-10;
        const scalar_type scale = 1.0 + damping_factor;

        // if have ndiis vectors
        if (errors_.size() == ndiis) { // holding max # of vectors already? drop the least recent error
          errors_.pop_front();
          EigenMatrixX Bcrop = B_.bottomRightCorner(ndiis-1,ndiis-1);
          Bcrop.conservativeResize(ndiis,ndiis);
          B_ = Bcrop;
        }

        // push error to the set
        errors_.push_back(error);
        const unsigned int nvec = errors_.size();

        // and compute the most recent elements of B, B(i,j) = <ei|ej>
        for (unsigned int i=0; i < nvec-1; i++)
          B_(i,nvec-1) = B_(nvec-1,i) = dot_product(errors_[i], errors_[nvec-1]);
        B_(nvec-1,nvec-1) = dot_product(errors_[nvec-1], errors_[nvec-1]);

        // compute extrapolation coefficients C_ and number of skipped vectors nskip_
        if (iter > start && (((iter - start) % ngroup) < ngroupdiis)) { // not the first iteration and need to extrapolate?

          scalar_type absdetA;
          nskip_ = 0; // how many oldest vectors to skip for the sake of conditioning?
                                        // try zero
          do {
            const unsigned int rank = nvec - nskip_ + 1; // size of matrix A

            // set up the DIIS linear system: A c = rhs
            EigenMatrixX A(rank, rank);
            C_.resize(rank);

            A.col(0).setConstant(-1.0);
            A.row(0).setConstant(-1.0);
            A(0,0) = 0.0;
            EigenVectorX rhs = EigenVectorX::Zero(rank);
            rhs[0] = -1.0;

            scalar_type norm = 1.0;
            if (std::abs(B_(nskip_,nskip_)) > zero_norm)
              norm = 1.0/std::abs(B_(nskip_,nskip_));

            A.block(1, 1, rank-1, rank-1) = B_.block(nskip_, nskip_, rank-1, rank-1) * norm;
            A.diagonal() *= scale;
            //for (unsigned int i=1; i < rank ; i++) {
            //  for (unsigned int j=1; j <= i ; j++) {
            //    A(i, j) = A(j, i) = B_(i+nskip-1, j+nskip-1) * norm;
            //    if (i==j) A(i, j) *= scale;
            //  }
            //}

#if 0
            std::cout << "DIIS: iter=" << iter << " nskip=" << nskip << " nvec=" << nvec << std::endl;
            std::cout << "DIIS: B=" << B_ << std::endl;
            std::cout << "DIIS: A=" << A << std::endl;
            std::cout << "DIIS: rhs=" << rhs << std::endl;
#endif

            // finally, solve the DIIS linear system
            Eigen::ColPivHouseholderQR<EigenMatrixX> A_QR = A.colPivHouseholderQr();
            C_ = A_QR.solve(rhs);
            absdetA = A_QR.absDeterminant();

            //std::cout << "DIIS: |A|=" << absdetA << " sol=" << c << std::endl;

            ++nskip_;

          } while (absdetA < zero_determinant && nskip_ < nvec); // while (system is poorly conditioned)

          // failed?
          if (absdetA < zero_determinant) {
            std::ostringstream oss;
            oss << "DIIS::extrapolate: poorly-conditioned system, |A| = " << absdetA;
            throw std::domain_error(oss.str());
          }
          --nskip_; // undo the last ++ :-(

          parameters_computed_ = true;
        }

      }

      void start_extrapolation() {
        if (start > iter) start = iter+1;
      }

      void reinitialize(const D* data = 0) {
        iter=0;
        if (data) {
          const bool do_mixing = (mixing_fraction != 0.0);
          if (do_mixing) x_extrap_.push_front(*data);
        }
      }

      const EigenVectorX &get_coeffs() {
        TA_USER_ASSERT(parameters_computed_ && C_.size() > 0,
                       "DIIS: empty coefficients, because they have not been computed");
        return C_;
      }

      unsigned int get_nskip() { return nskip_; }

      bool parameters_computed() { return parameters_computed_; }

    private:
      scalar_type error_;
      bool errorset_;

      unsigned int start;
      unsigned int ndiis;
      unsigned int iter;
      unsigned int ngroup;
      unsigned int ngroupdiis;
      scalar_type damping_factor;
      scalar_type mixing_fraction;

      EigenMatrixX B_; 
      EigenVectorX C_; 
      bool parameters_computed_; 
      unsigned int nskip_; 

      std::deque<D> x_; 
      std::deque<D> errors_; 
      std::deque<D> x_extrap_; 

      void set_error(scalar_type e) { error_ = e; errorset_ = true; }
      scalar_type error() { return error_; }

      void init() {
        iter = 0;

        B_ = EigenMatrixX::Zero(ndiis,ndiis);
        C_.resize(0);
        parameters_computed_ = false;
        nskip_ = 0;

        x_.clear();
        errors_.clear();
        x_extrap_.clear();
        //x_.resize(ndiis);
        //errors_.resize(ndiis);
        // x_extrap_ is bigger than the other because
        // it must hold data associated with the next iteration
        //x_extrap_.resize(diis+1);
      }

  }; // class DIIS

} // namespace TiledArray

#endif // TILEDARRAY_ALGEBRA_DIIS_H__INCLUDED