Unverified Commit cf7ca388 authored by Mccaskey, Alex's avatar Mccaskey, Alex Committed by GitHub

Merge pull request #236 from tnguyen-ornl/tnguyen/qLanczos

Implemented QLanczos algorithm plugin
parents 84eaaff9 942e7cd9
Pipeline #112777 passed with stage
in 15 minutes and 14 seconds
# Simple 1-qubit demonstration of the Quatum Imaginary Time Evolution algorithm
# Simple 1-qubit demonstration of the Quatum Imaginary Time Evolution/QLanczos algorithm
# Reference: https://www.nature.com/articles/s41567-019-0704-4
# Target H = 1/sqrt(2)(X + Z)
# Expected minimum value: -1.0
......@@ -12,9 +12,6 @@ qpu = xacc.getAccelerator('qpp')
# Construct the Hamiltonian as an XACC PauliOperator
ham = xacc.getObservable('pauli', '0.70710678118 X0 + 0.70710678118 Z0')
# We just need 1 qubit
buffer = xacc.qalloc(1)
# See Fig. 2 (e) of https://www.nature.com/articles/s41567-019-0704-4
# Horizontal axis: 0 -> 2.5
# The number of Trotter steps
......@@ -30,13 +27,27 @@ qite = xacc.getAlgorithm('qite', {
'step-size': stepSize,
'steps': nbSteps
})
# We just need 1 qubit
qiteBuffer = xacc.qalloc(1)
qiteResult = qite.execute(qiteBuffer)
# Create the QLanczos algorithm
qLanczos = xacc.getAlgorithm('QLanczos', {
'accelerator': qpu,
'observable': ham,
'step-size': stepSize,
'steps': nbSteps
})
qLanczosBuffer = xacc.qalloc(1)
qLanczosResult = qLanczos.execute(qLanczosBuffer)
result = qite.execute(buffer)
# Plot the energies
qiteEnergies = qiteBuffer.getInformation('exp-vals')
qLanczosEnergies = qLanczosBuffer.getInformation('exp-vals')
# Expected result: ~ -1
print('Min energy value = ', buffer.getInformation('opt-val'))
E = buffer.getInformation('exp-vals')
# Reproduce Fig. 2(e) of https://www.nature.com/articles/s41567-019-0704-4
plt.plot(np.arange(0, nbSteps + 1) * stepSize, E, 'ro', label = 'XACC QITE')
plt.plot(np.arange(0, nbSteps + 1) * stepSize, qiteEnergies, 'ro-', label = 'QITE')
plt.plot(np.arange(0, (nbSteps + 1)//2) * 2 * stepSize, qLanczosEnergies, 'k*-', label = 'QLanczos')
plt.legend()
plt.grid()
plt.show()
\ No newline at end of file
......@@ -33,8 +33,10 @@ public:
/**
*/
void Start(BundleContext context) {
auto c = std::make_shared<xacc::algorithm::QITE>();
context.RegisterService<xacc::Algorithm>(c);
auto qite = std::make_shared<xacc::algorithm::QITE>();
context.RegisterService<xacc::Algorithm>(qite);
auto qlanczos = std::make_shared<xacc::algorithm::QLanczos>();
context.RegisterService<xacc::Algorithm>(qlanczos);
}
/**
......
......@@ -134,6 +134,17 @@ arma::cx_mat createSMatrix(const std::vector<std::string>& in_pauliOps, const st
}
return S_Mat;
}
template<typename T>
std::vector<T> arange(T start, T stop, T step = 1)
{
std::vector<T> values;
for (T value = start; value < stop; value += step)
{
values.emplace_back(value);
}
return values;
}
}
using namespace xacc;
......@@ -303,11 +314,11 @@ double QITE::calcCurrentEnergy(int in_nbQubits) const
auto expval = buffers[i]->getExpectationValueZ();
energy += expval * coefficients[i];
}
std::cout << "Energy = " << energy << "\n";
xacc::info("Energy = " + std::to_string(energy));
return energy;
}
std::shared_ptr<Observable> QITE::calcAOps(const std::shared_ptr<AcceleratorBuffer>& in_buffer, std::shared_ptr<CompositeInstruction> in_kernel, std::shared_ptr<Observable> in_hmTerm) const
std::pair<double, std::shared_ptr<Observable>> QITE::calcAOps(const std::shared_ptr<AcceleratorBuffer>& in_buffer, std::shared_ptr<CompositeInstruction> in_kernel, std::shared_ptr<Observable> in_hmTerm) const
{
const auto pauliOps = generatePauliPermutation(in_buffer->size());
......@@ -409,7 +420,7 @@ std::shared_ptr<Observable> QITE::calcAOps(const std::shared_ptr<AcceleratorBuff
// which emulate the imaginary time evolution of the original observable.
std::shared_ptr<Observable> updatedAham = std::make_shared<xacc::quantum::PauliOperator>();
updatedAham->fromString(aObsStr);
return updatedAham;
return std::make_pair(std::sqrt(c), updatedAham);
}
void QITE::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const
......@@ -431,7 +442,7 @@ void QITE::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const
// Propagates the state via Trotter steps:
auto kernel = constructPropagateCircuit();
// Optimizes/calculates next A ops
auto nextAOps = calcAOps(buffer, kernel, hamOp);
auto [normVal, nextAOps] = calcAOps(buffer, kernel, hamOp);
m_approxOps.emplace_back(nextAOps);
m_energyAtStep.emplace_back(calcCurrentEnergy(buffer->size()));
}
......@@ -607,5 +618,199 @@ std::vector<double> QITE::execute(const std::shared_ptr<AcceleratorBuffer> buffe
xacc::error("This method is unsupported!");
return {};
}
bool QLanczos::initialize(const HeterogeneousMap& parameters)
{
// Calls the base (QITE) initialization
if (QITE::initialize(parameters))
{
// Gets some QLanczos-specific parameters if provided:
// Default params:
m_sLim = 0.75;
m_epsLim = 1e-12;
if (parameters.keyExists<double>("s-param"))
{
m_sLim = parameters.get<double>("s-param");
}
if (parameters.keyExists<double>("epsilon-param"))
{
m_epsLim = parameters.get<double>("epsilon-param");
}
return true;
}
return false;
}
void QLanczos::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const
{
std::vector<double> normAtStep;
std::vector<double> lanczosEnergy;
// Run on hardware/simulator using quantum gates/measure
// Initial energy
m_energyAtStep.emplace_back(calcCurrentEnergy(buffer->size()));
// Initial norm
normAtStep.emplace_back(1.0);
auto hamOp = std::make_shared<xacc::quantum::PauliOperator>();
for (const auto& hamTerm : m_observable->getNonIdentitySubTerms())
{
*hamOp = *hamOp + *(std::dynamic_pointer_cast<xacc::quantum::PauliOperator>(hamTerm));
}
// Time stepping:
for (int i = 0; i < m_nbSteps; ++i)
{
// Propagates the state via Trotter steps:
auto kernel = constructPropagateCircuit();
// Optimizes/calculates next A ops
auto [normVal, nextAOps] = calcAOps(buffer, kernel, hamOp);
m_approxOps.emplace_back(nextAOps);
m_energyAtStep.emplace_back(calcCurrentEnergy(buffer->size()));
normAtStep.emplace_back(normVal * normAtStep.back());
// Even steps:
if ((i % 2) == 0)
{
lanczosEnergy.emplace_back(calcQlanczosEnergy(normAtStep));
}
}
assert(m_energyAtStep.size() == m_nbSteps + 1);
// Last QLanczos energy value
buffer->addExtraInfo("opt-val", ExtraInfo(lanczosEnergy.back()));
// Also returns the full list of energy values
// at each QLanczos step.
buffer->addExtraInfo("exp-vals", ExtraInfo(lanczosEnergy));
}
double QLanczos::calcQlanczosEnergy(const std::vector<double>& normVec) const
{
const std::vector<size_t> lanczosSteps = arange(1UL, m_energyAtStep.size() + 1, 2UL);
const auto n = lanczosSteps.size();
// H and S matrices (Eq. 60)
arma::mat H(n, n, arma::fill::zeros);
arma::mat S(n, n, arma::fill::zeros);
int j = 0;
int k = 0;
// Iterate over l and l'
for (const auto& l : lanczosSteps)
{
int k = 0;
for (const auto& lp : lanczosSteps)
{
// Defining 2r = l + l'
const auto r = (l + lp) / 2;
// Eq. 61
S(j, k) = normVec[r] * normVec[r]/(normVec[l] * normVec[lp]);
// Eq. 62
H(j, k) = m_energyAtStep[r] * S(j, k);
k++;
}
j++;
}
const auto matFieldSampling = [](const arma::mat& in_H, const arma::mat& in_S, const std::vector<size_t> in_indexVec) {
arma::mat H(in_indexVec.size(), in_indexVec.size(), arma::fill::zeros);
arma::mat S(in_indexVec.size(), in_indexVec.size(), arma::fill::zeros);
assert(in_S.n_cols == in_H.n_cols && in_S.n_rows == in_H.n_rows);
for (size_t i = 0; i < in_H.n_rows; ++i)
{
for (size_t j = 0; j < in_H.n_cols; ++j)
{
if (xacc::container::contains(in_indexVec, i) && xacc::container::contains(in_indexVec, j))
{
const size_t rowIdx = std::distance(in_indexVec.begin(), std::find(in_indexVec.begin(), in_indexVec.end(), i));
const size_t colIdx = std::distance(in_indexVec.begin(), std::find(in_indexVec.begin(), in_indexVec.end(), j));
S(rowIdx, colIdx) = in_S(i, j);
H(rowIdx, colIdx) = in_H(i, j);
}
}
}
return std::make_pair(H, S);
};
// Regularize/stabilize H and S matrices with s and epsilon stabilization parameters.
// Returns the stabilized (H, S) matrices.
const auto regularizeMatrices = [&](double s, double eps) -> std::pair<arma::mat, arma::mat> {
std::vector<size_t> indexVec { 0 };
size_t ii = 0;
size_t jj = 0;
while (ii < H.n_rows && jj < (H.n_rows - 1))
{
for(jj = ii + 1; jj < H.n_rows; ++jj)
{
if(S(ii, jj) < s)
{
indexVec.emplace_back(jj);
break;
}
}
ii = indexVec.back();
}
if (!xacc::container::contains(indexVec, H.n_rows - 1))
{
indexVec.emplace_back(H.n_rows - 1);
}
auto [Hnew, Snew] = matFieldSampling(H, S, indexVec);
// Handles an edge case where Hnew and Snew are just single-element matrices;
// just returns those matrices.
if (indexVec.size() == 1)
{
assert(Hnew.n_elem == 1 && Snew.n_elem == 1);
// Just returns these matrices,
// no need to regularize any further.
return std::make_pair(Hnew, Snew);
}
// Truncates eigenvalues if less than epsilon
arma::vec sigma;
arma::mat V;
arma::eig_sym(sigma, V, Snew);
std::vector<size_t> indexVecEig;
for (size_t i = 0; i < sigma.n_elem; ++i)
{
if (sigma[i] > eps)
{
indexVecEig.emplace_back(i);
}
}
const auto eigenTransform = [](const arma::mat& in_mat, const arma::mat& in_eigenMat){
// Performs V_T * Matrix * V;
// where V is the eigenvector matrix.
arma::mat result = in_eigenMat.t() * in_mat;
result = result * in_eigenMat;
return result;
};
const arma::mat Snew2 = eigenTransform(Snew, V);
const arma::mat Hnew2 = eigenTransform(Hnew, V);
auto [Hnew3, Snew3] = matFieldSampling(Hnew2, Snew2, indexVecEig);
return std::make_pair(Hnew3, Snew3);
};
auto [Hreg, Sreg] = regularizeMatrices(m_sLim, m_epsLim);
arma::cx_vec eigval;
arma::cx_mat eigvec;
// Solves the generalized eigen val
arma::eig_pair(eigval, eigvec, Hreg, Sreg);
std::vector<double> energies;
for (size_t i = 0; i < eigval.n_elem; ++i)
{
// Energy values should be real
assert(std::abs(eigval(i).imag()) < 1e-9);
energies.emplace_back(eigval(i).real());
}
std::sort(energies.begin(), energies.end());
// Returns the smallest eigenvalue from the previous step.
return energies.front();
}
} // namespace algorithm
} // namespace xacc
\ No newline at end of file
......@@ -27,7 +27,7 @@ public:
const std::string name() const override { return "qite"; }
const std::string description() const override { return ""; }
DEFINE_ALGORITHM_CLONE(QITE)
private:
protected:
// Construct the Trotter propagate circuit up to current time step.
std::shared_ptr<CompositeInstruction> constructPropagateCircuit() const;
// Calculate the current energy, i.e.
......@@ -39,9 +39,10 @@ private:
// in_kernel: the kernel to evolve the system to this time step
// in_hmTerm: the H term to be approximate by the A term
// i.e. emulate the imaginary time evolution of that H term.
std::shared_ptr<Observable> calcAOps(const std::shared_ptr<AcceleratorBuffer>& in_buffer, std::shared_ptr<CompositeInstruction> in_kernel, std::shared_ptr<Observable> in_hmTerm) const;
// Returns the norm (as a double) and the A operator (Pauli observable)
std::pair<double, std::shared_ptr<Observable>> calcAOps(const std::shared_ptr<AcceleratorBuffer>& in_buffer, std::shared_ptr<CompositeInstruction> in_kernel, std::shared_ptr<Observable> in_hmTerm) const;
private:
protected:
// Number of Trotter steps
int m_nbSteps;
// dBeta, i.e. step size
......@@ -66,5 +67,22 @@ private:
// prepare the initial state.
int m_initialState;
};
// QLanczos Algorithm: extends QITE and
// typically provides better energy convergence.
class QLanczos : public QITE {
public:
bool initialize(const HeterogeneousMap &parameters) override;
const std::string name() const override { return "QLanczos"; }
const std::string description() const override { return ""; }
void execute(const std::shared_ptr<AcceleratorBuffer> buffer) const override;
DEFINE_ALGORITHM_CLONE(QLanczos)
private:
double calcQlanczosEnergy(const std::vector<double>& normVec) const;
private:
// Regularize parameters:
double m_sLim;
double m_epsLim;
};
} // namespace algorithm
} // namespace xacc
......@@ -12,4 +12,7 @@
# *******************************************************************************/
include_directories(${CMAKE_BINARY_DIR})
add_xacc_test(QITE)
target_link_libraries(QITETester xacc xacc-pauli)
\ No newline at end of file
target_link_libraries(QITETester xacc xacc-pauli)
add_xacc_test(QLanczos)
target_link_libraries(QLanczosTester xacc xacc-pauli)
\ No newline at end of file
/*******************************************************************************
* Copyright (c) 2020 UT-Battelle, LLC.
* All rights reserved. This program and the accompanying materials
* are made available under the terms of the Eclipse Public License v1.0
* and Eclipse Distribution License v1.0 which accompanies this
* distribution. The Eclipse Public License is available at
* http://www.eclipse.org/legal/epl-v10.html and the Eclipse Distribution
*License is available at https://eclipse.org/org/documents/edl-v10.php
*
* Contributors:
* Thien Nguyen - initial API and implementation
*******************************************************************************/
#include <gtest/gtest.h>
#include "xacc.hpp"
#include "xacc_service.hpp"
#include "PauliOperator.hpp"
TEST(QLanczosTester, checkSimple)
{
auto qLanczos = xacc::getService<xacc::Algorithm>("QLanczos");
std::shared_ptr<xacc::Observable> observable = std::make_shared<xacc::quantum::PauliOperator>();
observable->fromString("0.7071067811865475 X0 + 0.7071067811865475 Z0");
auto acc = xacc::getAccelerator("qpp");
const int nbSteps = 25;
const double stepSize = 0.1;
const bool initOk = qLanczos->initialize({
std::make_pair("accelerator", acc),
std::make_pair("steps", nbSteps),
std::make_pair("observable", observable),
std::make_pair("step-size", stepSize)
});
EXPECT_TRUE(initOk);
auto buffer = xacc::qalloc(1);
qLanczos->execute(buffer);
const double finalEnergy = (*buffer)["opt-val"].as<double>();
std::cout << "Final Energy: " << finalEnergy << "\n";
// Fig (2.e) of https://www.nature.com/articles/s41567-019-0704-4
// Minimal Energy = -1
EXPECT_NEAR(finalEnergy, -1.0, 1e-3);
const std::vector<double> energyValues = (*buffer)["exp-vals"].as<std::vector<double>>();
// Even steps only
EXPECT_LT(energyValues.size(), (nbSteps + 1) / 2 + 1);
// Should converge must faster (after 2 Lanczos steps)
// see Fig. 2(e)
for (int i = 1; i < energyValues.size(); ++i)
{
EXPECT_LT(energyValues[i], 0.9);
}
}
int main(int argc, char **argv) {
xacc::Initialize(argc, argv);
::testing::InitGoogleTest(&argc, argv);
auto ret = RUN_ALL_TESTS();
xacc::Finalize();
return ret;
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment