davidsonsolver.cc

Go to the documentation of this file.

/*
 * Copyright 2009-2020 The VOTCA Development Team (http://www.votca.org)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
 
// Standard includes
#include <iostream>
#include <stdexcept>
 
// Local VOTCA includes
#include "votca/xtp/davidsonsolver.h"
#include "votca/xtp/eigen.h"
 
namespace votca {
namespace xtp {
 
DavidsonSolver::DavidsonSolver(Logger &log) : log_(log) {}
 
void DavidsonSolver::printTiming(
    const std::chrono::time_point<std::chrono::system_clock> &start) const {
  XTP_LOG(Log::error, log_)
      << TimeStamp() << "-----------------------------------" << std::flush;
  std::chrono::time_point<std::chrono::system_clock> end =
      std::chrono::system_clock::now();
  std::chrono::duration<double> elapsed_time = end - start;
  XTP_LOG(Log::error, log_) << TimeStamp() << "- Davidson ran for "
                            << elapsed_time.count() << "secs." << std::flush;
  XTP_LOG(Log::error, log_)
      << TimeStamp() << "-----------------------------------" << std::flush;
}
 
void DavidsonSolver::checkOptions(Index operator_size) {
  //. search space exceeding the system size
  if (max_search_space_ > operator_size) {
    XTP_LOG(Log::error, log_)
        << TimeStamp() << " == Warning : Max search space ("
        << max_search_space_ << ") larger than system size (" << operator_size
        << ")" << std::flush;
 
    max_search_space_ = operator_size;
    XTP_LOG(Log::error, log_)
        << TimeStamp() << " == Warning : Max search space set to "
        << operator_size << std::flush;
 
    XTP_LOG(Log::error, log_)
        << TimeStamp()
        << " == Warning : If problems appear, try asking for less than "
        << Index(operator_size / 10) << " eigenvalues" << std::flush;
  }
}
 
void DavidsonSolver::printOptions(Index operator_size) const {
 
  XTP_LOG(Log::error, log_)
      << TimeStamp() << " Davidson Solver using " << OPENMP::getMaxThreads()
      << " threads." << std::flush;
  XTP_LOG(Log::error, log_)
      << TimeStamp() << " Tolerance : " << tol_ << std::flush;
 
  switch (this->davidson_correction_) {
    case CORR::DPR:
      XTP_LOG(Log::error, log_)
          << TimeStamp() << " DPR Correction" << std::flush;
      break;
    case CORR::OLSEN:
      XTP_LOG(Log::error, log_)
          << TimeStamp() << " Olsen Correction" << std::flush;
      break;
  }
 
  XTP_LOG(Log::error, log_) << TimeStamp() << " Matrix size : " << operator_size
                            << 'x' << operator_size << std::flush;
}
 
void DavidsonSolver::printIterationData(
    const DavidsonSolver::RitzEigenPair &rep,
    const DavidsonSolver::ProjectedSpace &proj, Index neigen) const {
 
  Index converged_roots = proj.root_converged.head(neigen).count();
  double percent_converged = 100 * double(converged_roots) / double(neigen);
  XTP_LOG(Log::error, log_)
      << TimeStamp()
      << boost::format(" %1$4d %2$12d \t %3$4.2e \t %4$5.2f%% converged") %
             i_iter_ % proj.search_space() %
             rep.res_norm().head(neigen).maxCoeff() % percent_converged
      << std::flush;
}
 
void DavidsonSolver::set_matrix_type(std::string mt) {
  if (mt == "HAM") {
    this->matrix_type_ = MATRIX_TYPE::HAM;
  } else if (mt == "SYMM") {
    this->matrix_type_ = MATRIX_TYPE::SYMM;
  } else {
    throw std::runtime_error(mt + " is not a valid Davidson matrix type");
  }
}
 
void DavidsonSolver::set_correction(std::string method) {
  if (method == "DPR") {
    this->davidson_correction_ = CORR::DPR;
  } else if (method == "OLSEN") {
    this->davidson_correction_ = CORR::OLSEN;
  } else {
    throw std::runtime_error(method +
                             " is not a valid Davidson correction method");
  }
}
 
void DavidsonSolver::set_tolerance(std::string tol) {
  if (tol == "loose") {
    this->tol_ = 1E-3;
  } else if (tol == "normal") {
    this->tol_ = 1E-4;
  } else if (tol == "strict") {
    this->tol_ = 1E-5;
  } else if (tol == "lapack") {
    this->tol_ = 1E-9;
  } else {
    throw std::runtime_error(tol + " is not a valid Davidson tolerance");
  }
}
 
void DavidsonSolver::set_size_update(std::string update_size) {
 
  if (update_size == "min") {
    this->davidson_update_ = UPDATE::MIN;
  } else if (update_size == "safe") {
    this->davidson_update_ = UPDATE::SAFE;
  } else if (update_size == "max") {
    this->davidson_update_ = UPDATE::MAX;
  } else {
    throw std::runtime_error(update_size + " is not a valid Davidson update");
  }
}
 
Index DavidsonSolver::getSizeUpdate(Index neigen) const {
  Index size_update;
  switch (this->davidson_update_) {
    case UPDATE::MIN:
      size_update = neigen;
      break;
    case UPDATE::SAFE:
      if (neigen < 20) {
        size_update = static_cast<Index>(1.5 * double(neigen));
      } else {
        size_update = neigen + 10;
      }
      break;
    case UPDATE::MAX:
      size_update = 2 * neigen;
      break;
    default:
      size_update = 2 * neigen;
      break;
  }
  return size_update;
}
 
ArrayXl DavidsonSolver::argsort(const Eigen::VectorXd &V) const {
  /* \brief return the index of the sorted vector */
  ArrayXl idx = ArrayXl::LinSpaced(V.rows(), 0, V.rows() - 1);
  std::sort(idx.data(), idx.data() + idx.size(),
            [&](Index i1, Index i2) { return V[i1] < V[i2]; });
  return idx;
}
 
Eigen::MatrixXd DavidsonSolver::setupInitialEigenvectors(
    Index size_initial_guess) const {
 
  Eigen::MatrixXd guess =
      Eigen::MatrixXd::Zero(Adiag_.size(), size_initial_guess);
  ArrayXl idx = DavidsonSolver::argsort(Adiag_);
 
  switch (this->matrix_type_) {
    case MATRIX_TYPE::SYMM:
      /* \brief Initialize the guess eigenvector so that they 'target' the
       * smallest diagonal elements */
      for (Index j = 0; j < size_initial_guess; j++) {
        guess(idx(j), j) = 1.0;
      }
      break;
 
    case MATRIX_TYPE::HAM:
      /* Initialize the guess eigenvector so that they 'target' the lowest
       * positive diagonal elements */
      Index ind0 = Adiag_.size() / 2;
      for (Index j = 0; j < size_initial_guess; j++) {
        guess(idx(ind0 + j), j) = 1.0;
      }
      break;
  }
  return guess;
}
DavidsonSolver::RitzEigenPair DavidsonSolver::getRitzEigenPairs(
    const ProjectedSpace &proj) const {
  // get the ritz vectors
  switch (this->matrix_type_) {
    case MATRIX_TYPE::SYMM: {
      return getRitz(proj);
    }
    case MATRIX_TYPE::HAM: {
      return getHarmonicRitz(proj);
    }
  }
  return RitzEigenPair();
}
 
DavidsonSolver::RitzEigenPair DavidsonSolver::getRitz(
    const DavidsonSolver::ProjectedSpace &proj) const {
 
  DavidsonSolver::RitzEigenPair rep;
  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(proj.T);
  if (es.info() != Eigen::ComputationInfo::Success) {
    std::cerr << "A\n" << proj.T << std::endl;
    throw std::runtime_error("Small hermitian eigenvalue problem failed.");
  }
  // we only need enough pairs for either extension of space or restart
  Index needed_pairs =
      std::min(proj.T.cols(), std::max(restart_size_, proj.size_update));
  rep.lambda = es.eigenvalues().head(needed_pairs);
  rep.U = es.eigenvectors().leftCols(needed_pairs);
 
  rep.q = proj.V * rep.U;                                       // Ritz vectors
  rep.res = proj.AV * rep.U - rep.q * rep.lambda.asDiagonal();  // residues
  return rep;
}
 
DavidsonSolver::RitzEigenPair DavidsonSolver::getHarmonicRitz(
    const ProjectedSpace &proj) const {
 
  /* Compute the Harmonic Ritz vector following
   * Computing Interior Eigenvalues of Large Matrices
   * Ronald B Morgan
   * LINEAR ALGEBRA AND ITS APPLICATIONS 154-156:289-309 (1991)
   * https://cpb-us-w2.wpmucdn.com/sites.baylor.edu/dist/e/71/files/2015/05/InterEvals-1vgdz91.pdf
   */
 
  RitzEigenPair rep;
  bool return_eigenvectors = true;
  Eigen::GeneralizedEigenSolver<Eigen::MatrixXd> ges(proj.T, proj.B,
                                                     return_eigenvectors);
  if (ges.info() != Eigen::ComputationInfo::Success) {
    std::cerr << "A\n" << proj.T << std::endl;
    std::cerr << "B\n" << proj.B << std::endl;
    throw std::runtime_error("Small generalized eigenvalue problem failed.");
  }
 
  std::vector<std::pair<Index, Index>> complex_pairs;
  for (Index i = 0; i < ges.eigenvalues().size(); i++) {
    if (ges.eigenvalues()(i).imag() != 0) {
      bool found_partner = false;
      for (auto &pair : complex_pairs) {
        if (pair.second > -1) {
          continue;
        } else {
          bool are_pair = (std::abs(ges.eigenvalues()(pair.first).real() -
                                    ges.eigenvalues()(i).real()) < 1e-9) &&
                          (std::abs(ges.eigenvalues()(pair.first).imag() +
                                    ges.eigenvalues()(i).imag()) < 1e-9);
          if (are_pair) {
            pair.second = i;
            found_partner = true;
          }
        }
      }
 
      if (!found_partner) {
        complex_pairs.emplace_back(i, -1);
      }
    }
  }
 
  for (const auto &pair : complex_pairs) {
    if (pair.second < 0) {
      throw std::runtime_error("Eigenvalue:" + std::to_string(pair.first) +
                               " is complex but has no partner.");
    }
  }
  if (!complex_pairs.empty()) {
    XTP_LOG(Log::warning, log_)
        << TimeStamp() << " Found " << complex_pairs.size()
        << " complex pairs in eigenvalue problem" << std::flush;
  }
  Eigen::VectorXd eigenvalues =
      Eigen::VectorXd::Zero(ges.eigenvalues().size() - complex_pairs.size());
  Eigen::MatrixXd eigenvectors =
      Eigen::MatrixXd::Zero(ges.eigenvectors().rows(),
                            ges.eigenvectors().cols() - complex_pairs.size());
 
  Index j = 0;
  for (Index i = 0; i < ges.eigenvalues().size(); i++) {
    bool is_second_in_complex_pair =
        std::find_if(complex_pairs.begin(), complex_pairs.end(),
                     [&](const std::pair<Index, Index> &pair) {
                       return pair.second == i;
                     }) != complex_pairs.end();
    if (is_second_in_complex_pair) {
      continue;
    } else {
      eigenvalues(j) = ges.eigenvalues()(i).real();
      eigenvectors.col(j) = ges.eigenvectors().col(i).real();
      eigenvectors.col(j).normalize();
      j++;
    }
  }
  // we only need enough pairs for either extension of space or restart
  Index needed_pairs =
      std::min(proj.T.cols(), std::max(restart_size_, proj.size_update));
  // Dropping one member of each complex pair can leave fewer candidates than
  // needed (e.g. right after a restart, when the space has restart_size_
  // columns). Then add the imaginary parts of the pairs' eigenvectors: with
  // the real parts they span the real invariant subspaces of the pairs.
  if (j < needed_pairs) {
    const Index missing =
        std::min(needed_pairs - j, Index(complex_pairs.size()));
    eigenvalues.conservativeResize(j + missing);
    eigenvectors.conservativeResize(Eigen::NoChange, j + missing);
    for (Index p = 0; p < missing; ++p) {
      const Index i = complex_pairs[p].first;
      eigenvalues(j) = ges.eigenvalues()(i).real();
      eigenvectors.col(j) = ges.eigenvectors().col(i).imag();
      eigenvectors.col(j).normalize();
      j++;
    }
  }
  ArrayXl idx =
      DavidsonSolver::argsort(eigenvalues).reverse().head(needed_pairs);
  // we need the largest values, because this is the inverse value, so
  // reverse list
 
  rep.U = DavidsonSolver::extract_vectors(eigenvectors, idx);
  rep.lambda = (rep.U.transpose() * proj.T * rep.U).diagonal();
  rep.q = proj.V * rep.U;                                       // Ritz vectors
  rep.res = proj.AV * rep.U - rep.q * rep.lambda.asDiagonal();  // residues
  return rep;
}
 
DavidsonSolver::ProjectedSpace DavidsonSolver::initProjectedSpace(
    Index neigen, Index size_initial_guess) const {
  DavidsonSolver::ProjectedSpace proj;
 
  // initial vector basis
  proj.V = DavidsonSolver::setupInitialEigenvectors(size_initial_guess);
 
  // update variables
  proj.size_update = DavidsonSolver::getSizeUpdate(neigen);
  proj.root_converged = ArrayXb::Constant(proj.size_update, false);
  return proj;
}
 
bool DavidsonSolver::checkConvergence(const DavidsonSolver::RitzEigenPair &rep,
                                      DavidsonSolver::ProjectedSpace &proj,
                                      Index neigen) const {
  // there can be fewer Ritz pairs than size_update while the search space is
  // still small (an initial guess with dependent columns dropped); roots
  // without a Ritz pair are not converged
  const Index n = std::min<Index>(proj.size_update, rep.res_norm().size());
  proj.root_converged = ArrayXb::Constant(proj.size_update, false);
  proj.root_converged.head(n) = (rep.res_norm().head(n) < tol_);
  return proj.root_converged.head(neigen).all();
}
 
Index DavidsonSolver::extendProjection(
    const DavidsonSolver::RitzEigenPair &rep,
    DavidsonSolver::ProjectedSpace &proj) const {
 
  // correction vectors only for roots with a Ritz pair (see
  // checkConvergence)
  const Index n = std::min<Index>(proj.size_update, rep.lambda.size());
  Index nupdate = (proj.root_converged.head(n) == false).count();
  Index oldsize = proj.V.cols();
  proj.V.conservativeResize(Eigen::NoChange, oldsize + nupdate);
 
  Index k = 0;
  for (Index j = 0; j < n; j++) {
    // skip the roots that have already converged
    if (proj.root_converged[j]) {
      continue;
    }
    // residue vector
    Eigen::VectorXd w =
        computeCorrectionVector(rep.q.col(j), rep.lambda(j), rep.res.col(j));
 
    // append the correction vector to the search space
    proj.V.col(oldsize + k) = w.normalized();
    k++;
  }
  const Index added = orthogonalize(proj.V, nupdate);
  if (added == 0) {
    throw std::runtime_error(
        "Davidson: no new search direction, all corrections are linearly "
        "dependent on the search space");
  }
  return added;
}
 
Eigen::MatrixXd DavidsonSolver::extract_vectors(const Eigen::MatrixXd &V,
                                                const ArrayXl &idx) const {
  Eigen::MatrixXd W = Eigen::MatrixXd::Zero(V.rows(), idx.size());
  for (Index i = 0; i < idx.size(); i++) {
    W.col(i) = V.col(idx(i));
  }
  return W;
}
 
Eigen::VectorXd DavidsonSolver::computeCorrectionVector(
    const Eigen::VectorXd &qj, double lambdaj,
    const Eigen::VectorXd &Aqj) const {
 
  /* compute correction vector with either DPR or OLSEN CORRECTION
   * For details on the method see :
   * Systematic Study of Selected Diagonalization Methods
   * for Configuration Interaction Matrices
   * M.L. Leininger et al .
   * Journal of Computational Chemistry Vol 22, No. 13 1574-1589 (2001)
   */
  Eigen::VectorXd correction;
  switch (this->davidson_correction_) {
    case CORR::DPR: {
      correction = dpr(Aqj, lambdaj);
      break;
    }
    case CORR::OLSEN: {
      correction = olsen(Aqj, qj, lambdaj);
      break;
    }
  }
  // make sure no nan values are there, instead we set them to zero
  return correction.unaryExpr(
      [](double v) { return std::isfinite(v) ? v : 0.0; });
}
 
Eigen::VectorXd DavidsonSolver::dpr(const Eigen::VectorXd &r,
                                    double lambda) const {
  /* \brief Compute the diagonal preconditoned residue :
  \delta = -r/(D - lambda)
  with |D - lambda| floored, so that a diagonal element at the Ritz value
  gives a large but finite component instead of inf/nan
   */
  constexpr double floor = 1e-8;
  return -r.array() / (Adiag_.array() - lambda).unaryExpr([](double d) {
    return (std::abs(d) < floor) ? (d < 0 ? -floor : floor) : d;
  });
}
 
Eigen::VectorXd DavidsonSolver::olsen(const Eigen::VectorXd &r,
                                      const Eigen::VectorXd &x,
                                      double lambda) const {
  /* \brief Compute the olsen correction :
  \delta = (D-\lambda)^{-1} (-r + \epsilon x)
  */
  Eigen::VectorXd delta = DavidsonSolver::dpr(r, lambda);
  double num = -x.transpose() * delta;
  double denom = -x.transpose() * dpr(x, lambda);
  if (!(std::abs(denom) > 1e-14 * std::abs(num)) || !std::isfinite(denom)) {
    return delta;  // no usable Olsen term: DPR
  }
  double eps = num / denom;
  delta += eps * x;
  return delta;
}
 
Index DavidsonSolver::orthogonalize(Eigen::MatrixXd &V, Index nupdate) const {
  return DavidsonSolver::gramschmidt(V, V.cols() - nupdate);
}
 
// Orthonormalises the columns from nstart on against the first nstart
// (assumed orthonormal) and against each other, classical Gram-Schmidt twice.
// Columns that are linearly dependent (less than 1e-10 of their norm left)
// are dropped; returns the number of columns kept from nstart on.
Index DavidsonSolver::gramschmidt(Eigen::MatrixXd &Q, Index nstart) const {
  const Index nupdate = Q.cols() - nstart;
  if (nupdate <= 0) {
    return 0;
  }
  const Eigen::VectorXd norms = Q.rightCols(nupdate).colwise().norm();
  auto project_out_old = [&](Index first, Index n) {
    if (nstart > 0 && n > 0) {
      Q.middleCols(first, n) -=
          Q.leftCols(nstart) *
          (Q.leftCols(nstart).transpose() * Q.middleCols(first, n));
    }
  };
  // against the existing vectors, as blocks
  project_out_old(nstart, nupdate);
  project_out_old(nstart, nupdate);
  // against each other, dropping dependent vectors
  Index kept = 0;
  for (Index j = 0; j < nupdate; ++j) {
    Eigen::VectorXd v = Q.col(nstart + j);
    for (Index pass = 0; pass < 2 && kept > 0; ++pass) {
      v -= Q.middleCols(nstart, kept) *
           (Q.middleCols(nstart, kept).transpose() * v);
    }
    const double norm = v.norm();
    if (!(norm > 1e-10 * norms(j))) {
      continue;
    }
    Q.col(nstart + kept) = v / norm;
    ++kept;
  }
  Q.conservativeResize(Eigen::NoChange, nstart + kept);
  // the normalisation can amplify what is left of the existing vectors
  project_out_old(nstart, kept);
  Q.rightCols(kept).colwise().normalize();
  if (kept < nupdate) {
    XTP_LOG(Log::info, log_)
        << TimeStamp() << " Davidson: dropped " << nupdate - kept
        << " linearly dependent search vector(s)" << std::flush;
  }
  return kept;
}
 
Eigen::MatrixXd DavidsonSolver::qr(const Eigen::MatrixXd &A) const {
 
  Index nrows = A.rows();
  Index ncols = A.cols();
  ncols = std::min(nrows, ncols);
  Eigen::MatrixXd I = Eigen::MatrixXd::Identity(nrows, ncols);
  Eigen::HouseholderQR<Eigen::MatrixXd> qr(A);
  return qr.householderQ() * I;
}
 
void DavidsonSolver::restart(const DavidsonSolver::RitzEigenPair &rep,
                             DavidsonSolver::ProjectedSpace &proj,
                             Index newvectors) const {
  Eigen::MatrixXd newV =
      Eigen::MatrixXd(proj.V.rows(), newvectors + restart_size_);
  newV.rightCols(newvectors) = proj.V.rightCols(newvectors);
  if (matrix_type_ == MATRIX_TYPE::SYMM) {
 
    newV.leftCols(restart_size_) = rep.q.leftCols(restart_size_);
    proj.AV *= rep.U.leftCols(restart_size_);  // corresponds to replacing
                                               // V with q.leftCols
  } else {
    Eigen::MatrixXd orthonormal =
        DavidsonSolver::qr(rep.U.leftCols(restart_size_));
    newV.leftCols(restart_size_) =
        proj.V.leftCols(proj.V.cols() - newvectors) * orthonormal;
    proj.AV *= orthonormal;
 
    proj.AAV *= orthonormal;
    proj.B = newV.leftCols(restart_size_).transpose() * proj.AAV;
  }
  proj.T = newV.leftCols(restart_size_).transpose() * proj.AV;
  proj.V = newV;
}
 
void DavidsonSolver::storeConvergedData(
    const DavidsonSolver::RitzEigenPair &rep, Index neigen) {
 
  DavidsonSolver::storeEigenPairs(rep, neigen);
  XTP_LOG(Log::error, log_) << TimeStamp() << " Davidson converged after "
                            << i_iter_ << " iterations." << std::flush;
  info_ = Eigen::ComputationInfo::Success;
}
 
void DavidsonSolver::storeNotConvergedData(
    const DavidsonSolver::RitzEigenPair &rep, const ArrayXb &root_converged,
    Index neigen) {
 
  DavidsonSolver::storeEigenPairs(rep, neigen);
 
  double percent_converged = 0;
 
  for (Index i = 0; i < neigen; i++) {
    if (!root_converged[i]) {
      eigenvalues_(i) = 0;
      eigenvectors_.col(i).setZero();
    } else {
      percent_converged += 1.;
    }
  }
  percent_converged /= double(neigen);
  percent_converged *= 100.;
  XTP_LOG(Log::error, log_)
      << TimeStamp() << "- Warning : Davidson "
      << boost::format("%1$5.2f%%") % percent_converged << " converged after "
      << i_iter_ << " iterations." << std::flush;
  info_ = Eigen::ComputationInfo::NoConvergence;
}
 
void DavidsonSolver::storeEigenPairs(const DavidsonSolver::RitzEigenPair &rep,
                                     Index neigen) {
  // store the eigenvalues/eigenvectors
  this->eigenvalues_ = rep.lambda.head(neigen);
  this->eigenvectors_ = rep.q.leftCols(neigen);
  this->eigenvectors_.colwise().normalize();
}
 
}  // namespace xtp
}  // namespace votca