Closes #557, in-place implementation of expval for LightningQubit

.github/CHANGELOG.md

@@ -15,6 +15,9 @@
 
 ### Improvements
 
+* Modify expval of named operators in LightningQubit for in-place computation of expectation value, to avoid creating an intermediate statevector
+  [(##565)] (https://github.com/PennyLaneAI/pennylane-lightning/pull/565)
+
 * Modify `setup.py` to use backend-specific build directory (`f"build_{backend}"`) to accelerate rebuilding backends in alternance.
   [(#540)] (https://github.com/PennyLaneAI/pennylane-lightning/pull/540)
 

pennylane_lightning/core/src/measurements/tests/Test_MeasurementsBase.cpp

@@ -162,12 +162,15 @@ template <typename TypeList> void testNamedObsExpval() {
         Measurements<StateVectorT> Measurer(statevector);
 
         std::vector<std::vector<size_t>> wires_list = {{0}, {1}, {2}};
-        std::vector<std::string> obs_name = {"PauliX", "PauliY", "PauliZ"};
+        std::vector<std::string> obs_name = {"PauliX", "PauliY", "PauliZ",
+                                             "Hadamard", "Identity"};
         // Expected results calculated with Pennylane default.qubit:
         std::vector<std::vector<PrecisionT>> exp_values_ref = {
             {0.49272486, 0.42073549, 0.28232124},
             {-0.64421768, -0.47942553, -0.29552020},
-            {0.58498357, 0.77015115, 0.91266780}};
+            {0.58498357, 0.77015115, 0.91266780},
+            {0.76205494, 0.84208402, 0.84498485},
+            {1.0, 1.0, 1.0}};
 
         for (size_t ind_obs = 0; ind_obs < obs_name.size(); ind_obs++) {
             DYNAMIC_SECTION(obs_name[ind_obs]

pennylane_lightning/core/src/simulators/lightning_qubit/measurements/ExpValFunctorsLQubit.hpp

@@ -0,0 +1,219 @@
+// Copyright 2018-2023 Xanadu Quantum Technologies Inc.
+
+// Licensed under the Apache License, Version 2.0 (the License);
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+
+// http://www.apache.org/licenses/LICENSE-2.0
+
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an AS IS BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+#pragma once
+
+#include "BitUtil.hpp"
+#include "Util.hpp"
+
+/// @cond DEV
+namespace {
+using namespace Pennylane::Util;
+} // namespace
+/// @endcond
+
+namespace Pennylane::LightningQubit::Functors {
+
+template <class PrecisionT> struct getExpectationValueIdentityFunctor {
+    const std::complex<PrecisionT> *arr;
+    size_t num_qubits;
+
+    getExpectationValueIdentityFunctor(
+        const std::complex<PrecisionT> *arr_,
+        [[maybe_unused]] std::size_t num_qubits_,
+        [[maybe_unused]] const std::vector<size_t> &wires) {
+        arr = arr_;
+        num_qubits = num_qubits_;
+    }
+
+    inline void operator()(PrecisionT &expval) const {
+        size_t k;
+#if defined(_OPENMP)
+#pragma omp parallel for default(none) shared(num_qubits, arr) private(k)      \
+    reduction(+ : expval)
+#endif
+        for (k = 0; k < exp2(num_qubits - 1); k++) {
+            expval += real(conj(arr[k] * arr[k]));
+        }
+    }
+};
+
+template <class PrecisionT> struct getExpectationValuePauliXFunctor {
+    const std::complex<PrecisionT> *arr;
+    size_t num_qubits;
+
+    size_t rev_wire;
+    size_t rev_wire_shift;
+    size_t wire_parity;
+    size_t wire_parity_inv;
+
+    getExpectationValuePauliXFunctor(const std::complex<PrecisionT> *arr_,
+                                     std::size_t num_qubits_,
+                                     const std::vector<size_t> &wires) {
+        arr = arr_;
+        num_qubits = num_qubits_;
+        rev_wire = num_qubits - wires[0] - 1;
+        rev_wire_shift = (static_cast<size_t>(1U) << rev_wire);
+        wire_parity = fillTrailingOnes(rev_wire);
+        wire_parity_inv = fillLeadingOnes(rev_wire + 1);
+    }
+
+    inline void operator()(PrecisionT &expval) const {
+        size_t k;
+        size_t i0;
+        size_t i1;
+#if defined(_OPENMP)
+#pragma omp parallel for default(none)                                         \
+    shared(num_qubits, wire_parity_inv, wire_parity, rev_wire_shift,           \
+               arr) private(k, i0, i1) reduction(+ : expval)
+#endif
+        for (k = 0; k < exp2(num_qubits - 1); k++) {
+            i0 = ((k << 1U) & wire_parity_inv) | (wire_parity & k);
+            i1 = i0 | rev_wire_shift;
+
+            expval +=
+                real(conj(arr[i0]) * arr[i1]) + real(conj(arr[i1]) * arr[i0]);
+        }
+    }
+};
+
+template <class PrecisionT> struct getExpectationValuePauliYFunctor {
+    const std::complex<PrecisionT> *arr;
+    size_t num_qubits;
+
+    size_t rev_wire;
+    size_t rev_wire_shift;
+    size_t wire_parity;
+    size_t wire_parity_inv;
+
+    getExpectationValuePauliYFunctor(const std::complex<PrecisionT> *arr_,
+                                     std::size_t num_qubits_,
+                                     const std::vector<size_t> &wires) {
+        arr = arr_;
+        num_qubits = num_qubits_;
+        rev_wire = num_qubits - wires[0] - 1;
+        rev_wire_shift = (static_cast<size_t>(1U) << rev_wire);
+        wire_parity = fillTrailingOnes(rev_wire);
+        wire_parity_inv = fillLeadingOnes(rev_wire + 1);
+    }
+
+    inline void operator()(PrecisionT &expval) const {
+        size_t k;
+        size_t i0;
+        size_t i1;
+        std::complex<PrecisionT> v0;
+        std::complex<PrecisionT> v1;
+#if defined(_OPENMP)
+#pragma omp parallel for default(none)                                         \
+    shared(num_qubits, wire_parity_inv, wire_parity, rev_wire_shift,           \
+               arr) private(k, i0, i1, v0, v1) reduction(+ : expval)
+#endif
+        for (k = 0; k < exp2(num_qubits - 1); k++) {
+            i0 = ((k << 1U) & wire_parity_inv) | (wire_parity & k);
+            i1 = i0 | rev_wire_shift;
+            v0 = arr[i0];
+            v1 = arr[i1];
+
+            expval += real(conj(arr[i0]) *
+                           std::complex<PrecisionT>{imag(v1), -real(v1)}) +
+                      real(conj(arr[i1]) *
+                           std::complex<PrecisionT>{-imag(v0), real(v0)});
+        }
+    }
+};
+
+template <class PrecisionT> struct getExpectationValuePauliZFunctor {
+    const std::complex<PrecisionT> *arr;
+    size_t num_qubits;
+
+    size_t rev_wire;
+    size_t rev_wire_shift;
+    size_t wire_parity;
+    size_t wire_parity_inv;
+
+    getExpectationValuePauliZFunctor(const std::complex<PrecisionT> *arr_,
+                                     std::size_t num_qubits_,
+                                     const std::vector<size_t> &wires) {
+        arr = arr_;
+        num_qubits = num_qubits_;
+        rev_wire = num_qubits - wires[0] - 1;
+        rev_wire_shift = (static_cast<size_t>(1U) << rev_wire);
+        wire_parity = fillTrailingOnes(rev_wire);
+        wire_parity_inv = fillLeadingOnes(rev_wire + 1);
+    }
+
+    inline void operator()(PrecisionT &expval) const {
+        size_t k;
+        size_t i0;
+        size_t i1;
+#if defined(_OPENMP)
+#pragma omp parallel for default(none)                                         \
+    shared(num_qubits, wire_parity_inv, wire_parity, rev_wire_shift,           \
+               arr) private(k, i0, i1) reduction(+ : expval)
+#endif
+        for (k = 0; k < exp2(num_qubits - 1); k++) {
+            i0 = ((k << 1U) & wire_parity_inv) | (wire_parity & k);
+            i1 = i0 | rev_wire_shift;
+
+            expval += real(conj(arr[i1]) * (-arr[i1])) +
+                      real(conj(arr[i0]) * (arr[i0]));
+        }
+    }
+};
+
+template <class PrecisionT> struct getExpectationValueHadamardFunctor {
+    const std::complex<PrecisionT> *arr;
+    size_t num_qubits;
+
+    size_t rev_wire;
+    size_t rev_wire_shift;
+    size_t wire_parity;
+    size_t wire_parity_inv;
+
+    PrecisionT isqrt2 = INVSQRT2<PrecisionT>();
+
+    getExpectationValueHadamardFunctor(const std::complex<PrecisionT> *arr_,
+                                       std::size_t num_qubits_,
+                                       const std::vector<size_t> &wires) {
+        arr = arr_;
+        num_qubits = num_qubits_;
+        rev_wire = num_qubits - wires[0] - 1;
+        rev_wire_shift = (static_cast<size_t>(1U) << rev_wire);
+        wire_parity = fillTrailingOnes(rev_wire);
+        wire_parity_inv = fillLeadingOnes(rev_wire + 1);
+    }
+
+    inline void operator()(PrecisionT &expval) const {
+        size_t k;
+        size_t i0;
+        size_t i1;
+        std::complex<PrecisionT> v0;
+        std::complex<PrecisionT> v1;
+#if defined(_OPENMP)
+#pragma omp parallel for default(none)                                         \
+    shared(num_qubits, wire_parity_inv, wire_parity, rev_wire_shift,           \
+               arr) private(k, i0, i1, v0, v1) reduction(+ : expval)
+#endif
+
+        for (k = 0; k < exp2(num_qubits - 1); k++) {
+            i0 = ((k << 1U) & wire_parity_inv) | (wire_parity & k);
+            i1 = i0 | rev_wire_shift;
+            v0 = arr[i0];
+            v1 = arr[i1];
+
+            expval += real(isqrt2 * (conj(arr[i0]) * (v0 + v1) +
+                                     conj(arr[i1]) * (v0 - v1)));
+        }
+    }
+};
+} // namespace Pennylane::LightningQubit::Functors

pennylane_lightning/core/src/simulators/lightning_qubit/measurements/MeasurementsLQubit.hpp

@@ -29,6 +29,7 @@
 #include <unordered_map>
 #include <vector>
 
+#include "ExpValFunctorsLQubit.hpp"
 #include "LinearAlgebra.hpp"
 #include "MeasurementsBase.hpp"
 #include "Observables.hpp"
@@ -44,6 +45,17 @@
 using namespace Pennylane::Observables;
 using Pennylane::LightningQubit::StateVectorLQubitManaged;
 using Pennylane::LightningQubit::Util::innerProdC;
+using namespace Pennylane::LightningQubit::Functors;
+using Pennylane::Util::INVSQRT2;
+enum class ExpValFunc : uint32_t {
+    BEGIN = 1,
+    Identity = 1,
+    PauliX,
+    PauliY,
+    PauliZ,
+    Hadamard,
+    END
+};
 } // namespace
 /// @endcond
 
@@ -67,7 +79,35 @@
 
   public:
     explicit Measurements(const StateVectorT &statevector)
-        : BaseType{statevector} {};
+        : BaseType{statevector} {
+        init_expval_funcs_();
+    };
+
+    /**
+     * @brief Templated method that returns the expectation value of named
+     * observables using in-place operations, without creating extra copy of
+     * the statevector.
+     *
+     * @tparam functor_t Expectation value functor class for in-place
+     * operations.
+     * @tparam nqubits Number of wires.
+     * @param wires Wires to which the observable is applied.
+     * @return expectation value of the observable.
+     */
+    template <template <class> class functor_t, size_t num_wires>
+    PrecisionT applyExpValNamedFunctor(const std::vector<size_t> &wires) {
+        if constexpr (num_wires > 0) {
+            PL_ASSERT(wires.size() == num_wires);
+        }
+
+        size_t num_qubits = this->_statevector.getNumQubits();
+        const std::complex<PrecisionT> *arr_data = this->_statevector.getData();
+
+        PrecisionT expval = 0.0;
+        functor_t<PrecisionT>(arr_data, num_qubits, wires)(expval);
+
+        return expval;
+    }
 
     /**
      * @brief Probabilities of each computational basis state.
@@ -165,17 +205,30 @@
      */
     PrecisionT expval(const std::string &operation,
                       const std::vector<size_t> &wires) {
-        // Copying the original state vector, for the application of the
-        // observable operator.
-        StateVectorLQubitManaged<PrecisionT> operator_statevector(
-            this->_statevector);
-
-        operator_statevector.applyOperation(operation, wires);
-
-        ComplexT expected_value = innerProdC(this->_statevector.getData(),
-                                             operator_statevector.getData(),
-                                             this->_statevector.getLength());
-        return std::real(expected_value);
+        // In-place calculation of expval without creating duplicate of the
+        // statevector.
+        switch (expval_funcs_[operation]) {
+        case ExpValFunc::Identity:
+            return applyExpValNamedFunctor<getExpectationValueIdentityFunctor,
+                                           0>(wires);
+        case ExpValFunc::PauliX:
+            return applyExpValNamedFunctor<getExpectationValuePauliXFunctor, 1>(
+                wires);
+        case ExpValFunc::PauliY:
+            return applyExpValNamedFunctor<getExpectationValuePauliYFunctor, 1>(
+                wires);
+        case ExpValFunc::PauliZ:
+            return applyExpValNamedFunctor<getExpectationValuePauliZFunctor, 1>(
+                wires);
+        case ExpValFunc::Hadamard:
+            return applyExpValNamedFunctor<getExpectationValueHadamardFunctor,
+                                           1>(wires);
+        default:
+            PL_ABORT(
+                std::string("Expval does not exist for named observable ") +
+                operation);
+            break;
+        }
     };
 
     /**
@@ -661,5 +714,21 @@
         }
         return init_idx;
     }
+
+    std::unordered_map<std::string, ExpValFunc> expval_funcs_;
+
+    // clang-format off
+        /**
+        * @brief Register generator operations in the generators_indices_ attribute:
+        *        an unordered_map mapping strings to GateOperation enumeration keywords.
+        */
+        void init_expval_funcs_() {
+            expval_funcs_["Identity"] = ExpValFunc::Identity;
+            expval_funcs_["PauliX"]   = ExpValFunc::PauliX;
+            expval_funcs_["PauliY"]   = ExpValFunc::PauliY;
+            expval_funcs_["PauliZ"]   = ExpValFunc::PauliZ;
+            expval_funcs_["Hadamard"] = ExpValFunc::Hadamard;
+        }
+    // clang-format on
 }; // class Measurements
 } // namespace Pennylane::LightningQubit::Measures