diff --git a/libs/solvers/include/cudaq/solvers/adapt.h b/libs/solvers/include/cudaq/solvers/adapt.h index 4fc58665..cef0598f 100644 --- a/libs/solvers/include/cudaq/solvers/adapt.h +++ b/libs/solvers/include/cudaq/solvers/adapt.h @@ -70,7 +70,23 @@ class adapt_impl : public extension_point { /// @param pool Pool of operators /// @param optimizer Optimization algorithm /// @param gradient Gradient calculation method - /// @param options Additional options for the algorithm + /// @param options Additional options for the algorithm. Supported Keys: + /// - "max_iter" (int): Maximum number of iterations [default: 30] + /// - "grad_norm_tolerance" (double): Convergence tolerance for gradient norm + /// [default: 1e-5] + /// - "grad_norm_diff_tolerance" (double): Tolerance for difference between + /// gradient norms [default: 1e-5] + /// - "threshold_energy" (double): Energy convergence threshold [default: + /// 1e-6] + /// - "initial_theta" (double): Initial value for theta parameter [default: + /// 0.0] + /// - "verbose" (bool): Enable detailed output logging [default: false] + /// - "shots" (int): Number of measurement shots (-1 for exact simulation) + /// [default: -1] + /// - "dynamic_start" (string): Optimization mode for the theta parameters at + /// each iteration. + /// It can be either "warm", or "cold". [default: "cold"] + /// - "tol" (double): Tolerance for optimization [default: 1e-12] /// @return Result of the ADAPT-VQE algorithm virtual result run(const cudaq::qkernel &)> &initState, const spin_op &H, const std::vector &pool, @@ -88,7 +104,22 @@ class adapt_impl : public extension_point { /// @param initialState Initial state preparation quantum kernel /// @param H Hamiltonian operator /// @param poolList Pool of operators -/// @param options Additional options for the algorithm +/// @param options Additional options for the algorithm. Supported Keys: +/// - "max_iter" (int): Maximum number of iterations [default: 30] +/// - "grad_norm_tolerance" (double): Convergence tolerance for gradient norm +/// [default: 1e-5] +/// - "grad_norm_diff_tolerance" (double): Tolerance for difference between +/// gradient norms [default: 1e-5] +/// - "threshold_energy" (double): Energy convergence threshold [default: 1e-6] +/// - "initial_theta" (double): Initial value for theta parameter [default: +/// 0.0] +/// - "verbose" (bool): Enable detailed output logging [default: false] +/// - "shots" (int): Number of measurement shots (-1 for exact simulation) +/// [default: -1] +/// - "dynamic_start" (string): Optimization mode for the theta parameters at +/// each iteration. +/// It can be either "warm", or "cold". [default: "cold"] +/// - "tol" (double): Tolerance for optimization [default: 1e-12] /// @return Result of the ADAPT-VQE algorithm static inline adapt::result adapt_vqe(const cudaq::qkernel &)> &initialState, @@ -106,7 +137,22 @@ adapt_vqe(const cudaq::qkernel &)> &initialState, /// @param H Hamiltonian operator /// @param poolList Pool of operators /// @param optimizer Custom optimization algorithm -/// @param options Additional options for the algorithm +/// @param options Additional options for the algorithm. Supported Keys: +/// - "max_iter" (int): Maximum number of iterations [default: 30] +/// - "grad_norm_tolerance" (double): Convergence tolerance for gradient norm +/// [default: 1e-5] +/// - "grad_norm_diff_tolerance" (double): Tolerance for difference between +/// gradient norms [default: 1e-5] +/// - "threshold_energy" (double): Energy convergence threshold [default: 1e-6] +/// - "initial_theta" (double): Initial value for theta parameter [default: +/// 0.0] +/// - "verbose" (bool): Enable detailed output logging [default: false] +/// - "shots" (int): Number of measurement shots (-1 for exact simulation) +/// [default: -1] +/// - "dynamic_start" (string): Optimization mode for the theta parameters at +/// each iteration. +/// It can be either "warm", or "cold". [default: "cold"] +/// - "tol" (double): Tolerance for optimization [default: 1e-12] /// @return Result of the ADAPT-VQE algorithm static inline adapt::result adapt_vqe(const cudaq::qkernel &)> &initialState, @@ -125,7 +171,22 @@ adapt_vqe(const cudaq::qkernel &)> &initialState, /// @param poolList Pool of operators /// @param optimizer Custom optimization algorithm /// @param gradient Gradient calculation method -/// @param options Additional options for the algorithm +/// @param options Additional options for the algorithm. Supported Keys: +/// - "max_iter" (int): Maximum number of iterations [default: 30] +/// - "grad_norm_tolerance" (double): Convergence tolerance for gradient norm +/// [default: 1e-5] +/// - "grad_norm_diff_tolerance" (double): Tolerance for difference between +/// gradient norms [default: 1e-5] +/// - "threshold_energy" (double): Energy convergence threshold [default: 1e-6] +/// - "initial_theta" (double): Initial value for theta parameter [default: +/// 0.0] +/// - "verbose" (bool): Enable detailed output logging [default: false] +/// - "shots" (int): Number of measurement shots (-1 for exact simulation) +/// [default: -1] +/// - "dynamic_start" (string): Optimization mode for the theta parameters at +/// each iteration. +/// It can be either "warm", or "cold". [default: "cold"] +/// - "tol" (double): Tolerance for optimization [default: 1e-12] /// @return Result of the ADAPT-VQE algorithm static inline adapt::result adapt_vqe(const cudaq::qkernel &)> &initialState, diff --git a/libs/solvers/include/cudaq/solvers/operators/operator_pools/qaoa_operator_pool.h b/libs/solvers/include/cudaq/solvers/operators/operator_pools/qaoa_operator_pool.h index 67224380..f839d843 100644 --- a/libs/solvers/include/cudaq/solvers/operators/operator_pools/qaoa_operator_pool.h +++ b/libs/solvers/include/cudaq/solvers/operators/operator_pools/qaoa_operator_pool.h @@ -12,16 +12,33 @@ namespace cudaq::solvers { -/// @brief Test +/// @brief QAOA operator pool is a class with generates a pool +/// of QAOA operators for use in quantum algorithms, like Adapt-QAOA. +/// @details This class extends the operator_pool interface +/// therefore inherits the extension_point template, allowing for +/// runtime extensibility. class qaoa_pool : public operator_pool { public: + /// @brief Generate a vector of spin operators based on the provided + /// configuration. + /// @param config A heterogeneous map containing configuration parameters for + /// operator generation. + /// @return A vector of cudaq::spin_op objects representing the generated + /// operator pool. std::vector generate(const heterogeneous_map &config) const override; + + /// @brief Call to macro for defining the creator function for the qaoa_pool + /// extension + /// @details This function is used by the extension point mechanism to create + /// instances of the qaoa_pool class. The extension is registered with a name CUDAQ_EXTENSION_CUSTOM_CREATOR_FUNCTION_WITH_NAME( qaoa_pool, "qaoa", static std::unique_ptr create() { return std::make_unique(); }) }; + +/// @brief Register the qaoa_pool extension type with the CUDA-Q framework CUDAQ_REGISTER_TYPE(qaoa_pool) } // namespace cudaq::solvers diff --git a/libs/solvers/include/cudaq/solvers/operators/operator_pools/spin_complement_gsd.h b/libs/solvers/include/cudaq/solvers/operators/operator_pools/spin_complement_gsd.h index 1d2c0919..ad309e48 100644 --- a/libs/solvers/include/cudaq/solvers/operators/operator_pools/spin_complement_gsd.h +++ b/libs/solvers/include/cudaq/solvers/operators/operator_pools/spin_complement_gsd.h @@ -15,13 +15,31 @@ namespace cudaq::solvers { // adapted from // https://github.com/mayhallgroup/adapt-vqe/blob/master/src/operator_pools.py -/// @brief Test +/// @brief spin_complement_gsd operator pool is a class which generates a pool +/// of operators with the spin complement ground state degeneracy method for use +/// in quantum algorithms, like Adapt-VQE. +/// @details This class extends the operator_pool interface +/// therefore inherits extension_point template, allowing for +/// runtime extensibility. class spin_complement_gsd : public operator_pool { public: + /// @brief Generate a vector of spin operators based on the provided + /// configuration. + /// @param config A heterogeneous map containing configuration parameters for + /// operator generation. + /// @return A vector of cudaq::spin_op objects representing the generated + /// operator pool. std::vector generate(const heterogeneous_map &config) const override; + + /// @brief Call to macro for defining the creator function for the + /// spin_complement_gsd extension + /// @details This function is used by the extension point mechanism to create + /// instances of the spin_complement_gsd class. CUDAQ_EXTENSION_CREATOR_FUNCTION(operator_pool, spin_complement_gsd) }; +/// @brief Register the spin_complement_gsd extension type with the CUDA-Q +/// framework CUDAQ_REGISTER_TYPE(spin_complement_gsd) } // namespace cudaq::solvers \ No newline at end of file diff --git a/libs/solvers/include/cudaq/solvers/operators/operator_pools/uccsd_operator_pool.h b/libs/solvers/include/cudaq/solvers/operators/operator_pools/uccsd_operator_pool.h index 02387b1e..c2df7e1d 100644 --- a/libs/solvers/include/cudaq/solvers/operators/operator_pools/uccsd_operator_pool.h +++ b/libs/solvers/include/cudaq/solvers/operators/operator_pools/uccsd_operator_pool.h @@ -12,14 +12,32 @@ namespace cudaq::solvers { -/// @brief Test +/// @brief UCCSD operator pool is a class with generates a pool +/// of UCCSD operators for use in quantum algorithms, like Adapt-VQE. +/// @details This class extends the operator_pool interface +/// therefore inherits the extension_point template, allowing for +/// runtime extensibility. class uccsd : public operator_pool { public: + /// @brief Generate a vector of spin operators based on the provided + /// configuration. + /// @details The UCCSD operator pool is generated with an imaginary factor 'i + /// in the coefficients of the operators, which simplifies the use for running + /// Adapt-VQE. + /// @param config A heterogeneous map containing configuration parameters for + /// operator generation. + /// @return A vector of cudaq::spin_op objects representing the generated + /// operator pool. std::vector generate(const heterogeneous_map &config) const override; + + /// @brief Call to macro for defining the creator function for an extension + /// @details This function is used by the extension point mechanism to create + /// instances of the uccsd class. CUDAQ_EXTENSION_CREATOR_FUNCTION(operator_pool, uccsd) }; +/// @brief Register the uccsd extension type with the CUDA-Q framework CUDAQ_REGISTER_TYPE(uccsd) } // namespace cudaq::solvers diff --git a/libs/solvers/include/cudaq/solvers/vqe.h b/libs/solvers/include/cudaq/solvers/vqe.h index a0c30748..904bc982 100644 --- a/libs/solvers/include/cudaq/solvers/vqe.h +++ b/libs/solvers/include/cudaq/solvers/vqe.h @@ -37,7 +37,8 @@ struct vqe_result { /// @param optimizer Optimization algorithm to use /// @param gradient Gradient computation method /// @param initial_parameters Initial parameters for the optimization -/// @param options Additional options for the VQE algorithm +/// @param options Additional options for the VQE algorithm. Available options +/// - "tol" (double): Tolerance for optimization [default: 1e-12] /// @return VQE result containing optimal energy, parameters, and iteration data template requires std::invocable> diff --git a/libs/solvers/lib/adapt/adapt_simulator.cpp b/libs/solvers/lib/adapt/adapt_simulator.cpp index 45b742e2..3f528b24 100644 --- a/libs/solvers/lib/adapt/adapt_simulator.cpp +++ b/libs/solvers/lib/adapt/adapt_simulator.cpp @@ -6,6 +6,8 @@ * the terms of the Apache License 2.0 which accompanies this distribution. * ******************************************************************************/ +#include + #include "common/Logger.h" #include "cudaq.h" @@ -14,7 +16,6 @@ #include "cudaq/solvers/adapt/adapt_simulator.h" #include "cudaq/solvers/vqe.h" -#include #include namespace cudaq::solvers::adapt { @@ -30,8 +31,18 @@ simulator::run(const cudaq::qkernel &)> &initialState, std::vector pauliWords; std::vector thetas, coefficients; + std::vector poolIndices; std::vector chosenOps; - auto tol = options.get("grad_norm_tolerance", 1e-5); + double latestEnergy = std::numeric_limits::max(); + double ediff = std::numeric_limits::max(); + + int maxIter = options.get("max_iter", 30); + auto grad_norm_tolerance = options.get("grad_norm_tolerance", 1e-5); + auto tolNormDiff = options.get("grad_norm_diff_tolerance", 1e-5); + auto thresholdE = options.get("threshold_energy", 1e-6); + auto initTheta = options.get("initial_theta", 0.0); + auto mutable_options = options; + auto numQubits = H.num_qubits(); // Assumes each rank can see numQpus, models a distributed // architecture where each rank is a compute node, and each node @@ -51,9 +62,39 @@ simulator::run(const cudaq::qkernel &)> &initialState, std::size_t remainder = total_elements % numRanks; std::size_t start = rank * elements_per_rank + std::min(rank, remainder); std::size_t end = start + elements_per_rank + (rank < remainder ? 1 : 0); - for (int i = start; i < end; i++) { - auto op = pool[i]; - commutators.emplace_back(H * op - op * H); + + // Check if operator has only imaginary coefficients + // checking the first one is enough, we assume the pool is homogeneous + const auto &c = pool[0].begin()->get_coefficient(); + bool isImaginary = + (std::abs(c.real()) <= 1e-9) && (std::abs(c.imag()) > 1e-9); + auto coeff = (!isImaginary) ? std::complex{0.0, 1.0} + : std::complex{1.0, 0.0}; + + auto cleanCommutator = + [](cudaq::spin_op commutator) -> std::tuple { + cudaq::spin_op cleaned; + std::size_t numTerms = 0; + commutator.for_each_term([&](const auto &term) { + if (term.get_coefficient().real() != 0.0 || + term.get_coefficient().imag() != 0.0) { + if (numTerms == 0) { + cleaned = term; + numTerms++; + } else { + cleaned += term; + numTerms++; + } + } + }); + return std::make_tuple(cleaned, numTerms); + }; + + for (auto &op : pool) { + auto commutator = H * op - op * H; + auto [cleanedCom, numTerms] = cleanCommutator(commutator); + if (numTerms > 0) + commutators.push_back(coeff * cleanedCom); } nlohmann::json initInfo = {{"num-qpus", numQpus}, @@ -70,28 +111,44 @@ simulator::run(const cudaq::qkernel &)> &initialState, // Start of with the initial |psi_n> cudaq::state state = get_state(adapt_kernel, numQubits, initialState, thetas, - coefficients, pauliWords); - std::size_t count = 0; + coefficients, pauliWords, poolIndices); + + int step = 0; while (true) { + if (options.get("verbose", false)) + printf("Step %d\n", step); + if (step >= maxIter) { + std::cerr + << "Warning: Timed out, number of iteration steps exceeds maxIter!" + << std::endl; + break; + } + step++; // Step 1 - compute vector std::vector gradients; + std::vector results; double gradNorm = 0.0; std::vector resultHandles; - for (std::size_t i = 0, qpuCounter = 0; i < commutators.size(); i++) { - if (rank == 0) - cudaq::info("Compute commutator {}", i); - if (qpuCounter % numQpus == 0) - qpuCounter = 0; - resultHandles.emplace_back( - observe_async(qpuCounter++, prepare_state, commutators[i], state)); + if (numQpus == 1) { + for (std::size_t i = 0; i < commutators.size(); i++) { + cudaq::info("Compute commutator {}", i); + results.emplace_back(observe(prepare_state, commutators[i], state)); + } + } else { + for (std::size_t i = 0, qpuCounter = 0; i < commutators.size(); i++) { + if (rank == 0) + cudaq::info("Compute commutator {}", i); + if (qpuCounter % numQpus == 0) + qpuCounter = 0; + resultHandles.emplace_back( + observe_async(qpuCounter++, prepare_state, commutators[i], state)); + } + for (auto &handle : resultHandles) + results.emplace_back(handle.get()); } - std::vector results; - for (auto &handle : resultHandles) - results.emplace_back(handle.get()); - // Get the gradient results std::transform(results.begin(), results.end(), std::back_inserter(gradients), @@ -101,15 +158,16 @@ simulator::run(const cudaq::qkernel &)> &initialState, double norm = 0.0; for (auto &g : gradients) norm += g * g; + norm = std::sqrt(norm); // All ranks have a norm, need to reduce that across all if (mpi::is_initialized()) norm = cudaq::mpi::all_reduce(norm, std::plus()); - // All ranks have a max gradient and index auto iter = std::max_element(gradients.begin(), gradients.end()); double maxGrad = *iter; auto maxOpIdx = std::distance(gradients.begin(), iter); + if (mpi::is_initialized()) { std::vector allMaxOpIndices(numRanks); std::vector allMaxGrads(numRanks); @@ -140,16 +198,22 @@ simulator::run(const cudaq::qkernel &)> &initialState, } // Convergence is reached if gradient values are small - if (std::sqrt(std::fabs(norm)) < tol || std::fabs(lastNorm - norm) < tol) + if (norm < grad_norm_tolerance || + std::fabs(lastNorm - norm) < tolNormDiff || ediff < thresholdE) break; // Use the operator from the pool auto op = pool[maxOpIdx]; + if (!isImaginary) + op = std::complex{0.0, 1.0} * pool[maxOpIdx]; + chosenOps.push_back(op); - thetas.push_back(0.0); + thetas.push_back(initTheta); + for (auto o : op) { pauliWords.emplace_back(o.to_string(false)); coefficients.push_back(o.get_coefficient().imag()); + poolIndices.push_back(maxOpIdx); } optim::optimizable_function objective; @@ -159,13 +223,9 @@ simulator::run(const cudaq::qkernel &)> &initialState, objective = [&, thetas, coefficients](const std::vector &x, std::vector &dx) mutable { auto res = cudaq::observe(adapt_kernel, H, numQubits, initialState, x, - coefficients, pauliWords); + coefficients, pauliWords, poolIndices); if (options.get("verbose", false)) printf(" = %.12lf\n", res.expectation()); - // data.emplace_back(x, res, observe_execution_type::function); - // for (auto datum : gradient->data) - // data.push_back(datum); - return res.expectation(); }; } else { @@ -178,31 +238,32 @@ simulator::run(const cudaq::qkernel &)> &initialState, [&, thetas, coefficients, pauliWords](const std::vector xx) { std::apply([&](auto &&...new_args) { adapt_kernel(new_args...); }, std::forward_as_tuple(numQubits, initialState, xx, - coefficients, pauliWords)); + coefficients, pauliWords, + poolIndices)); }, H); objective = [&, thetas, coefficients](const std::vector &x, std::vector &dx) mutable { // FIXME get shots in here... auto res = cudaq::observe(adapt_kernel, H, numQubits, initialState, x, - coefficients, pauliWords); + coefficients, pauliWords, poolIndices); if (options.get("verbose", false)) printf(" = %.12lf\n", res.expectation()); defaultGradient->compute(x, dx, res.expectation(), options.get("shots", -1)); - - // data.emplace_back(x, res, observe_execution_type::function); - // for (auto datum : gradient->data) - // data.push_back(datum); - return res.expectation(); }; } - // FIXME fix the const_cast + if (options.contains("dynamic_start")) { + if (options.get("dynamic_start") == "warm") + mutable_options.insert("initial_parameters", thetas); + } + auto [groundEnergy, optParams] = - const_cast(optimizer).optimize(thetas.size(), - objective, options); + const_cast(optimizer).optimize( + thetas.size(), objective, mutable_options); + // Set the new optimzal parameters thetas = optParams; energy = groundEnergy; @@ -210,7 +271,10 @@ simulator::run(const cudaq::qkernel &)> &initialState, // Set the norm for the next iteration's check lastNorm = norm; state = get_state(adapt_kernel, numQubits, initialState, thetas, - coefficients, pauliWords); + coefficients, pauliWords, poolIndices); + + ediff = std::fabs(latestEnergy - groundEnergy); + latestEnergy = groundEnergy; } return std::make_tuple(energy, thetas, chosenOps); diff --git a/libs/solvers/lib/adapt/device/adapt.cpp b/libs/solvers/lib/adapt/device/adapt.cpp index c5895859..28946dc3 100644 --- a/libs/solvers/lib/adapt/device/adapt.cpp +++ b/libs/solvers/lib/adapt/device/adapt.cpp @@ -7,6 +7,7 @@ ******************************************************************************/ #include "cudaq.h" +#include namespace cudaq { @@ -15,13 +16,20 @@ adapt_kernel(std::size_t numQubits, const cudaq::qkernel &)> &statePrep, const std::vector &thetas, const std::vector &coefficients, - const std::vector &trotterOpList) { + const std::vector &trotterOpList, + const std::vector poolIndices) { cudaq::qvector q(numQubits); statePrep(q); - for (std::size_t i = 0; i < thetas.size(); i++) - for (std::size_t j = 0; j < trotterOpList.size(); j++) { - auto &term = trotterOpList[j]; - exp_pauli(thetas[i] * coefficients[j], q, term); + + auto i = 0; + for (std::size_t j = 0; j < trotterOpList.size();) { + auto index = poolIndices[j]; + auto &term = trotterOpList[j]; + exp_pauli(thetas[i] * coefficients[j], q, term); + j++; + if (index != poolIndices[j]) { + i++; } + } } } // namespace cudaq diff --git a/libs/solvers/lib/adapt/device/adapt.h b/libs/solvers/lib/adapt/device/adapt.h index 2363b75b..7d807b3c 100644 --- a/libs/solvers/lib/adapt/device/adapt.h +++ b/libs/solvers/lib/adapt/device/adapt.h @@ -27,6 +27,8 @@ namespace cudaq { /// @param coefficients Vector of coefficients for the Hamiltonian terms /// @param trotterOpList Vector of Pauli words representing the Trotter /// operators +/// @param poolIndices Vector of indices indicating which terms belong to which +/// operator /// /// @note This is a CUDA-Q quantum kernel function and should be executed within /// the CUDA-Q framework. It applies the ADAPT-VQE circuit construction based on @@ -37,6 +39,7 @@ void adapt_kernel(std::size_t numQubits, const cudaq::qkernel &)> &statePrep, const std::vector &thetas, const std::vector &coefficients, - const std::vector &trotterOpList); + const std::vector &trotterOpList, + const std::vector poolIndices); } // namespace cudaq diff --git a/libs/solvers/lib/operators/operator_pools/uccsd_operator_pool.cpp b/libs/solvers/lib/operators/operator_pools/uccsd_operator_pool.cpp index 2f90b438..2f4690cf 100644 --- a/libs/solvers/lib/operators/operator_pools/uccsd_operator_pool.cpp +++ b/libs/solvers/lib/operators/operator_pools/uccsd_operator_pool.cpp @@ -5,7 +5,6 @@ * This source code and the accompanying materials are made available under * * the terms of the Apache License 2.0 which accompanies this distribution. * ******************************************************************************/ - #include "cudaq/solvers/operators/operator_pools/uccsd_operator_pool.h" #include "cudaq/solvers/stateprep/uccsd.h" @@ -38,9 +37,9 @@ uccsd::generate(const heterogeneous_map &config) const { cudaq::spin_op o(numQubits); for (std::size_t i = p + 1; i < q; i++) o *= cudaq::spin::z(i); - - ops.emplace_back(cudaq::spin::y(p) * o * cudaq::spin::x(q)); - ops.emplace_back(cudaq::spin::x(p) * o * cudaq::spin::y(q)); + std::complex c = {0.5, 0}; + ops.emplace_back(c * cudaq::spin::y(p) * o * cudaq::spin::x(q) - + c * cudaq::spin::x(p) * o * cudaq::spin::y(q)); }; auto addDoublesExcitation = [numQubits](std::vector &ops, @@ -65,9 +64,7 @@ uccsd::generate(const heterogeneous_map &config) const { j_occ = q; a_virt = s; b_virt = r; - } else if - - (p > q && r < s) { + } else if (p > q && r < s) { i_occ = q; j_occ = p; a_virt = r; @@ -79,30 +76,26 @@ uccsd::generate(const heterogeneous_map &config) const { for (std::size_t i = a_virt + 1; i < b_virt; i++) parity_b *= cudaq::spin::z(i); - ops.emplace_back(cudaq::spin::x(i_occ) * parity_a * cudaq::spin::x(j_occ) * - cudaq::spin::x(a_virt) * parity_b * - cudaq::spin::y(b_virt)); - ops.emplace_back(cudaq::spin::x(i_occ) * parity_a * cudaq::spin::x(j_occ) * - cudaq::spin::y(a_virt) * parity_b * - cudaq::spin::x(b_virt)); - ops.emplace_back(cudaq::spin::x(i_occ) * parity_a * cudaq::spin::y(j_occ) * - cudaq::spin::y(a_virt) * parity_b * - cudaq::spin::y(b_virt)); - ops.emplace_back(cudaq::spin::y(i_occ) * parity_a * cudaq::spin::x(j_occ) * - cudaq::spin::y(a_virt) * parity_b * - cudaq::spin::y(b_virt)); - ops.emplace_back(cudaq::spin::x(i_occ) * parity_a * cudaq::spin::y(j_occ) * - cudaq::spin::x(a_virt) * parity_b * - cudaq::spin::x(b_virt)); - ops.emplace_back(cudaq::spin::y(i_occ) * parity_a * cudaq::spin::x(j_occ) * - cudaq::spin::x(a_virt) * parity_b * - cudaq::spin::x(b_virt)); - ops.emplace_back(cudaq::spin::y(i_occ) * parity_a * cudaq::spin::y(j_occ) * - cudaq::spin::x(a_virt) * parity_b * - cudaq::spin::y(b_virt)); - ops.emplace_back(cudaq::spin::y(i_occ) * parity_a * cudaq::spin::y(j_occ) * - cudaq::spin::y(a_virt) * parity_b * - cudaq::spin::x(b_virt)); + auto op_term_temp = cudaq::spin::x(i_occ) * parity_a * + cudaq::spin::x(j_occ) * cudaq::spin::x(a_virt) * + parity_b * cudaq::spin::y(b_virt); + op_term_temp += cudaq::spin::x(i_occ) * parity_a * cudaq::spin::x(j_occ) * + cudaq::spin::y(a_virt) * parity_b * cudaq::spin::x(b_virt); + op_term_temp += cudaq::spin::x(i_occ) * parity_a * cudaq::spin::y(j_occ) * + cudaq::spin::y(a_virt) * parity_b * cudaq::spin::y(b_virt); + op_term_temp += cudaq::spin::y(i_occ) * parity_a * cudaq::spin::x(j_occ) * + cudaq::spin::y(a_virt) * parity_b * cudaq::spin::y(b_virt); + op_term_temp -= cudaq::spin::x(i_occ) * parity_a * cudaq::spin::y(j_occ) * + cudaq::spin::x(a_virt) * parity_b * cudaq::spin::x(b_virt); + op_term_temp -= cudaq::spin::y(i_occ) * parity_a * cudaq::spin::x(j_occ) * + cudaq::spin::x(a_virt) * parity_b * cudaq::spin::x(b_virt); + op_term_temp -= cudaq::spin::y(i_occ) * parity_a * cudaq::spin::y(j_occ) * + cudaq::spin::x(a_virt) * parity_b * cudaq::spin::y(b_virt); + op_term_temp -= cudaq::spin::y(i_occ) * parity_a * cudaq::spin::y(j_occ) * + cudaq::spin::y(a_virt) * parity_b * cudaq::spin::x(b_virt); + + std::complex c = {0.125, 0}; + ops.emplace_back(c * op_term_temp); }; for (auto &sa : singlesAlpha) diff --git a/libs/solvers/lib/optimizers/lbfgs/lbfgs.cpp b/libs/solvers/lib/optimizers/lbfgs/lbfgs.cpp index 52bca68a..f6573d35 100644 --- a/libs/solvers/lib/optimizers/lbfgs/lbfgs.cpp +++ b/libs/solvers/lib/optimizers/lbfgs/lbfgs.cpp @@ -16,7 +16,7 @@ optimization_result lbfgs::optimize(std::size_t dim, history.clear(); cudaq::optim::LBFGSObjective f( opt_function, options.get("initial_parameters", std::vector(dim)), - options.get("function_tolerance", 1e-12), + options.get("tol", 1e-12), options.get("max_iterations", std::numeric_limits::max()), history, options.get("verbose", false)); return f.run(dim); diff --git a/libs/solvers/python/bindings/solvers/py_solvers.cpp b/libs/solvers/python/bindings/solvers/py_solvers.cpp index 67c45d5a..2e9a2d6f 100644 --- a/libs/solvers/python/bindings/solvers/py_solvers.cpp +++ b/libs/solvers/python/bindings/solvers/py_solvers.cpp @@ -807,6 +807,8 @@ The exact set of possible types may depend on the specific optimization algorith optOptions.insert( "max_iterations", cudaqx::getValueOr(options, "max_iterations", -1)); + // in case the privded optimizer is not a scipy one + optOptions.insert("tol", getValueOr(options, "tol", 1e-12)); optOptions.insert("verbose", cudaqx::getValueOr(options, "verbose", false)); @@ -876,6 +878,7 @@ options : dict - verbose : bool, optional Whether to print verbose output. Default is False. - optimizer : str, optional Name of the classical optimizer to use. Default is 'cobyla'. - gradient : str, optional Method for gradient computation (for gradient-based optimizers). Default is 'parameter_shift'. + - tol (double): Tolerance value for the optimizer. Default 1e-12. Returns: -------- @@ -930,8 +933,25 @@ RuntimeError auto *p = reinterpret_cast(fptr); cudaq::registry::__cudaq_registerLinkableKernel(p, baseName.c_str(), p); heterogeneous_map optOptions; + optOptions.insert("max_iter", getValueOr(options, "max_iter", 30)); + optOptions.insert( + "grad_norm_tolerance", + getValueOr(options, "grad_norm_tolerance", 1e-5)); + optOptions.insert( + "grad_norm_diff_tolerance", + getValueOr(options, "grad_norm_diff_tolerance", 1e-5)); + optOptions.insert( + "threshold_energy", + getValueOr(options, "threshold_energy", 1e-6)); + optOptions.insert("initial_theta", + getValueOr(options, "initial_theta", 0.0)); optOptions.insert("verbose", getValueOr(options, "verbose", false)); + optOptions.insert("shots", getValueOr(options, "shots", -1)); + optOptions.insert("tol", getValueOr(options, "tol", 1e-12)); + optOptions.insert( + "dynamic_start", + getValueOr(options, "dynamic_start", "cold")); // Handle the case where the user has provided a SciPy optimizer if (options.contains("optimizer") && @@ -967,6 +987,19 @@ RuntimeError Keyword Args: optimizer (str): Optional name of the optimizer to use. Defaults to cobyla. gradient (str): Optional name of the gradient method to use. Defaults to empty. + + Options Dictionary: + The following keys are supported in the options dictionary: + - max_iter (int): Maximum number of iterations. Default: 30 + - grad_norm_tolerance (float): Convergence tolerance for gradient norm. Default: 1e-5 + - grad_norm_diff_tolerance (float): Tolerance for difference between gradient norms. Default: 1e-5 + - threshold_energy (float): Energy convergence threshold. Default: 1e-6 + - initial_theta (float): Initial value for theta parameter. Default: 0.0 + - verbose (bool): Enable detailed output logging. Default: False + - shots (int): Number of measurement shots (-1 for exact simulation). Default: -1 + - tol (double): Tolerance value for the optimizer. Default 1e-12 + - dynamic_start (string): Optimization mode for the theta parameters at each iteration. It can be either "warm", or "cold". Default: "cold" + Returns: The result of the ADAPT-VQE optimization. diff --git a/libs/solvers/python/tests/test_adapt.py b/libs/solvers/python/tests/test_adapt.py index 9a7da027..e639da1a 100644 --- a/libs/solvers/python/tests/test_adapt.py +++ b/libs/solvers/python/tests/test_adapt.py @@ -33,11 +33,23 @@ def initState(q: cudaq.qview): print(energy) assert np.isclose(energy, -1.137, atol=1e-3) - energy, thetas, ops = solvers.adapt_vqe(initState, - molecule.hamiltonian, - operators, - optimizer='lbfgs', - gradient='central_difference') + energy, thetas, ops = solvers.adapt_vqe( + initState, + molecule.hamiltonian, + operators, + optimizer='lbfgs', + gradient='central_difference', + options={ + 'max_iter': 5, # Maximum iterations + 'grad_norm_tolerance': 1e-6, # Gradient norm convergence tolerance + 'grad_norm_diff_tolerance': + 1e-6, # Gradient norm difference tolerance + 'threshold_energy': 1e-7, # Energy convergence threshold + 'initial_theta': 0.1, # Initial theta parameter + 'verbose': True, # Enable detailed logging + 'shots': 10000 # Number of measurement shots + }) + print(energy) assert np.isclose(energy, -1.137, atol=1e-3) diff --git a/libs/solvers/python/tests/test_operator_pools.py b/libs/solvers/python/tests/test_operator_pools.py index f7b5570c..3bf70c02 100644 --- a/libs/solvers/python/tests/test_operator_pools.py +++ b/libs/solvers/python/tests/test_operator_pools.py @@ -16,7 +16,7 @@ def test_generate_with_default_config(): num_qubits=4, num_electrons=2) assert operators - assert len(operators) == 2 * 2 + 1 * 8 + assert len(operators) == 2 + 1 for op in operators: assert op.get_qubit_count() == 4 @@ -28,12 +28,13 @@ def test_generate_with_custom_coefficients(): num_electrons=2) assert operators - assert len(operators) == (2 * 2 + 1 * 8) + assert len(operators) == (2 + 1) for i, op in enumerate(operators): assert op.get_qubit_count() == 4 - expected_coeff = 1.0 - assert np.isclose(op.get_coefficient().real, expected_coeff) + expected_coeff = [0.5, 0.125] + for term in op: + assert (abs(term.get_coefficient().real) in expected_coeff) def test_generate_with_odd_electrons(): @@ -43,7 +44,7 @@ def test_generate_with_odd_electrons(): spin=1) assert operators - assert len(operators) == 2 * 4 + 4 * 8 + assert len(operators) == 2 * 2 + 4 for op in operators: assert op.get_qubit_count() == 6 @@ -55,47 +56,60 @@ def test_generate_with_large_system(): num_electrons=10) assert operators - assert len(operators) > 875 + assert len(operators) == 875 for op in operators: assert op.get_qubit_count() == 20 def test_uccsd_operator_pool_correctness(): - # Generate the UCCSD operator pool pool = solvers.get_operator_pool("uccsd", num_qubits=4, num_electrons=2) - # Convert SpinOperators to strings - pool_strings = [op.to_string(False) for op in pool] + temp_data = [[], [], []] + data_counter = 0 + for op in pool: + op.for_each_term(lambda term: temp_data[data_counter].append( + (term.to_string(False), term.get_coefficient()))) + data_counter += 1 + + # Assert + expected_operators = [["XZYI", "YZXI"], ["IXZY", "IYZX"], + [ + "YYYX", "YXXX", "XXYX", "YYXY", "XYYY", "XXXY", + "YXYY", "XYXX" + ]] + expected_coefficients = [[complex(-0.5, 0), + complex(0.5, 0)], + [complex(-0.5, 0), + complex(0.5, 0)], + [ + complex(-0.125, 0), + complex(-0.125, 0), + complex(0.125, 0), + complex(-0.125, 0), + complex(0.125, 0), + complex(0.125, 0), + complex(0.125, 0), + complex(-0.125, 0) + ]] - # Expected result - expected_pool = [ - 'YZXI', 'XZYI', 'IYZX', 'IXZY', 'XXXY', 'XXYX', 'XYYY', 'YXYY', 'XYXX', - 'YXXX', 'YYXY', 'YYYX' - ] - - # Assert that the generated pool matches the expected result - assert pool_strings == expected_pool, f"Expected {expected_pool}, but got {pool_strings}" - - # Additional checks - assert len(pool) == len( - expected_pool - ), f"Expected {len(expected_pool)} operators, but got {len(pool)}" - - # Check that all operators have the correct length (4 qubits) - for op_string in pool_strings: - assert len( - op_string - ) == 4, f"Operator {op_string} does not have the expected length of 4" - - # Check that all operators contain only valid characters (I, X, Y, Z) valid_chars = set('IXYZ') - for op_string in pool_strings: - assert set(op_string).issubset( - valid_chars), f"Operator {op_string} contains invalid characters" - - -def test_generate_with_invalid_config(): - # Missing required parameters - with pytest.raises(RuntimeError): - pool = solvers.get_operator_pool("uccsd") + assert len(temp_data) == len( + expected_operators + ), f"Number of generated operators ({len(temp_data)}) does not match expected count ({len(expected_operators)})" + + for i in range(len(temp_data)): + for j in range(len(temp_data[i])): + op_string = temp_data[i][j][0] + op_coeff = temp_data[i][j][1] + # Check operator length + assert len( + op_string + ) == 4, f"Operator {op_string} does not have the expected length of 4" + index = expected_operators[i].index(op_string) + + assert op_coeff == expected_coefficients[i][index], \ + f"Coefficient mismatch at index {i}, {index}: expected {expected_coefficients[i][index]}, got {op_coeff}" + + assert set(op_string).issubset(valid_chars), \ + f"Operator {op_string} contains invalid characters" diff --git a/libs/solvers/python/tests/test_vqe.py b/libs/solvers/python/tests/test_vqe.py index 8b86c3c8..e1b6cac6 100644 --- a/libs/solvers/python/tests/test_vqe.py +++ b/libs/solvers/python/tests/test_vqe.py @@ -29,11 +29,12 @@ def ansatz(theta: float): hamiltonian = 5.907 - 2.1433 * spin.x(0) * spin.x(1) - 2.1433 * spin.y( 0) * spin.y(1) + .21829 * spin.z(0) - 6.125 * spin.z(1) - # Can specify optimizer and gradient + # Can specify optimizer and gradient and tol energy, params, all_data = solvers.vqe(lambda thetas: ansatz(thetas[0]), hamiltonian, [0.], optimizer='lbfgs', - gradient='parameter_shift') + gradient='parameter_shift', + tol=1e-7) assert np.isclose(-1.74, energy, atol=1e-2) all_data[0].result.dump() diff --git a/libs/solvers/unittests/nvqpp/test_kernels.cpp b/libs/solvers/unittests/nvqpp/test_kernels.cpp index 255d3685..10a60384 100644 --- a/libs/solvers/unittests/nvqpp/test_kernels.cpp +++ b/libs/solvers/unittests/nvqpp/test_kernels.cpp @@ -7,6 +7,22 @@ __qpu__ void hartreeFock2Electrons(cudaq::qvector<> &q) { x(q[i]); } +__qpu__ void statePrepNElectrons(cudaq::qvector<> &q, + std::size_t numElectrons) { + for (std::size_t i = 0; i < numElectrons; i++) + x(q[i]); +} + +__qpu__ void statePrep4Electrons(cudaq::qvector<> &q) { + for (std::size_t i = 0; i < 4; i++) + x(q[i]); +} + +__qpu__ void statePrep6Electrons(cudaq::qvector<> &q) { + for (std::size_t i = 0; i < 6; i++) + x(q[i]); +} + __qpu__ void ansatz(std::vector theta) { cudaq::qvector q(2); x(q[0]); diff --git a/libs/solvers/unittests/nvqpp/test_kernels.h b/libs/solvers/unittests/nvqpp/test_kernels.h index b4ccd5b5..e21ad514 100644 --- a/libs/solvers/unittests/nvqpp/test_kernels.h +++ b/libs/solvers/unittests/nvqpp/test_kernels.h @@ -1,6 +1,9 @@ #include "cudaq.h" __qpu__ void hartreeFock2Electrons(cudaq::qvector<> &q); +__qpu__ void statePrepNElectrons(cudaq::qvector<> &q, std::size_t numElectrons); +__qpu__ void statePrep4Electrons(cudaq::qvector<> &q); +__qpu__ void statePrep6Electrons(cudaq::qvector<> &q); __qpu__ void ansatz(std::vector theta); __qpu__ void ansatzNonStdSignature(double theta, int N); __qpu__ void callUccsdStatePrep(std::vector params); diff --git a/libs/solvers/unittests/test_adapt.cpp b/libs/solvers/unittests/test_adapt.cpp index f04baf07..b7309c29 100644 --- a/libs/solvers/unittests/test_adapt.cpp +++ b/libs/solvers/unittests/test_adapt.cpp @@ -1,5 +1,5 @@ /******************************************************************************* - * Copyright (c) 2024 NVIDIA Corporation & Affiliates. * + * Copyright (c) 2024 - 2025 NVIDIA Corporation & Affiliates. * * All rights reserved. * * * * This source code and the accompanying materials are made available under * @@ -12,37 +12,483 @@ #include "cudaq.h" #include "nvqpp/test_kernels.h" #include "cudaq/solvers/adapt.h" +#include "cudaq/solvers/operators.h" -std::vector h2_data{ - 3, 1, 1, 3, 0.0454063, 0, 2, 0, 0, 0, 0.17028, 0, - 0, 0, 2, 0, -0.220041, -0, 1, 3, 3, 1, 0.0454063, 0, - 0, 0, 0, 0, -0.106477, 0, 0, 2, 0, 0, 0.17028, 0, - 0, 0, 0, 2, -0.220041, -0, 3, 3, 1, 1, -0.0454063, -0, - 2, 2, 0, 0, 0.168336, 0, 2, 0, 2, 0, 0.1202, 0, - 0, 2, 0, 2, 0.1202, 0, 2, 0, 0, 2, 0.165607, 0, - 0, 2, 2, 0, 0.165607, 0, 0, 0, 2, 2, 0.174073, 0, - 1, 1, 3, 3, -0.0454063, -0, 15}; +bool check_gpu_available() { + static int gpu_available = -1; // -1 means unknown; 0 = no; 1 = yes + if (gpu_available < 0) { + if (std::system("nvidia-smi > /dev/null 2>&1") != 0) + gpu_available = 0; + else + gpu_available = 1; + } + return static_cast(gpu_available); +} + +class SolversTester : public ::testing::Test { +protected: + static cudaq::spin_op h; + static cudaq::spin_op hamli; + static cudaq::spin_op hamhh; + static cudaq::spin_op hambeh2; + + static void SetUpTestSuite() { + std::vector h2_data{ + 3, 1, 1, 3, 0.0454063, 0, 2, 0, 0, 0, 0.17028, 0, + 0, 0, 2, 0, -0.220041, -0, 1, 3, 3, 1, 0.0454063, 0, + 0, 0, 0, 0, -0.106477, 0, 0, 2, 0, 0, 0.17028, 0, + 0, 0, 0, 2, -0.220041, -0, 3, 3, 1, 1, -0.0454063, -0, + 2, 2, 0, 0, 0.168336, 0, 2, 0, 2, 0, 0.1202, 0, + 0, 2, 0, 2, 0.1202, 0, 2, 0, 0, 2, 0.165607, 0, + 0, 2, 2, 0, 0.165607, 0, 0, 0, 2, 2, 0.174073, 0, + 1, 1, 3, 3, -0.0454063, -0, 15}; + h = cudaq::spin_op(h2_data, 4); + + cudaq::solvers::molecular_geometry geometryLiH = { + {"Li", {0.3925, 0., 0.}}, {"H", {-1.1774, 0., 0.0}}}; + cudaq::solvers::molecular_geometry geometryHH = {{"H", {0., 0., 0.}}, + {"H", {0., 0., .7474}}}; + cudaq::solvers::molecular_geometry geometryBeH2 = {{"Be", {0.0, 0.0, 0.0}}, + {"H", {0.0, 0.0, 1.3}}, + {"H", {0.0, 0.0, -1.3}}}; + auto hh = cudaq::solvers::create_molecule( + geometryHH, "sto-3g", 0, 0, + {.casci = true, .ccsd = true, .verbose = true}); + auto lih = cudaq::solvers::create_molecule( + geometryLiH, "sto-3g", 0, 0, + {.casci = true, .ccsd = true, .verbose = true}); + auto beh2 = cudaq::solvers::create_molecule( + geometryBeH2, "sto-3g", 0, 0, + {.casci = true, .ccsd = true, .verbose = true}); -TEST(SolversTester, checkSimpleAdapt) { - cudaq::spin_op h(h2_data, 4); + hamli = lih.hamiltonian; + hamhh = hh.hamiltonian; + hambeh2 = beh2.hamiltonian; + } +}; + +cudaq::spin_op SolversTester::h; +cudaq::spin_op SolversTester::hamli; +cudaq::spin_op SolversTester::hamhh; +cudaq::spin_op SolversTester::hambeh2; + +TEST_F(SolversTester, checkSimpleAdapt_H2) { auto pool = cudaq::solvers::operator_pool::get("spin_complement_gsd"); auto poolList = pool->generate({{"num-orbitals", h.num_qubits() / 2}}); - auto [energy, thetas, ops] = cudaq::solvers::adapt_vqe( - hartreeFock2Electrons, h, poolList, - {{"grad_norm_tolerance", 1e-3}, {"verbose", true}}); + auto [energy, thetas, ops] = + cudaq::solvers::adapt_vqe(hartreeFock2Electrons, h, poolList, + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-6}}); EXPECT_NEAR(energy, -1.13, 1e-2); } -TEST(SolversTester, checkSimpleAdaptGradient) { - cudaq::spin_op h(h2_data, 4); +TEST_F(SolversTester, checkSimpleAdapt_H2Sto3g) { auto pool = cudaq::solvers::operator_pool::get("spin_complement_gsd"); + auto poolList = pool->generate({{"num-orbitals", hamhh.num_qubits() / 2}}); + auto [energy, thetas, ops] = + cudaq::solvers::adapt_vqe(hartreeFock2Electrons, hamhh, poolList, + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-6}}); + EXPECT_NEAR(energy, -1.13, 1e-2); +} + +TEST_F(SolversTester, checkSimpleAdaptGradient_H2) { + auto pool = cudaq::solvers::operator_pool::get("spin_complement_gsd"); + auto opt = cudaq::optim::optimizer::get("lbfgs"); auto poolList = pool->generate({{"num-orbitals", h.num_qubits() / 2}}); + auto [energy, thetas, ops] = cudaq::solvers::adapt_vqe( + hartreeFock2Electrons, h, poolList, *opt, "central_difference", + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-6}}); + EXPECT_NEAR(energy, -1.13, 1e-2); + for (std::size_t i = 0; i < thetas.size(); i++) + printf("%lf -> %s\n", thetas[i], ops[i].to_string().c_str()); +} + +TEST_F(SolversTester, checkSimpleAdaptGradient_H2Sto3g) { + auto pool = cudaq::solvers::operator_pool::get("spin_complement_gsd"); auto opt = cudaq::optim::optimizer::get("lbfgs"); + auto poolList = pool->generate({{"num-orbitals", hamhh.num_qubits() / 2}}); + auto [energy, thetas, ops] = cudaq::solvers::adapt_vqe( + hartreeFock2Electrons, hamhh, poolList, *opt, "central_difference", + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-6}}); + EXPECT_NEAR(energy, -1.13, 1e-2); + for (std::size_t i = 0; i < thetas.size(); i++) + printf("%lf -> %s\n", thetas[i], ops[i].to_string().c_str()); +} + +TEST_F(SolversTester, checkSimpleAdaptUCCSD_H2) { + auto pool = cudaq::solvers::operator_pool::get("uccsd"); + heterogeneous_map config; + config.insert("num-qubits", h.num_qubits()); + config.insert("num-electrons", 2); + auto poolList = pool->generate(config); + auto [energy, thetas, ops] = + cudaq::solvers::adapt_vqe(hartreeFock2Electrons, h, poolList, + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-6}}); + EXPECT_NEAR(energy, -1.13, 1e-2); +} + +TEST_F(SolversTester, checkSimpleAdaptUCCSD_H2Sto3g) { + auto pool = cudaq::solvers::operator_pool::get("uccsd"); + heterogeneous_map config; + config.insert("num-qubits", hamhh.num_qubits()); + config.insert("num-electrons", 2); + auto poolList = pool->generate(config); + auto [energy, thetas, ops] = + cudaq::solvers::adapt_vqe(hartreeFock2Electrons, hamhh, poolList, + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-6}}); + EXPECT_NEAR(energy, -1.13, 1e-2); +} + +TEST_F(SolversTester, checkSimpleAdaptGradientUCCSD_H2) { + auto pool = cudaq::solvers::operator_pool::get("uccsd"); + auto opt = cudaq::optim::optimizer::get("lbfgs"); + heterogeneous_map config; + config.insert("num-qubits", h.num_qubits()); + config.insert("num-electrons", 2); + auto poolList = pool->generate(config); + auto [energy, thetas, ops] = cudaq::solvers::adapt_vqe( + hartreeFock2Electrons, h, poolList, *opt, "central_difference", + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-6}}); + EXPECT_NEAR(energy, -1.13, 1e-2); + + for (std::size_t i = 0; i < thetas.size(); i++) + printf("%lf -> %s\n", thetas[i], ops[i].to_string().c_str()); +} + +TEST_F(SolversTester, checkSimpleAdaptGradientUCCSD_H2Sto3g) { + auto pool = cudaq::solvers::operator_pool::get("uccsd"); + auto opt = cudaq::optim::optimizer::get("lbfgs"); + heterogeneous_map config; + config.insert("num-qubits", hamhh.num_qubits()); + config.insert("num-electrons", 2); + auto poolList = pool->generate(config); + auto [energy, thetas, ops] = cudaq::solvers::adapt_vqe( + hartreeFock2Electrons, hamhh, poolList, *opt, "central_difference", + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-6}}); + EXPECT_NEAR(energy, -1.13, 1e-2); + for (std::size_t i = 0; i < thetas.size(); i++) + printf("%lf -> %s\n", thetas[i], ops[i].to_string().c_str()); +} + +TEST_F(SolversTester, checkSimpleAdaptUCCSD_H2_warm) { + auto pool = cudaq::solvers::operator_pool::get("uccsd"); + heterogeneous_map config; + config.insert("num-qubits", h.num_qubits()); + config.insert("num-electrons", 2); + auto poolList = pool->generate(config); + auto [energy, thetas, ops] = + cudaq::solvers::adapt_vqe(hartreeFock2Electrons, h, poolList, + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-6}, + {"dynamic_start", "warm"}}); + EXPECT_NEAR(energy, -1.13, 1e-2); +} + +TEST_F(SolversTester, checkSimpleAdaptUCCSD_H2Sto3g_warm) { + auto pool = cudaq::solvers::operator_pool::get("uccsd"); + heterogeneous_map config; + config.insert("num-qubits", hamhh.num_qubits()); + config.insert("num-electrons", 2); + auto poolList = pool->generate(config); + auto [energy, thetas, ops] = + cudaq::solvers::adapt_vqe(hartreeFock2Electrons, hamhh, poolList, + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-6}, + {"dynamic_start", "warm"}}); + EXPECT_NEAR(energy, -1.13, 1e-2); +} + +TEST_F(SolversTester, checkSimpleAdaptGradientUCCSD_H2_warm) { + auto pool = cudaq::solvers::operator_pool::get("uccsd"); + auto opt = cudaq::optim::optimizer::get("lbfgs"); + heterogeneous_map config; + config.insert("num-qubits", h.num_qubits()); + config.insert("num-electrons", 2); + auto poolList = pool->generate(config); auto [energy, thetas, ops] = cudaq::solvers::adapt_vqe( hartreeFock2Electrons, h, poolList, *opt, "central_difference", - {{"grad_norm_tolerance", 1e-3}, {"verbose", true}}); + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-6}, + {"dynamic_start", "warm"}}); EXPECT_NEAR(energy, -1.13, 1e-2); + for (std::size_t i = 0; i < thetas.size(); i++) + printf("%lf -> %s\n", thetas[i], ops[i].to_string().c_str()); +} + +TEST_F(SolversTester, checkSimpleAdaptGradientUCCSD_H2Sto3g_warm) { + auto pool = cudaq::solvers::operator_pool::get("uccsd"); + auto opt = cudaq::optim::optimizer::get("lbfgs"); + heterogeneous_map config; + config.insert("num-qubits", hamhh.num_qubits()); + config.insert("num-electrons", 2); + auto poolList = pool->generate(config); + auto [energy, thetas, ops] = cudaq::solvers::adapt_vqe( + hartreeFock2Electrons, hamhh, poolList, *opt, "central_difference", + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-6}, + {"dynamic_start", "warm"}}); + EXPECT_NEAR(energy, -1.13, 1e-2); + for (std::size_t i = 0; i < thetas.size(); i++) + printf("%lf -> %s\n", thetas[i], ops[i].to_string().c_str()); +} + +TEST_F(SolversTester, checkSimpleAdapt_LiHSto3g) { + if (!check_gpu_available()) + GTEST_SKIP() << "No GPU available, skipping test because CPU is slow"; + + auto pool = cudaq::solvers::operator_pool::get("spin_complement_gsd"); + auto poolList = pool->generate({{"num-orbitals", hamli.num_qubits() / 2}}); + auto [energy, thetas, ops] = + cudaq::solvers::adapt_vqe(statePrep4Electrons, hamli, poolList, + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-5}}); + EXPECT_NEAR(energy, -7.88, 1e-2); +} + +TEST_F(SolversTester, checkSimpleAdaptGradient_LiHSto3g) { + if (!check_gpu_available()) + GTEST_SKIP() << "No GPU available, skipping test because CPU is slow"; + + auto pool = cudaq::solvers::operator_pool::get("spin_complement_gsd"); + auto opt = cudaq::optim::optimizer::get("lbfgs"); + auto poolList = pool->generate({{"num-orbitals", hamli.num_qubits() / 2}}); + auto [energy, thetas, ops] = cudaq::solvers::adapt_vqe( + statePrep4Electrons, hamli, poolList, *opt, "central_difference", + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-5}}); + EXPECT_NEAR(energy, -7.88, 1e-2); + for (std::size_t i = 0; i < thetas.size(); i++) + printf("%lf -> %s\n", thetas[i], ops[i].to_string().c_str()); +} + +TEST_F(SolversTester, checkSimpleAdaptUCCSD_LiHSto3g) { + if (!check_gpu_available()) + GTEST_SKIP() << "No GPU available, skipping test because CPU is slow"; + + auto pool = cudaq::solvers::operator_pool::get("uccsd"); + heterogeneous_map config; + config.insert("num-qubits", hamli.num_qubits()); + config.insert("num-electrons", 4); + auto poolList = pool->generate(config); + auto [energy, thetas, ops] = + cudaq::solvers::adapt_vqe(statePrep4Electrons, hamli, poolList, + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-5}}); + EXPECT_NEAR(energy, -7.88, 1e-2); +} + +TEST_F(SolversTester, checkSimpleAdaptGradientUCCSD_LiHSto3g) { + if (!check_gpu_available()) + GTEST_SKIP() << "No GPU available, skipping test because CPU is slow"; + + auto pool = cudaq::solvers::operator_pool::get("uccsd"); + auto opt = cudaq::optim::optimizer::get("lbfgs"); + heterogeneous_map config; + config.insert("num-qubits", hamli.num_qubits()); + config.insert("num-electrons", 4); + auto poolList = pool->generate(config); + auto [energy, thetas, ops] = cudaq::solvers::adapt_vqe( + statePrep4Electrons, hamli, poolList, *opt, "central_difference", + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-5}}); + EXPECT_NEAR(energy, -7.88, 1e-2); + for (std::size_t i = 0; i < thetas.size(); i++) + printf("%lf -> %s\n", thetas[i], ops[i].to_string().c_str()); +} + +TEST_F(SolversTester, checkSimpleAdaptGradientUCCSD_LiHSto3g_warm) { + if (!check_gpu_available()) + GTEST_SKIP() << "No GPU available, skipping test because CPU is slow"; + + auto pool = cudaq::solvers::operator_pool::get("uccsd"); + auto opt = cudaq::optim::optimizer::get("lbfgs"); + heterogeneous_map config; + config.insert("num-qubits", hamli.num_qubits()); + config.insert("num-electrons", 4); + auto poolList = pool->generate(config); + auto [energy, thetas, ops] = cudaq::solvers::adapt_vqe( + statePrep4Electrons, hamli, poolList, *opt, "central_difference", + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-5}, + {"dynamic_start", "warm"}}); + EXPECT_NEAR(energy, -7.88, 1e-2); + for (std::size_t i = 0; i < thetas.size(); i++) + printf("%lf -> %s\n", thetas[i], ops[i].to_string().c_str()); +} + +TEST_F(SolversTester, checkSimpleAdapt_BeH2Sto3g) { + if (!check_gpu_available()) + GTEST_SKIP() << "No GPU available, skipping test because CPU is slow"; + + auto pool = cudaq::solvers::operator_pool::get("spin_complement_gsd"); + auto poolList = pool->generate({{"num-orbitals", hambeh2.num_qubits() / 2}}); + auto [energy, thetas, ops] = + cudaq::solvers::adapt_vqe(statePrep6Electrons, hambeh2, poolList, + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-5}}); + EXPECT_NEAR(energy, -15.59, 1e-2); +} + +TEST_F(SolversTester, checkSimpleAdaptGradient_BeH2Sto3g) { + if (!check_gpu_available()) + GTEST_SKIP() << "No GPU available, skipping test because CPU is slow"; + + auto pool = cudaq::solvers::operator_pool::get("spin_complement_gsd"); + auto opt = cudaq::optim::optimizer::get("lbfgs"); + auto poolList = pool->generate({{"num-orbitals", hambeh2.num_qubits() / 2}}); + auto [energy, thetas, ops] = cudaq::solvers::adapt_vqe( + statePrep6Electrons, hambeh2, poolList, *opt, "central_difference", + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-5}}); + EXPECT_NEAR(energy, -15.59, 1e-2); + for (std::size_t i = 0; i < thetas.size(); i++) + printf("%lf -> %s\n", thetas[i], ops[i].to_string().c_str()); +} + +TEST_F(SolversTester, checkSimpleAdaptUCCSD_BeH2Sto3g) { + if (!check_gpu_available()) + GTEST_SKIP() << "No GPU available, skipping test because CPU is slow"; + + auto pool = cudaq::solvers::operator_pool::get("uccsd"); + heterogeneous_map config; + config.insert("num-qubits", hambeh2.num_qubits()); + config.insert("num-electrons", 6); + auto poolList = pool->generate(config); + auto [energy, thetas, ops] = + cudaq::solvers::adapt_vqe(statePrep6Electrons, hambeh2, poolList, + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-5}}); + EXPECT_NEAR(energy, -15.59, 1e-2); +} + +TEST_F(SolversTester, checkSimpleAdaptGradientUCCSD_BeH2Sto3g) { + if (!check_gpu_available()) + GTEST_SKIP() << "No GPU available, skipping test because CPU is slow"; + + auto pool = cudaq::solvers::operator_pool::get("uccsd"); + auto opt = cudaq::optim::optimizer::get("lbfgs"); + heterogeneous_map config; + config.insert("num-qubits", hambeh2.num_qubits()); + config.insert("num-electrons", 6); + auto poolList = pool->generate(config); + auto [energy, thetas, ops] = cudaq::solvers::adapt_vqe( + statePrep6Electrons, hambeh2, poolList, *opt, "central_difference", + {{"grad_norm_tolerance", 1e-3}, + {"verbose", true}, + {"max_iter", 15}, + {"grad_norm_diff_tolerance", 1e-5}, + {"threshold_energy", 5e-6}, + {"initial_theta", 0.0}, + {"tol", 1e-5}}); + EXPECT_NEAR(energy, -15.59, 1e-2); for (std::size_t i = 0; i < thetas.size(); i++) printf("%lf -> %s\n", thetas[i], ops[i].to_string().c_str()); } \ No newline at end of file diff --git a/libs/solvers/unittests/test_operator_pools.cpp b/libs/solvers/unittests/test_operator_pools.cpp index 8fb9fe60..f9a48cec 100644 --- a/libs/solvers/unittests/test_operator_pools.cpp +++ b/libs/solvers/unittests/test_operator_pools.cpp @@ -7,7 +7,11 @@ ******************************************************************************/ #include "cudaq/solvers/operators/operator_pool.h" +#include +#include #include +#include +#include using namespace cudaqx; @@ -19,7 +23,7 @@ TEST(UCCSDTest, GenerateWithDefaultConfig) { auto operators = pool->generate(config); ASSERT_FALSE(operators.empty()); - EXPECT_EQ(operators.size(), 2 * 2 + 1 * 8); + EXPECT_EQ(operators.size(), 2 + 1); for (const auto &op : operators) { EXPECT_EQ(op.num_qubits(), 4); @@ -30,7 +34,7 @@ TEST(UCCSDTest, GenerateFromAPIFunction) { auto operators = cudaq::solvers::get_operator_pool( "uccsd", {{"num-qubits", 4}, {"num-electrons", 2}}); ASSERT_FALSE(operators.empty()); - EXPECT_EQ(operators.size(), 2 * 2 + 1 * 8); + EXPECT_EQ(operators.size(), 2 + 1); for (const auto &op : operators) { EXPECT_EQ(op.num_qubits(), 4); @@ -46,11 +50,23 @@ TEST(UCCSDTest, GenerateWithCustomCoefficients) { auto operators = pool->generate(config); ASSERT_FALSE(operators.empty()); - EXPECT_EQ(operators.size(), (2 * 2 + 1 * 8)); + EXPECT_EQ(operators.size(), 2 + 1); + std::vector> temp_coeffs; for (size_t i = 0; i < operators.size(); ++i) { EXPECT_EQ(operators[i].num_qubits(), 4); - EXPECT_DOUBLE_EQ(1.0, operators[i].get_coefficient().real()); + + operators[i].for_each_term([&](const auto &term) { + temp_coeffs.push_back(term.get_coefficient()); + }); + } + + for (size_t j = 0; j < temp_coeffs.size(); ++j) { + double real_part = temp_coeffs[j].real(); + EXPECT_TRUE(real_part == 0.5 || real_part == -0.5 || real_part == 0.125 || + real_part == -0.125) + << "Coefficient at index " << j + << " has unexpected value: " << real_part; } } @@ -64,7 +80,7 @@ TEST(UCCSDTest, GenerateWithOddElectrons) { auto operators = pool->generate(config); ASSERT_FALSE(operators.empty()); - EXPECT_EQ(operators.size(), 2 * 4 + 4 * 8); + EXPECT_EQ(operators.size(), 2 * 2 + 4); for (const auto &op : operators) EXPECT_EQ(op.num_qubits(), 6); @@ -79,7 +95,7 @@ TEST(UCCSDTest, GenerateWithLargeSystem) { auto operators = pool->generate(config); ASSERT_FALSE(operators.empty()); - EXPECT_GT(operators.size(), 875); + EXPECT_EQ(operators.size(), 875); for (const auto &op : operators) { EXPECT_EQ(op.num_qubits(), 20); @@ -97,32 +113,59 @@ TEST(UccsdOperatorPoolTest, GeneratesCorrectOperators) { auto operators = pool->generate(config); // Convert SpinOperators to strings - std::vector operator_strings; + std::vector> terms_strings; + std::vector>> coefficients; for (const auto &op : operators) { - operator_strings.push_back(op.to_string(false)); + std::vector> temp_coeffs; + std::vector string_rep; + op.for_each_term([&](const auto &term) { + string_rep.push_back(term.to_string(false)); + temp_coeffs.push_back(term.get_coefficient()); + }); + terms_strings.push_back(string_rep); + coefficients.push_back(temp_coeffs); } // Assert - std::vector expected_operators = { - "YZXI", "XZYI", "IYZX", "IXZY", "XXXY", "XXYX", - "XYYY", "YXYY", "XYXX", "YXXX", "YYXY", "YYYX"}; - - ASSERT_EQ(operator_strings.size(), expected_operators.size()) + std::vector> expected_operators = { + {"XZYI", "YZXI"}, + {"IXZY", "IYZX"}, + {"YYYX", "YXXX", "XXYX", "YYXY", "XYYY", "XXXY", "YXYY", "XYXX"}}; + + std::vector>> expected_coefficients = { + {std::complex(-0.5, 0), std::complex(0.5, 0)}, + {std::complex(-0.5, 0), std::complex(0.5, 0)}, + {std::complex(-0.125, 0), std::complex(-0.125, 0), + std::complex(0.125, 0), std::complex(-0.125, 0), + std::complex(0.125, 0), std::complex(0.125, 0), + std::complex(0.125, 0), std::complex(-0.125, 0)}}; + EXPECT_EQ(terms_strings.size(), expected_operators.size()) << "Number of generated operators does not match expected count"; for (size_t i = 0; i < expected_operators.size(); ++i) { - EXPECT_EQ(operator_strings[i], expected_operators[i]) - << "Mismatch at index " << i; + for (size_t j = 0; j < expected_operators[i].size(); ++j) { + + EXPECT_EQ(expected_operators[i][j].length(), 4) + << "Operator " << expected_operators[i][j] + << " does not have the expected length of 4"; + + EXPECT_EQ(terms_strings[i][j], expected_operators[i][j]) + << "Mismatch at index " << i << ", " << j; + std::cout << coefficients[i][j] << std::endl; + std::cout << expected_coefficients[i][j] << std::endl; + EXPECT_EQ(coefficients[i][j], expected_coefficients[i][j]) + << "Mismatch at index " << i << ", " << j; + } } // Additional checks - for (const auto &op_string : operator_strings) { - EXPECT_EQ(op_string.length(), 4) - << "Operator " << op_string - << " does not have the expected length of 4"; - - EXPECT_TRUE(op_string.find_first_not_of("IXYZ") == std::string::npos) - << "Operator " << op_string << " contains invalid characters"; + for (size_t k = 0; k < terms_strings.size(); ++k) { + for (size_t l = 0; l < terms_strings[k].size(); ++l) { + std::string term_string = terms_strings[k][l]; + EXPECT_EQ(term_string.size(), 4) + << "Operator " << term_string + << " does not have the expected length of 4"; + } } } diff --git a/libs/solvers/unittests/test_optimizers.cpp b/libs/solvers/unittests/test_optimizers.cpp index 600fe859..024a7e19 100644 --- a/libs/solvers/unittests/test_optimizers.cpp +++ b/libs/solvers/unittests/test_optimizers.cpp @@ -54,8 +54,7 @@ TEST(OptimTester, checkLBFGS) { { auto nIters = optimizer->history.size(); // Try to set the tolerance - auto [opt, params] = - optimizer->optimize(2, f, {{"function_tolerance", 1e-3}}); + auto [opt, params] = optimizer->optimize(2, f, {{"tol", 1e-3}}); EXPECT_TRUE(optimizer->history.size() < nIters); }