Program Listing for File UnitaryEvolver.hpp

↰ Return to documentation for file (include/Suzuki-Trotter-Evolver/UnitaryEvolver.hpp)

#include <complex>
#include <cmath>
#include <vector>
#include <Eigen/Eigenvalues>
#include <Eigen/SparseCore>
using std::complex;
using std::vector;
using Eigen::Dynamic;
using std::pow;
using std::conj;

namespace Suzuki_Trotter_Evolver {

template<int n = Dynamic, int m = Dynamic> using DMatrix =
    Eigen::Matrix<complex<double>, n, m>;

typedef Eigen::SparseMatrix<complex<double>> SMatrix;

template<int n_ctrl = Dynamic,
         int dim = Dynamic,
         typename Matrix = DMatrix<dim, dim>>
struct UnitaryEvolver {
    private: const complex<double> minus_i = complex<double>(0.0, -1.0);
    public:

    static const int dim_x_n_ctrl = (dim == Dynamic || n_ctrl == Dynamic)
                       ? Dynamic
                       : dim * n_ctrl;
    size_t length;

    Eigen::Array<complex<double>, dim, 1> d0;

    vector<Eigen::Array<complex<double>, dim, 1>> ds;

    Matrix u0;

    Matrix u0_inverse;

    vector<Matrix> us;

    vector<Matrix> us_individual;

    vector<Matrix> us_inverse_individual;

    vector<Matrix> hs;

    Matrix u0_inverse_u_last;

    template<typename T = Matrix>
    UnitaryEvolver(std::enable_if_t<std::is_same<T, DMatrix<dim, dim>>::value,
                   DMatrix<dim, dim>> drift_hamiltonian,
                   DMatrix<dim_x_n_ctrl, dim> control_hamiltonians
                  ) {
        typedef Eigen::SelfAdjointEigenSolver<DMatrix<dim, dim>> EigSolver;
        EigSolver eigs0(drift_hamiltonian);
        d0 = minus_i*eigs0.eigenvalues().array();
        u0 = eigs0.eigenvectors();
        u0_inverse = u0.adjoint();

        Eigen::Index l = drift_hamiltonian.rows();
        length = control_hamiltonians.rows()/control_hamiltonians.cols();

        if (length == 0) {
            u0_inverse_u_last = DMatrix<dim, dim>::Identity(l, l);
            us.push_back(u0);
        } else {
            hs.push_back((DMatrix<dim, dim>)(control_hamiltonians.block(0, 0, l, l)));
            EigSolver eigs(hs.back());
            ds.push_back((Eigen::Array<complex<double>, dim, 1>)(
                minus_i*eigs.eigenvalues().array()));
            DMatrix<dim, dim> u = eigs.eigenvectors();
            us_individual.push_back(u);
            us_inverse_individual.push_back(u.adjoint());
            us.push_back(u.adjoint() * u0);
            us.push_back(u); // append twice
            for (size_t i = 1; i < length; i++) {
                hs.push_back((DMatrix<dim, dim>)(control_hamiltonians.block(i*l, 0, l, l)));
                eigs = EigSolver(hs.back());
                ds.push_back((Eigen::Array<complex<double>, dim, 1>)(
                    minus_i*eigs.eigenvalues().array()));
                u = eigs.eigenvectors();
                us_individual.push_back(u);
                us_inverse_individual.push_back(u.adjoint());
                us.back() = u.adjoint() * us.back();
                us.push_back(u);
            }
            u0_inverse_u_last.noalias() = u0_inverse * u;
        }
    };

    // Sparse initialisation
    template<typename T = Matrix>
    UnitaryEvolver(std::enable_if_t<std::is_same<T, SMatrix>::value,
                   DMatrix<dim, dim>> drift_hamiltonian,
                   DMatrix<dim_x_n_ctrl, dim> control_hamiltonians) {
        typedef Eigen::SelfAdjointEigenSolver<DMatrix<dim, dim>> EigSolver;
        EigSolver eigs0(drift_hamiltonian);
        d0 = minus_i*eigs0.eigenvalues().array();
        u0 = eigs0.eigenvectors().sparseView();
        u0_inverse = u0.adjoint();

        Eigen::Index l = drift_hamiltonian.rows();
        length = control_hamiltonians.rows()/control_hamiltonians.cols();
        if (length == 0) {
            u0_inverse_u_last = DMatrix<dim, dim>::Identity(l, l).sparseView();
            us.push_back(u0);
        } else {
            hs.push_back(control_hamiltonians.block(0, 0, l, l).sparseView());
            EigSolver eigs((DMatrix<dim, dim>)(hs.back()));
            ds.push_back((Eigen::Array<complex<double>, dim, 1>)(
                                           minus_i*eigs.eigenvalues().array()));
            Matrix u = eigs.eigenvectors().sparseView();
            us_individual.push_back(u);
            us_inverse_individual.push_back(u.adjoint());
            us.push_back((u.adjoint() * u0).pruned());
            us.push_back(u); // append twice
            for (size_t i = 1; i < length; i++) {
                hs.push_back(control_hamiltonians.block(i*l, 0, l, l).sparseView());
                eigs = EigSolver((DMatrix<dim, dim>)(hs.back()));
                ds.push_back((Eigen::Array<complex<double>, dim, 1>)(
                                           minus_i*eigs.eigenvalues().array()));
                u = eigs.eigenvectors().sparseView();
                us_individual.push_back(u);
                us_inverse_individual.push_back(u.adjoint());
                us.back() = (u.adjoint() * us.back()).pruned();
                us.push_back(u);
            }
            u0_inverse_u_last = u0_inverse * u;
        }
    };

    UnitaryEvolver(size_t l,
                   Eigen::Array<complex<double>, dim, 1> d0,
                   vector<Eigen::Array<complex<double>, dim, 1>> ds,
                   Matrix u0,
                   Matrix u0_inverse,
                   vector<Matrix> us,
                   vector<Matrix> us_individual,
                   vector<Matrix> us_inverse_individual,
                   vector<Matrix> hs,
                   Matrix u0_inverse_u_last)
        : length{l},
          d0{d0},
          ds{ds},
          u0{u0},
          u0_inverse{u0_inverse},
          us{us},
          us_individual{us_individual},
          us_inverse_individual{us_inverse_individual},
          hs{hs},
          u0_inverse_u_last{u0_inverse_u_last}
        {};

    DMatrix<dim, 1> propagate(DMatrix<Dynamic, n_ctrl> ctrl_amp,
                              DMatrix<dim, 1> state,
                              double dt) {
        size_t steps = ctrl_amp.rows();
        DMatrix<dim, 1> exp_d0 = (d0*dt).exp().matrix();

        ctrl_amp *= dt;
        state = u0_inverse * state;
        state = state.cwiseProduct(exp_d0);
        for (size_t j = 0; j < length; j++) {
            state = us[j] * state;
            state = state.cwiseProduct((ds[j]*ctrl_amp(0, j)).exp().matrix());
        }
        for (size_t i = 1; i < steps; i++) {
            state = u0_inverse_u_last * state;
            state = state.cwiseProduct(exp_d0);
            for (size_t j = 0; j < length; j++) {
                state = us[j] * state;
                state = state.cwiseProduct(
                    (ds[j]*ctrl_amp(i, j)).exp().matrix());
            }
        }
        state = us[length] * state;
        return state;
    };

    template<int l = Dynamic>
    DMatrix<dim, l> propagate_collection(DMatrix<Dynamic, n_ctrl> ctrl_amp,
                                         DMatrix<dim, l> states,
                                         double dt) {
        size_t steps = ctrl_amp.rows();
        Eigen::Array<complex<double>, dim, 1> exp_d0 = (d0*dt).exp();

        ctrl_amp *= dt;
        states = u0_inverse*states;
        states = states.array().colwise() * exp_d0;
        for (size_t j = 0; j < length; j++) {
            states = us[j] * states;
            states = states.array().colwise() * (ds[j]*ctrl_amp(0, j)).exp();
        }
        for (size_t i = 1; i < steps; i++) {
            states = u0_inverse_u_last * states;
            states = states.array().colwise() * exp_d0;
            for (size_t j = 0; j < length; j++) {
                states = us[j] * states;
                states = states.array().colwise()
                        *(ds[j]*ctrl_amp(i, j)).exp();
            }
        }
        states = us[length] * states;
        return states;
    };

    DMatrix<dim, Dynamic> propagate_all(DMatrix<Dynamic, n_ctrl> ctrl_amp,
                                      DMatrix<dim, 1> state,
                                      double dt) {
        size_t steps = ctrl_amp.rows();
        DMatrix<dim, 1> exp_d0 = (d0*dt).exp().matrix();
        DMatrix<dim, Dynamic> phi(state.rows(), steps+1);
        ctrl_amp *= dt;
        phi.col(0) = state;
        phi.col(1) = u0_inverse * state;
        phi.col(1) = phi.col(1).cwiseProduct(exp_d0);
        for (size_t j = 0; j < length; j++) {
            phi.col(1) = us[j] * phi.col(1);
            phi.col(1) = phi.col(1).cwiseProduct(
                (ds[j]*ctrl_amp(0, j)).exp().matrix());
        }
        phi.col(1) = us[length] * phi.col(1);
        for (size_t i = 1, k = 2; i < steps; i++, k++) {
            phi.col(k) = u0_inverse * phi.col(i);
            phi.col(k) = phi.col(k).cwiseProduct(exp_d0);
            for (size_t j = 0; j < length; j++) {
                phi.col(k) = us[j] * phi.col(k);
                phi.col(k) = phi.col(k).cwiseProduct(
                    (ds[j]*ctrl_amp(i, j)).exp().matrix());
            }
            phi.col(k) = us[length] * phi.col(k);
        }
        return phi;
    };

    complex<double> evolved_expectation_value(DMatrix<Dynamic, n_ctrl> ctrl_amp,
                                              DMatrix<dim, 1> state,
                                              double dt,
                                              DMatrix<dim, dim> observable) {
        state = propagate(ctrl_amp, state, dt);
        return state.dot(observable * state); // dot conjugates and transposes
    };

    DMatrix<Dynamic, 1> evolved_expectation_value_all(
                                              DMatrix<Dynamic, n_ctrl> ctrl_amp,
                                              DMatrix<dim, 1> state,
                                              double dt,
                                              DMatrix<dim, dim> observable) {
        DMatrix<dim, Dynamic> phi = propagate_all(ctrl_amp, state, dt);
        return (phi.adjoint() * observable * phi).diagonal();
    };

    std::tuple<complex<double>, Eigen::Matrix<double, Dynamic, n_ctrl>>
    switching_function(DMatrix<Dynamic, n_ctrl> ctrl_amp,
                       DMatrix<dim, 1> state,
                       double dt,
                       DMatrix<dim, dim> cost) {
        // forward propagation
        size_t steps = ctrl_amp.rows();
        DMatrix<dim, 1> exp_d0 = (d0*dt).exp().matrix();
        ctrl_amp *= dt;
        state = u0_inverse * state;
        state = state.cwiseProduct(exp_d0);
        for (size_t j = 0; j < length; j++) {
            state = us[j] * state;
            state = state.cwiseProduct((ds[j]*ctrl_amp(0, j)).exp().matrix());
        }
        for (size_t i = 1; i < steps; i++) {
            state = u0_inverse_u_last * state;
            state = state.cwiseProduct(exp_d0);
            for (size_t j = 0; j < length; j++) {
                state = us[j] * state;
                state = state.cwiseProduct(
                                         (ds[j]*ctrl_amp(i, j)).exp().matrix());
            }
        }
        state = us[length] * state;
        // Lambda
        DMatrix<dim, 1> lambda = cost * state;
        // Energy
        complex<double> E = state.dot(lambda);
        // Back propagation
        DMatrix<dim, 1> exp_d0_c = exp_d0.conjugate();
        ctrl_amp *= -1;
        DMatrix<Dynamic, n_ctrl> out = DMatrix<Dynamic, n_ctrl>::Zero(steps,
                                                                      length);
        for (size_t k = steps-1; true;) {
            for (size_t j = length-1; true;) {
                out(k, j) = lambda.dot(hs[j] * state);
                state = us_inverse_individual[j] * state;
                state = state.cwiseProduct(
                                         (ds[j]*ctrl_amp(k, j)).exp().matrix());
                state = us_individual[j] * state;
                lambda = us_inverse_individual[j] * lambda;
                lambda = lambda.cwiseProduct(
                    (ds[j]*ctrl_amp(k, j)).exp().matrix());
                lambda = us_individual[j] * lambda;
                if (j == 0) {
                    break;
                }
                j--;
            }
            if (k == 0) {
                break;
            }
            state = u0_inverse * state;
            state = state.cwiseProduct(exp_d0_c);
            state = u0 * state;
            lambda = u0_inverse * lambda;
            lambda = lambda.cwiseProduct(exp_d0_c);
            lambda = u0 * lambda;
            k--;
        }
        return
            std::tuple<complex<double>, Eigen::Matrix<double, Dynamic, n_ctrl>>(
                E,
                2*out.imag()
        );
    };
};
}