Line data    Source code

       1             : // SPDX-FileCopyrightText: 2024 Pairinteraction Developers
       2             : // SPDX-License-Identifier: LGPL-3.0-or-later
       3             : 
       4             : #include "pairinteraction/interfaces/DiagonalizerInterface.hpp"
       5             : 
       6             : #include "pairinteraction/enums/FloatType.hpp"
       7             : #include "pairinteraction/utils/eigen_assertion.hpp"
       8             : #include "pairinteraction/utils/eigen_compat.hpp"
       9             : 
      10             : #include <Eigen/Dense>
      11             : #include <Eigen/SparseCore>
      12             : #include <optional>
      13             : #include <spdlog/spdlog.h>
      14             : 
      15             : namespace pairinteraction {
      16             : template <typename Scalar>
      17          45 : DiagonalizerInterface<Scalar>::DiagonalizerInterface(FloatType float_type)
      18          45 :     : float_type(float_type) {}
      19             : 
      20             : template <typename Scalar>
      21             : template <typename ScalarLim>
      22             : Eigen::MatrixX<ScalarLim>
      23         340 : DiagonalizerInterface<Scalar>::subtract_mean(const Eigen::MatrixX<Scalar> &matrix, real_t &shift,
      24             :                                              double rtol) const {
      25             :     using real_lim_t = typename traits::NumTraits<ScalarLim>::real_t;
      26         340 :     int dim = matrix.rows();
      27             : 
      28         340 :     assert(rtol > 0);
      29         340 :     assert(rtol < 1);
      30             : 
      31             :     // Subtract the mean of the diagonal elements from the diagonal
      32         340 :     shift = matrix.diagonal().real().mean();
      33         340 :     Eigen::MatrixX<ScalarLim> shifted_matrix =
      34         332 :         (matrix - shift * Eigen::MatrixX<Scalar>::Identity(dim, dim)).template cast<ScalarLim>();
      35             : 
      36             :     // Ensure the accuracy the results
      37         340 :     double floating_point_error = 5 * std::numeric_limits<real_lim_t>::epsilon();
      38             : 
      39         340 :     if (floating_point_error > rtol) {
      40           0 :         SPDLOG_WARN(
      41             :             "Because the floating point precision is too low, the eigenvalues cannot be calculated "
      42             :             "accurately. The estimated floating point error ({} * ||H||) is larger than the "
      43             :             "specified tolerance ({} * ||H||). Try to use a 'float_type' with higher precision or "
      44             :             "a larger 'rtol'.",
      45             :             floating_point_error, rtol);
      46             :     }
      47             : 
      48             :     // Note that this check is not as handwavy as it seems at first glance:
      49             :     //
      50             :     // Let us first analyze the error of the eigenvalues due to the limited floating point
      51             :     // precision. Under the assumption that the diagonalization routine is backward stable and the
      52             :     // we know that the calculated eigenvalues and eigenvectors correspond to the diagonalization of
      53             :     // a matrix H+E, where H is the original matrix and E a small perturbation with
      54             :     // ||E|| <= c*eps*||H||. Here, for an hermitian matrix, c is a modest constant and eps is the
      55             :     // machine epsilon. Taking this result, we can use the Bauer-Fike theorem, simplified for an
      56             :     // hermitian matrix H, to bound the error of the eigenvalues:
      57             :     //
      58             :     //     |lambda - lambda_exact| <= c*eps*||H||,
      59             :     //
      60             :     // where lambda is the calculated eigenvalue and lambda_exact is the closest exact eigenvalue of
      61             :     // H. Thus, to ensure that the relative tolerance rtol can be met, we need c*eps*||H|| <=
      62             :     // rtol*||H||, which we check in the code above.
      63             :     //
      64             :     // If we know that for a user, it is permissible that |lambda - lambda_exact| ~ rtol*||H||, we
      65             :     // assume that the eigenvectors are allowed to have an error as well. To obtain a reasonable
      66             :     // bound for the strength of this error, we demand that the error of an eigenvalue, which would
      67             :     // result from calculating the eigenvalue by transforming the Hamiltonian with the
      68             :     // eigenvector v, should be smaller than rtol * ||H||. For v = v_exact + w, we obtain:
      69             :     //
      70             :     //     |v^dag * H * v - lambda_exact| <= 2*||w||*||H||
      71             :     //
      72             :     // Thus, to keep this error below rtol * ||H||, we need ||w|| <= rtol/2.
      73             :     //
      74             :     // In the diagonalizer classes derived from the DiagonalizerInterface, we will make use of the
      75             :     // fact that a small error of the eigenvectors is tolerable to increase the sparsity of the
      76             :     // eigenvector matrix by setting every entry smaller than rtol/(2*sqrt(dim(H))) to zero.
      77             :     // The factor 1/sqrt(dim(H)) ensures that the bound is met also if multiple entries are
      78             :     // set to zero.
      79             : 
      80         680 :     return shifted_matrix;
      81           0 : }
      82             : 
      83             : template <typename Scalar>
      84             : template <typename RealLim>
      85             : Eigen::VectorX<typename DiagonalizerInterface<Scalar>::real_t>
      86         340 : DiagonalizerInterface<Scalar>::add_mean(const Eigen::VectorX<RealLim> &shifted_eigenvalues,
      87             :                                         real_t shift) const {
      88             :     // Add the mean of the diagonal elements back to the eigenvalues
      89         340 :     Eigen::VectorX<real_t> eigenvalues = shifted_eigenvalues.template cast<real_t>();
      90         340 :     eigenvalues.array() += shift;
      91         340 :     return eigenvalues;
      92           0 : }
      93             : 
      94             : template <typename Scalar>
      95             : EigenSystemH<Scalar>
      96          24 : DiagonalizerInterface<Scalar>::eigh(const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> &matrix,
      97             :                                     std::optional<real_t> min_eigenvalue,
      98             :                                     std::optional<real_t> max_eigenvalue, double rtol) const {
      99          24 :     int dim = matrix.rows();
     100             : 
     101          24 :     auto eigensys = eigh(matrix, rtol);
     102             : 
     103             :     auto *it_begin =
     104          24 :         std::lower_bound(eigensys.eigenvalues.data(), eigensys.eigenvalues.data() + dim,
     105          24 :                          min_eigenvalue.value_or(std::numeric_limits<real_t>::lowest() / 2));
     106             :     auto *it_end =
     107          24 :         std::upper_bound(eigensys.eigenvalues.data(), eigensys.eigenvalues.data() + dim,
     108          24 :                          max_eigenvalue.value_or(std::numeric_limits<real_t>::max() / 2));
     109          24 :     eigensys.eigenvectors = eigensys.eigenvectors
     110          24 :                                 .block(0, std::distance(eigensys.eigenvalues.data(), it_begin), dim,
     111             :                                        std::distance(it_begin, it_end))
     112             :                                 .eval();
     113          24 :     eigensys.eigenvalues = eigensys.eigenvalues
     114          24 :                                .segment(std::distance(eigensys.eigenvalues.data(), it_begin),
     115             :                                         std::distance(it_begin, it_end))
     116             :                                .eval();
     117             : 
     118          24 :     return eigensys;
     119           0 : }
     120             : 
     121             : // Explicit instantiations
     122             : template class DiagonalizerInterface<double>;
     123             : template class DiagonalizerInterface<std::complex<double>>;
     124             : 
     125             : template Eigen::MatrixX<float>
     126             : DiagonalizerInterface<double>::subtract_mean(const Eigen::MatrixX<double> &matrix, double &shift,
     127             :                                              double rtol) const;
     128             : template Eigen::MatrixX<std::complex<float>>
     129             : DiagonalizerInterface<std::complex<double>>::subtract_mean(
     130             :     const Eigen::MatrixX<std::complex<double>> &matrix, double &shift, double rtol) const;
     131             : 
     132             : template Eigen::VectorX<double>
     133             : DiagonalizerInterface<double>::add_mean(const Eigen::VectorX<float> &shifted_eigenvalues,
     134             :                                         double shift) const;
     135             : template Eigen::VectorX<double> DiagonalizerInterface<std::complex<double>>::add_mean(
     136             :     const Eigen::VectorX<float> &shifted_eigenvalues, double shift) const;
     137             : 
     138             : template Eigen::MatrixX<double>
     139             : DiagonalizerInterface<double>::subtract_mean(const Eigen::MatrixX<double> &matrix, double &shift,
     140             :                                              double rtol) const;
     141             : template Eigen::MatrixX<std::complex<double>>
     142             : DiagonalizerInterface<std::complex<double>>::subtract_mean(
     143             :     const Eigen::MatrixX<std::complex<double>> &matrix, double &shift, double rtol) const;
     144             : 
     145             : template Eigen::VectorX<double>
     146             : DiagonalizerInterface<double>::add_mean(const Eigen::VectorX<double> &shifted_eigenvalues,
     147             :                                         double shift) const;
     148             : template Eigen::VectorX<double> DiagonalizerInterface<std::complex<double>>::add_mean(
     149             :     const Eigen::VectorX<double> &shifted_eigenvalues, double shift) const;
     150             : } // namespace pairinteraction