diff --git a/docs/sphinx/api/solvers/cpp_api.rst b/docs/sphinx/api/solvers/cpp_api.rst index fb40c647..3e47c586 100644 --- a/docs/sphinx/api/solvers/cpp_api.rst +++ b/docs/sphinx/api/solvers/cpp_api.rst @@ -7,6 +7,8 @@ CUDA-Q Solvers C++ API .. doxygenclass:: cudaq::solvers::spin_complement_gsd .. doxygenclass:: cudaq::solvers::uccsd .. doxygenclass:: cudaq::solvers::uccgsd +.. doxygenclass:: cudaq::solvers::upccgsd + .. doxygenclass:: cudaq::solvers::qaoa_pool .. doxygenfunction:: cudaq::solvers::get_operator_pool @@ -70,6 +72,8 @@ CUDA-Q Solvers C++ API .. doxygenfunction:: cudaq::solvers::stateprep::uccsd(cudaq::qview<>, const std::vector&, std::size_t) .. doxygenfunction:: cudaq::solvers::stateprep::get_uccgsd_pauli_lists .. doxygenfunction:: cudaq::solvers::stateprep::uccgsd(cudaq::qview<>, const std::vector&, const std::vector>&, const std::vector>&) +.. doxygenfunction:: cudaq::solvers::stateprep::get_upccgsd_pauli_lists +.. doxygenfunction:: cudaq::solvers::stateprep::upccgsd(cudaq::qview<>, const std::vector&, const std::vector>&, const std::vector>&) .. doxygenstruct:: cudaq::solvers::qaoa_result diff --git a/docs/sphinx/api/solvers/python_api.rst b/docs/sphinx/api/solvers/python_api.rst index 8e82a1fd..b1c37374 100644 --- a/docs/sphinx/api/solvers/python_api.rst +++ b/docs/sphinx/api/solvers/python_api.rst @@ -28,7 +28,11 @@ CUDA-Q Solvers Python API .. autofunction:: cudaq_solvers.stateprep.get_num_uccsd_parameters .. autofunction:: cudaq_solvers.stateprep.get_uccsd_excitations .. autofunction:: cudaq_solvers.stateprep.get_uccgsd_pauli_lists +.. autofunction:: cudaq_solvers.stateprep.get_upccgsd_pauli_lists + .. autofunction:: cudaq_solvers.stateprep.uccgsd +.. autofunction:: cudaq_solvers.stateprep.upccgsd + .. autofunction:: cudaq_solvers.get_num_qaoa_parameters diff --git a/docs/sphinx/components/solvers/introduction.rst b/docs/sphinx/components/solvers/introduction.rst index 26b8287c..2299973c 100644 --- a/docs/sphinx/components/solvers/introduction.rst +++ b/docs/sphinx/components/solvers/introduction.rst @@ -509,6 +509,12 @@ CUDA-QX provides several pre-built operator pools for ADAPT-VQE: single and double excitations, regardless of their occupied/virtual status in the reference state. Capable of capturing both dynamic and static (strong) correlation but at the cost of increased circuit depth and parameter count. +* **upccgsd**: UCC generalized singles and paired doubles. + A structured, lower-depth variant of UCCGSD in which double excitations are restricted to paired electron + transfers between spatial orbitals, following the UpCC (pair-cluster) construction. + While less expressive than full UCCGSD, it retains generalized spin-preserving singles + and the most physically relevant pair-correlation channels, substantially reducing circuit depth and parameter count. + This makes UpCCGSD attractive for larger systems or hardware-constrained regimes where UCCGSD is too costly. * **qaoa**: QAOA mixer excitation operators It generates all possible single-qubit X and Y terms, along with all possible two-qubit interaction terms (XX, YY, XY, YX, XZ, ZX, YZ, ZY) across every pair of qubits. @@ -537,6 +543,10 @@ CUDA-QX provides several pre-built operator pools for ADAPT-VQE: "uccgsd", num_orbitals=molecule.n_orbitals ) + upccgsd_ops = solvers.get_operator_pool( + "upccgsd", + num_orbitals=molecule.n_orbitals +) Available Ansatz ^^^^^^^^^^^^^^^^^^ @@ -545,6 +555,7 @@ CUDA-QX provides several state preparations ansatz for VQE. * **uccsd**: UCCSD operators * **uccgsd**: UCC generalized singles and doubles +* **upccgsd**: UCC generalized singles and paired doubles .. code-block:: python @@ -585,6 +596,26 @@ CUDA-QX provides several state preparations ansatz for VQE. for i in range(numElectrons): x(q[i]) solvers.stateprep.uccgsd(q, thetas, pauliWordsList, coefficientsList) + + # Using UpCCGSD ansatz + geometry = [('H', (0., 0., 0.)), ('H', (0., 0., .7474))] + molecule = solvers.create_molecule(geometry, 'sto-3g', 0, 0, casci=True) + + numQubits = molecule.n_orbitals * 2 + numElectrons = molecule.n_electrons + + # Get grouped Pauli words and coefficients from UpCCGSD pool + pauliWordsList, coefficientsList = solvers.stateprep.get_upccgsd_pauli_lists( + numQubits, only_singles=False, only_doubles=False) + + @cudaq.kernel + def ansatz(numQubits: int, numElectrons: int, thetas: list[float], + pauliWordsList: list[list[cudaq.pauli_word]], + coefficientsList: list[list[float]]): + q = cudaq.qvector(numQubits) + for i in range(numElectrons): + x(q[i]) + solvers.stateprep.upccgsd(q, thetas, pauliWordsList, coefficientsList) Algorithm Parameters diff --git a/libs/solvers/include/cudaq/solvers/operators/operator_pools/upccgsd_operator_pool.h b/libs/solvers/include/cudaq/solvers/operators/operator_pools/upccgsd_operator_pool.h new file mode 100644 index 00000000..fd8f3c83 --- /dev/null +++ b/libs/solvers/include/cudaq/solvers/operators/operator_pools/upccgsd_operator_pool.h @@ -0,0 +1,47 @@ +/****************************************************************-*- C++ -*-**** + * Copyright (c) 2025 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +#pragma once + +#include "../operator_pool.h" + +namespace cudaq::solvers { + +/// @brief UpCCGSD operator pool is a class which generates a pool +/// of UpCCGSD 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. The UpCCGSD pool restricts generalized +/// doubles to paired αβ→αβ excitations while retaining spin-preserving +/// generalized singles. +class upccgsd : public operator_pool { + +public: + /// @brief Generate a vector of spin operators based on the provided + /// configuration. + /// @details The UpCCGSD 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, including the number of spatial orbitals + /// ("num-orbitals" or "num_orbitals"). + /// @return A vector of cudaq::spin_op objects representing the generated + /// UpCCGSD 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 upccgsd class. + CUDAQ_EXTENSION_CREATOR_FUNCTION(operator_pool, upccgsd) +}; + +/// @brief Register the upccgsd extension type with the CUDA-Q framework +CUDAQ_REGISTER_TYPE(upccgsd) + +} // namespace cudaq::solvers diff --git a/libs/solvers/include/cudaq/solvers/stateprep/upccgsd.h b/libs/solvers/include/cudaq/solvers/stateprep/upccgsd.h new file mode 100644 index 00000000..422202cb --- /dev/null +++ b/libs/solvers/include/cudaq/solvers/stateprep/upccgsd.h @@ -0,0 +1,42 @@ +/****************************************************************-*- C++ -*-**** + * Copyright (c) 2025 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ +#pragma once + +#include "cudaq.h" +#include +#include + +namespace cudaq::solvers::stateprep { + +/// @brief Generate UpCCGSD operator pool (Python-style unique singles/doubles) +/// and extract Pauli words and coefficients. +/// @details The UpCCGSD pool is constructed using spin-preserving generalized +/// singles and paired αβ→αβ generalized doubles in an interleaved spin-orbital +/// ordering (α0, β0, α1, β1, ..., α_{N-1}, β_{N-1}). +/// @param norbitals Number of spin orbitals (qubits) +/// @param only_doubles If true, include only double excitations +/// @return Pair of lists: [Pauli words grouped by excitation], +/// [coefficients grouped by excitation] +std::pair>, + std::vector>> +get_upccgsd_pauli_lists(std::size_t norbitals,bool only_doubles = false); + +/// \pure_device_kernel +/// +/// @brief Apply UpCCGSD ansatz to a qubit register using grouped Pauli words +/// and coefficients. +/// @param qubits Qubit register +/// @param thetas Vector of rotation angles (one per excitation group) +/// @param pauliWordsList Pauli words grouped by excitation +/// @param coefficientsList Coefficients grouped by excitation +__qpu__ void +upccgsd(cudaq::qview<> qubits, const std::vector &thetas, + const std::vector> &pauliWordsList, + const std::vector> &coefficientsList); + +} // namespace cudaq::solvers::stateprep diff --git a/libs/solvers/lib/CMakeLists.txt b/libs/solvers/lib/CMakeLists.txt index d8137cbf..da23edc6 100644 --- a/libs/solvers/lib/CMakeLists.txt +++ b/libs/solvers/lib/CMakeLists.txt @@ -24,6 +24,7 @@ add_library(cudaq-solvers SHARED operators/operator_pools/uccsd_operator_pool.cpp operators/operator_pools/qaoa_operator_pool.cpp operators/operator_pools/uccgsd_operator_pool.cpp + operators/operator_pools/upccgsd_operator_pool.cpp operators/uccgsd_excitation_utils.cpp version.cpp ) diff --git a/libs/solvers/lib/operators/operator_pools/upccgsd_operator_pool.cpp b/libs/solvers/lib/operators/operator_pools/upccgsd_operator_pool.cpp new file mode 100644 index 00000000..1db72a1a --- /dev/null +++ b/libs/solvers/lib/operators/operator_pools/upccgsd_operator_pool.cpp @@ -0,0 +1,116 @@ +/******************************************************************************* + * Copyright (c) 2025 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * 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/upccgsd_operator_pool.h" +#include "cudaq/solvers/operators/uccgsd_excitation_utils.h" + +using namespace cudaqx; + +namespace cudaq::solvers { + +std::vector +upccgsd::generate(const heterogeneous_map &config) const { + + // Number of spatial orbitals (same convention as UCCGSD operator pool) + const std::size_t numOrbitals = + config.get({"num-orbitals", "num_orbitals"}); + if (numOrbitals == 0) + throw std::invalid_argument("num-orbitals must be > 0"); + + // Interleaved spin-orbitals: (α0, β0, α1, β1, ..., α_{N-1}, β_{N-1}) + const std::size_t numQubits = 2 * numOrbitals; + + std::vector ops; + + + // 1) Spin-preserving generalized singles on α and β + // + // We use generate_uccgsd_singles(numQubits) to get all unordered pairs + // (p,q) with p > q over spin-orbitals, then restrict to same-spin pairs: + // - interleaved mapping means: + // even indices -> α spin + // odd indices -> β spin + // - spin-preserving => p % 2 == q % 2 + // + // Each selected pair is mapped to a JW single-excitation operator via + // addUCCGSDSingleExcitation. + auto singles = generate_uccgsd_singles(numQubits); + for (const auto &[p, q] : singles) { + // keep only spin-preserving single excitations (α→α or β→β) + if ((p % 2) == (q % 2)) { + addUCCGSDSingleExcitation(ops, p, q); + } + } + + // 2) Paired generalized doubles (UpCC-style) in interleaved ordering + // + // UpCCGSD pair doubles correspond to exciting an αβ pair on spatial orbital p + // to an αβ pair on spatial orbital q: + // + // | pα, pβ > -> | qα, qβ > + // + // in a GENERALIZED sense ("G" in UpCCGSD): we consider all unordered pairs + // of distinct spatial orbitals p < q. + // + // Interleaved spin-orbital indices: + // spatial orbital o: + // o_α -> 2*o + // o_β -> 2*o + 1 + // + // For a given spatial pair (p, q), p < q: + // p_alpha = 2*p + // p_beta = 2*p + 1 + // q_alpha = 2*q + // q_beta = 2*q + 1 + // + // We then form two spin-orbital *pairs* to pass to addUCCGSDDoubleExcitation: + // + // creation pair (P, Q) = sorted(q_alpha, q_beta) with P > Q + // annihilation pair (R, S) = sorted(p_alpha, p_beta) with R > S + // + // This satisfies the helper's requirement p > q and r > s, and reuses + // the same JW mapping as the UCCGSD double-excitation utility, but now + // restricted to UpCC-style αβ pair doubles. + + for (std::size_t p = 0; p < numOrbitals; ++p) { + for (std::size_t q = p + 1; q < numOrbitals; ++q) { + // Interleaved spin-orbital indices for p and q + const std::size_t p_alpha = 2 * p; + const std::size_t p_beta = 2 * p + 1; + const std::size_t q_alpha = 2 * q; + const std::size_t q_beta = 2 * q + 1; + + // Sort each pair so that first > second, as required by the helper. + std::size_t P, Qp; + if (q_beta > q_alpha) { + P = q_beta; + Qp = q_alpha; + } else { + P = q_alpha; + Qp = q_beta; + } + + std::size_t R, S; + if (p_beta > p_alpha) { + R = p_beta; + S = p_alpha; + } else { + R = p_alpha; + S = p_beta; + } + + // Use the same double-excitation JW mapping helper as UCCGSD, + // but restricted to these UpCC-style pair indices. + addUCCGSDDoubleExcitation(ops, P, Qp, R, S); + } + } + + return ops; +} + +} // namespace cudaq::solvers diff --git a/libs/solvers/lib/stateprep/CMakeLists.txt b/libs/solvers/lib/stateprep/CMakeLists.txt index 4b8ba5fc..e4777473 100644 --- a/libs/solvers/lib/stateprep/CMakeLists.txt +++ b/libs/solvers/lib/stateprep/CMakeLists.txt @@ -10,4 +10,5 @@ cudaqx_add_device_code(cudaq-solvers SOURCES uccsd.cpp uccgsd.cpp + upccgsd.cpp ) diff --git a/libs/solvers/lib/stateprep/upccgsd.cpp b/libs/solvers/lib/stateprep/upccgsd.cpp new file mode 100644 index 00000000..25b68b3b --- /dev/null +++ b/libs/solvers/lib/stateprep/upccgsd.cpp @@ -0,0 +1,159 @@ +/******************************************************************************* + * Copyright (c) 2025 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * 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/stateprep/upccgsd.h" +#include "cudaq/solvers/operators/uccgsd_excitation_utils.h" +#include +#include +#include + +namespace cudaq::solvers::stateprep { + +// Wrapper helpers that reuse the shared JW excitation utilities + + +static void addUpCCGSDSingleExcitation(std::vector &ops, + std::size_t p, std::size_t q) { + // Reuse the UCCGSD single-excitation builder + cudaq::solvers::addUCCGSDSingleExcitation(ops, p, q); +} + +static void addUpCCGSDDoubleExcitation(std::vector &ops, + std::size_t p, std::size_t q, + std::size_t r, std::size_t s) { + // Reuse the UCCGSD double-excitation builder + cudaq::solvers::addUCCGSDDoubleExcitation(ops, p, q, r, s); +} + +// UpCCGSD unique singles: spin-preserving generalized singles +// norbitals is interpreted as the number of spin orbitals. +std::vector> +upccgsd_unique_singles(std::size_t norbitals) { + // Start from all generalized singles over spin orbitals + auto allSingles = cudaq::solvers::generate_uccgsd_singles(norbitals); + + std::vector> filtered; + filtered.reserve(allSingles.size()); + + for (auto [p, q] : allSingles) { + // Interleaved mapping: even -> α, odd -> β + // Keep only spin-preserving singles (α→α or β→β) + if ((p % 2) == (q % 2)) + filtered.emplace_back(p, q); + } + + return filtered; +} + +// UpCCGSD unique doubles: paired αβ@p -> αβ@q generalized doubles +// norbitals is the number of spin orbitals (qubits) arranged interleaved as +// (α0, β0, α1, β1, ..., α_{N-1}, β_{N-1}), so N_spatial = norbitals / 2. +std::vector, + std::pair>> +upccgsd_unique_doubles(std::size_t norbitals) { + if (norbitals % 2 != 0) + throw std::invalid_argument( + "upccgsd_unique_doubles expects an even number of spin orbitals."); + + const std::size_t numSpatialOrbitals = norbitals / 2; + + std::vector, + std::pair>> + doubles; + + for (std::size_t p = 0; p < numSpatialOrbitals; ++p) { + for (std::size_t q = p + 1; q < numSpatialOrbitals; ++q) { + // Interleaved spin-orbital indices for spatial orbitals p and q + const std::size_t p_alpha = 2 * p; + const std::size_t p_beta = 2 * p + 1; + const std::size_t q_alpha = 2 * q; + const std::size_t q_beta = 2 * q + 1; + + // Sort each pair so that first > second, as required by the helper. + std::size_t P, Qp; + if (q_beta > q_alpha) { + P = q_beta; + Qp = q_alpha; + } else { + P = q_alpha; + Qp = q_beta; + } + + std::size_t R, S; + if (p_beta > p_alpha) { + R = p_beta; + S = p_alpha; + } else { + R = p_alpha; + S = p_beta; + } + + doubles.push_back({{P, Qp}, {R, S}}); + } + } + + return doubles; +} + +// Python-style UpCCGSD pool: build Pauli-word lists + coefficients +// norbitals = # spin orbitals / qubits (same as uccgsd stateprep). +std::pair>, + std::vector>> +get_upccgsd_pauli_lists(std::size_t norbitals, bool only_doubles) { + + std::vector ops; + + if (only_doubles) { + // Only UpCCGSD paired doubles + for (auto pair : upccgsd_unique_doubles(norbitals)) { + auto [pq, rs] = pair; + addUpCCGSDDoubleExcitation(ops, pq.first, pq.second, rs.first, rs.second); + } + } else { + // Default: add all UpCCGSD singles... + for (auto [p, q] : upccgsd_unique_singles(norbitals)) + addUpCCGSDSingleExcitation(ops, p, q); + // ...and all UpCCGSD paired doubles + for (auto pair : upccgsd_unique_doubles(norbitals)) { + auto [pq, rs] = pair; + addUpCCGSDDoubleExcitation(ops, pq.first, pq.second, rs.first, rs.second); + } + } + + std::vector> pauliWordsList; + std::vector> coefficientsList; + + for (const auto &op : ops) { + std::vector words; + std::vector coeffs; + for (const auto &term : op) { + words.push_back(term.get_pauli_word(norbitals)); + coeffs.push_back(term.evaluate_coefficient().real()); + } + pauliWordsList.push_back(std::move(words)); + coefficientsList.push_back(std::move(coeffs)); + } + + return {pauliWordsList, coefficientsList}; +} +// QPU kernel: exp(i θ_k H_k) with H_k expanded into Pauli words +__qpu__ void +upccgsd(cudaq::qview<> qubits, const std::vector &thetas, + const std::vector> &pauliWordsList, + const std::vector> &coefficientsList) { + for (std::size_t i = 0; i < pauliWordsList.size(); ++i) { + // Use the same theta for all terms in this excitation + double theta = thetas[i]; + const auto &words = pauliWordsList[i]; + const auto &coeffs = coefficientsList[i]; + for (std::size_t j = 0; j < words.size(); ++j) { + exp_pauli(theta * coeffs[j], qubits, words[j]); + } + } +} + +} // namespace cudaq::solvers::stateprep diff --git a/libs/solvers/python/bindings/solvers/py_solvers.cpp b/libs/solvers/python/bindings/solvers/py_solvers.cpp index acb39da2..dc44c72c 100644 --- a/libs/solvers/python/bindings/solvers/py_solvers.cpp +++ b/libs/solvers/python/bindings/solvers/py_solvers.cpp @@ -16,6 +16,8 @@ #include "cudaq/solvers/adapt.h" #include "cudaq/solvers/qaoa.h" #include "cudaq/solvers/stateprep/uccgsd.h" +#include "cudaq/solvers/stateprep/upccgsd.h" + #include "cudaq/solvers/stateprep/uccsd.h" #include "cudaq/solvers/version.h" #include "cudaq/solvers/vqe.h" @@ -172,6 +174,15 @@ void addStatePrepKernels(py::module &mod) { "Takes as input the qubits, grouped rotational parameters, grouped Pauli " "words, " "and grouped coefficients."); + cudaq::python::addDeviceKernelInterop< + cudaq::qview<>, const std::vector &, + const std::vector> &, + const std::vector> &>( + mod, "stateprep", "upccgsd", + "Unitary Coupled Cluster Generalized Singles and Paired Doubles Ansatz. " + "Takes as input the qubits, grouped rotational parameters, grouped Pauli " + "words, " + "and grouped coefficients."); cudaq::python::addDeviceKernelInterop, double, std::size_t, std::size_t>( mod, "stateprep", "single_excitation", @@ -202,6 +213,21 @@ void addStatePrepKernels(py::module &mod) { R"( Generate UCCGSD operator pool (Python-style unique singles/doubles) and extract Pauli words and coefficients grouped by excitation. + Args: + num_qubits (int): Number of spin orbitals (qubits). + only_singles (bool): If True, only single excitations. + only_doubles (bool): If True, only double excitations. + + Returns: + Tuple[List[List[PauliWord]], List[List[float]]]: Pauli words and coefficients grouped by excitation. + )"); + stateprep.def("get_upccgsd_pauli_lists", + &cudaq::solvers::stateprep::get_upccgsd_pauli_lists, + py::arg("num_qubits"), + py::arg("only_doubles") = false, + R"( + Generate UpCCGSD operator pool (Python-style unique singles/doubles) and extract Pauli words and coefficients grouped by excitation. + Args: num_qubits (int): Number of spin orbitals (qubits). only_singles (bool): If True, only single excitations. diff --git a/libs/solvers/python/tests/test_upccgsd.py b/libs/solvers/python/tests/test_upccgsd.py new file mode 100644 index 00000000..fc1d5edf --- /dev/null +++ b/libs/solvers/python/tests/test_upccgsd.py @@ -0,0 +1,96 @@ +import os + +import pytest +import numpy as np + +import cudaq +import cudaq_solvers as solvers +from scipy.optimize import minimize +import subprocess + +def is_nvidia_gpu_available(): + """Check if NVIDIA GPU is available using nvidia-smi command.""" + try: + result = subprocess.run(["nvidia-smi"], + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + text=True) + if result.returncode == 0 and "GPU" in result.stdout: + return True + except FileNotFoundError: + # The nvidia-smi command is not found, indicating no NVIDIA GPU drivers + return False + return False + +def test_solvers_upccgsd_exc_list(): + N = 20 # or whatever molecule.n_orbitals * 2 would be + pauliWordsList, coefficientsList = solvers.stateprep.get_upccgsd_pauli_lists( + N, only_doubles=False + ) + parameter_count = len(coefficientsList) + M = N/2 + ideal_count = (3/2) * M * (M-1) + assert parameter_count == ideal_count + pauliWordsList, coefficientsList = solvers.stateprep.get_upccgsd_pauli_lists( + N, only_doubles=True + ) + parameter_count = len(coefficientsList) + M = N/2 + ideal_count = (1/2) * M * (M-1) + assert parameter_count == ideal_count + + + + + +def test_solvers_vqe_upccgsd_h2(): + + geometry = [('H', (0., 0., 0.)), ('H', (0., 0., .7474))] + molecule = solvers.create_molecule(geometry, 'sto-3g', 0, 0, casci=True) + + numQubits = molecule.n_orbitals * 2 + numElectrons = molecule.n_electrons + + # Get grouped Pauli words and coefficients from UpCCGSD pool + pauliWordsList, coefficientsList = solvers.stateprep.get_upccgsd_pauli_lists( + numQubits, only_doubles=False) + + # Number of theta parameters = number of excitation groups + parameter_count = len(coefficientsList) + assert parameter_count == 3 + + @cudaq.kernel + def ansatz(numQubits: int, numElectrons: int, thetas: list[float], + pauliWordsList: list[list[cudaq.pauli_word]], + coefficientsList: list[list[float]]): + q = cudaq.qvector(numQubits) + for i in range(numElectrons): + x(q[i]) + # Apply UpCCGSD circuit with grouped thetas + solvers.stateprep.upccgsd(q, thetas, pauliWordsList, coefficientsList) + + x0 = [0.0 for _ in range(parameter_count)] + + def cost(theta): + + theta = theta.tolist() + + energy = cudaq.observe(ansatz, molecule.hamiltonian, numQubits, + numElectrons, theta, pauliWordsList, + coefficientsList).expectation() + return energy + + res = minimize(cost, + x0, + method='COBYLA', + options={ + 'maxiter': 1000, + 'rhobeg': 0.1, + 'disp': True + }) + energy = res.fun + + hf_energy = molecule.energies.get("HF", None) + + if hf_energy is not None: + assert energy <= hf_energy + 1e-4 \ No newline at end of file diff --git a/libs/solvers/unittests/CMakeLists.txt b/libs/solvers/unittests/CMakeLists.txt index a6d4c365..c4a7a642 100644 --- a/libs/solvers/unittests/CMakeLists.txt +++ b/libs/solvers/unittests/CMakeLists.txt @@ -74,6 +74,11 @@ target_link_libraries(test_uccgsd_operator_pool PRIVATE GTest::gtest_main cudaq- add_dependencies(CUDAQXSolversUnitTests test_uccgsd_operator_pool) gtest_discover_tests(test_uccgsd_operator_pool) +add_executable(test_upccgsd_operator_pool test_upccgsd_operator_pool.cpp) +target_link_libraries(test_upccgsd_operator_pool PRIVATE GTest::gtest_main cudaq-solvers cudaq::cudaq) +add_dependencies(CUDAQXSolversUnitTests test_upccgsd_operator_pool) +gtest_discover_tests(test_upccgsd_operator_pool) + add_executable(test_qaoa test_qaoa.cpp) target_link_libraries(test_qaoa PRIVATE GTest::gtest_main cudaq-solvers test-kernels cudaq::cudaq) gtest_discover_tests(test_qaoa) diff --git a/libs/solvers/unittests/test_upccgsd_operator_pool.cpp b/libs/solvers/unittests/test_upccgsd_operator_pool.cpp new file mode 100644 index 00000000..fc237034 --- /dev/null +++ b/libs/solvers/unittests/test_upccgsd_operator_pool.cpp @@ -0,0 +1,367 @@ +/******************************************************************************* + * Copyright (c) 2025 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * 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_pool.h" +#include "cudaq/solvers/stateprep/upccgsd.h" +#include +#include +#include +#include + +using namespace cudaqx; + +// Helper function to count operators by degree (singles vs doubles) +std::pair +countSinglesAndDoubles(const std::vector &ops) { + size_t singles = 0, doubles = 0; + for (const auto &op : ops) { + // Heuristic: count X/Y in the string representation. + std::string op_str = op.to_string(); + + size_t xy_count = 0; + for (size_t i = 0; i + 1 < op_str.length(); ++i) { + if ((op_str[i] == 'X' || op_str[i] == 'Y') && + std::isdigit(op_str[i + 1])) { + xy_count++; + } + } + + // Same threshold logic as before: singles are "small", doubles "large". + if (xy_count <= 10) + singles++; + else + doubles++; + } + return {singles, doubles}; +} + +// ============================================================================ +// Test 1: Correct Number of Operators +// ============================================================================ +TEST(UPCCGSDOperatorPoolTest, CorrectNumberOfOperators) { + auto pool = cudaq::solvers::operator_pool::get("upccgsd"); + + // Recall: M = numOrbitals, n_qubits = 2M + // singles = M (M - 1) + // doubles = M (M - 1) / 2 + // total = 3/2 M (M - 1) + + // Test case 1: 2 orbitals (4 qubits) + { + heterogeneous_map config; + config.insert("num-orbitals", 2); + auto ops = pool->generate(config); + + std::size_t M = 2; + std::size_t expected_singles = M * (M - 1); // 2 + std::size_t expected_doubles = M * (M - 1) / 2; // 1 + std::size_t expected_total = expected_singles + expected_doubles; // 3 + + EXPECT_EQ(ops.size(), expected_total) + << "For 2 orbitals: expected " << expected_total << " operators"; + } + + // Test case 2: 3 orbitals (6 qubits) + { + heterogeneous_map config; + config.insert("num-orbitals", 3); + auto ops = pool->generate(config); + + std::size_t M = 3; + std::size_t expected_singles = M * (M - 1); // 6 + std::size_t expected_doubles = M * (M - 1) / 2; // 3 + std::size_t expected_total = expected_singles + expected_doubles; // 9 + + EXPECT_EQ(ops.size(), expected_total) + << "For 3 orbitals: expected " << expected_total << " operators"; + } + + // Test case 3: 4 orbitals (8 qubits) + { + heterogeneous_map config; + config.insert("num-orbitals", 4); + auto ops = pool->generate(config); + + std::size_t M = 4; + std::size_t expected_singles = M * (M - 1); // 12 + std::size_t expected_doubles = M * (M - 1) / 2; // 6 + std::size_t expected_total = expected_singles + expected_doubles; // 18 + + EXPECT_EQ(ops.size(), expected_total) + << "For 4 orbitals: expected " << expected_total << " operators"; + } +} + +// ============================================================================ +// Test 2: No Duplicate Operators +// ============================================================================ +TEST(UPCCGSDOperatorPoolTest, NoDuplicateOperators) { + auto pool = cudaq::solvers::operator_pool::get("upccgsd"); + + heterogeneous_map config; + config.insert("num-orbitals", 3); // 6 qubits + auto ops = pool->generate(config); + + std::set qubit_set; + for (size_t i = 0; i < 6; ++i) + qubit_set.insert(i); + + std::set unique_ops; + for (auto op : ops) { + op.canonicalize(qubit_set); + unique_ops.insert(op.to_string()); + } + + EXPECT_EQ(unique_ops.size(), ops.size()) + << "Found duplicate operators! Expected " << ops.size() + << " unique operators but got " << unique_ops.size(); +} + +// ============================================================================ +// Test 3: Verify Singles and Doubles Separately +// ============================================================================ +TEST(UPCCGSDOperatorPoolTest, CorrectSinglesAndDoublesCount) { + auto pool = cudaq::solvers::operator_pool::get("upccgsd"); + + heterogeneous_map config; + config.insert("num-orbitals", 3); // 6 qubits + auto ops = pool->generate(config); + + auto [singles_count, doubles_count] = countSinglesAndDoubles(ops); + + std::size_t M = 3; + std::size_t expected_singles = M * (M - 1); // 6 + std::size_t expected_doubles = M * (M - 1) / 2; // 3 + + EXPECT_EQ(singles_count, expected_singles) + << "Expected " << expected_singles << " single excitations"; + EXPECT_EQ(doubles_count, expected_doubles) + << "Expected " << expected_doubles << " double excitations"; +} + +// ============================================================================ +// Test 4: Verify All Operators Are Non-Empty +// ============================================================================ +TEST(UPCCGSDOperatorPoolTest, AllOperatorsNonEmpty) { + auto pool = cudaq::solvers::operator_pool::get("upccgsd"); + heterogeneous_map config; + config.insert("num-orbitals", 2); + auto ops = pool->generate(config); + + for (size_t i = 0; i < ops.size(); ++i) { + EXPECT_GT(ops[i].num_terms(), 0) + << "Operator " << i << " is empty (has 0 terms)"; + for (const auto &term : ops[i]) + EXPECT_FALSE(term.is_identity()) + << "Term " << term.to_string() << " is identity"; + } +} + +// ============================================================================ +// Test 5: Verify Operator Coefficients +// (same coefficients as UCCGSD, just fewer operators) +// ============================================================================ +TEST(UPCCGSDOperatorPoolTest, VerifyOperatorCoefficients) { + auto pool = cudaq::solvers::operator_pool::get("upccgsd"); + heterogeneous_map config; + config.insert("num-orbitals", 2); + auto ops = pool->generate(config); + + std::vector coefficients; + for (const auto &op : ops) { + for (const auto &term : op) { + double coeff = std::abs(term.evaluate_coefficient().real()); + coefficients.push_back(coeff); + } + } + + for (size_t i = 0; i < coefficients.size(); ++i) { + bool is_valid = (std::abs(coefficients[i] - 0.5) < 1e-10) || + (std::abs(coefficients[i] - 0.125) < 1e-10); + EXPECT_TRUE(is_valid) << "Coefficient " << i + << " has unexpected value: " << coefficients[i] + << " (expected 0.5 or 0.125)"; + } +} + +// ============================================================================ +// Test 6: Consistency with Stateprep Version +// ============================================================================ +TEST(UPCCGSDOperatorPoolTest, ConsistentWithStateprep) { + auto pool = cudaq::solvers::operator_pool::get("upccgsd"); + heterogeneous_map config; + config.insert("num-orbitals", 2); // 4 qubits + auto pool_ops = pool->generate(config); + + // stateprep uses "norbitals" = #spin orbitals = #qubits + std::size_t n_qubits = 4; + auto [pauli_lists, coeff_lists] = + cudaq::solvers::stateprep::get_upccgsd_pauli_lists(n_qubits, false); + + EXPECT_EQ(pool_ops.size(), pauli_lists.size()) + << "Operator pool and stateprep should generate same number of " + "operators. " + << "Pool: " << pool_ops.size() << ", Stateprep: " << pauli_lists.size(); + + for (size_t i = 0; i < std::min(pool_ops.size(), pauli_lists.size()); ++i) { + EXPECT_EQ(pool_ops[i].num_terms(), pauli_lists[i].size()) + << "Operator " << i << " has different number of terms. " + << "Pool: " << pool_ops[i].num_terms() + << ", Stateprep: " << pauli_lists[i].size(); + } +} + + +// ============================================================================ +// Test 7: Doubles Only Option (via stateprep) +// ============================================================================ +TEST(UPCCGSDOperatorPoolTest, DoublesOnlyGeneration) { + std::size_t M = 3; // spatial orbitals + std::size_t n_qubits = 2*M; // spin orbitals + + auto [pauli_lists, coeff_lists] = + cudaq::solvers::stateprep::get_upccgsd_pauli_lists( + n_qubits, /*only_doubles=*/true); + + std::size_t expected_doubles = M * (M - 1) / 2; // paired αβ→αβ + EXPECT_EQ(pauli_lists.size(), expected_doubles) + << "Doubles-only should generate " << expected_doubles << " operators"; +} + +// ============================================================================ +// Test 8: Scaling Test - Verify Formula for Multiple Sizes +// ============================================================================ +TEST(UPCCGSDOperatorPoolTest, ScalingBehavior) { + auto pool = cudaq::solvers::operator_pool::get("upccgsd"); + + // (M = numOrbitals, expected total count) + std::vector> test_cases = { + {2, 3}, // M=2: singles=2, doubles=1 + {3, 9}, // M=3: singles=6, doubles=3 + {4, 18}, // M=4: singles=12, doubles=6 + {5, 30}, // M=5: singles=20, doubles=10 + }; + + for (auto [M, expected_count] : test_cases) { + heterogeneous_map config; + config.insert("num-orbitals", M); + auto ops = pool->generate(config); + + std::size_t singles = M * (M - 1); + std::size_t doubles = M * (M - 1) / 2; + std::size_t total = singles + doubles; + + EXPECT_EQ(ops.size(), total) + << "For M=" << M << " orbitals: expected " + << total << " operators, got " << ops.size(); + EXPECT_EQ(total, expected_count) + << "Formula mismatch for M=" << M << " orbitals"; + } +} + +// ============================================================================ +// Test 9: Edge Case - Minimal System +// ============================================================================ +TEST(UPCCGSDOperatorPoolTest, MinimalSystem) { + auto pool = cudaq::solvers::operator_pool::get("upccgsd"); + + // With 1 orbital (2 qubits), there are no spin-preserving generalized + // singles and no paired doubles, so the pool should be empty. + { + heterogeneous_map config; + config.insert("num-orbitals", 1); + auto ops = pool->generate(config); + EXPECT_EQ(ops.size(), 0) + << "With 1 orbital, UpCCGSD should generate 0 operators"; + } + + // With 2 orbitals, we should have 2 singles + 1 double = 3 operators. + { + heterogeneous_map config; + config.insert("num-orbitals", 2); + auto ops = pool->generate(config); + EXPECT_EQ(ops.size(), 3) + << "With 2 orbitals, UpCCGSD should generate 3 operators"; + } +} + +// ============================================================================ +// Test 10: Verify Hermiticity / Anti-Hermiticity Structure +// ============================================================================ +TEST(UPCCGSDOperatorPoolTest, OperatorsAreHermitianGenerators) { + auto pool = cudaq::solvers::operator_pool::get("upccgsd"); + heterogeneous_map config; + config.insert("num-orbitals", 2); + auto ops = pool->generate(config); + + for (size_t i = 0; i < ops.size(); ++i) { + auto matrix_G = ops[i].to_matrix(); + + // Check G is Hermitian + auto adjoint_G = ops[i].to_matrix(); + for (size_t row = 0; row < matrix_G.rows(); ++row) + for (size_t col = 0; col < matrix_G.cols(); ++col) + adjoint_G(row, col) = std::conj(matrix_G(col, row)); + + bool is_hermitian = true; + for (size_t row = 0; row < matrix_G.rows(); ++row) { + for (size_t col = 0; col < matrix_G.cols(); ++col) { + auto diff = adjoint_G(row, col) - matrix_G(row, col); + if (std::abs(diff) > 1e-10) { + is_hermitian = false; + break; + } + } + if (!is_hermitian) + break; + } + + EXPECT_TRUE(is_hermitian) << "Operator " << i << " (G) is not Hermitian"; + } +} + +// ============================================================================ +// Test 11: Verify Operator Pool Returns Correct Type +// ============================================================================ +TEST(UPCCGSDOperatorPoolTest, ReturnsSpinOperators) { + auto pool = cudaq::solvers::operator_pool::get("upccgsd"); + ASSERT_NE(pool, nullptr) << "Failed to get upccgsd operator pool"; + + heterogeneous_map config; + config.insert("num-orbitals", 2); + auto ops = pool->generate(config); + + for (size_t i = 0; i < ops.size(); ++i) { + EXPECT_NO_THROW(ops[i].to_string()) + << "Operator " << i << " is not a valid spin_op"; + } +} + +// ============================================================================ +// Test 12: Performance Test - Larger System +// ============================================================================ +TEST(UPCCGSDOperatorPoolTest, LargeSystemPerformance) { + auto pool = cudaq::solvers::operator_pool::get("upccgsd"); + heterogeneous_map config; + config.insert("num-orbitals", 8); // 16 qubits + + auto start = std::chrono::high_resolution_clock::now(); + auto ops = pool->generate(config); + auto end = std::chrono::high_resolution_clock::now(); + + auto duration = + std::chrono::duration_cast(end - start); + + EXPECT_LT(duration.count(), 5000) + << "Generation took too long: " << duration.count() << "ms"; + + std::size_t M = 8; + std::size_t expected = + M * (M - 1) + M * (M - 1) / 2; // singles + doubles + EXPECT_EQ(ops.size(), expected); +}