diff --git a/runtime/common/NoiseModel.cpp b/runtime/common/NoiseModel.cpp index aa276f68284..3ecdc7cc70a 100644 --- a/runtime/common/NoiseModel.cpp +++ b/runtime/common/NoiseModel.cpp @@ -10,8 +10,73 @@ #include "Logger.h" #include "common/CustomOp.h" #include "common/EigenDense.h" +#include +#include + namespace cudaq { +// Helper to check whether a matrix is a scaled unitary matrix, i.e., `k * U` +// where U is a unitary matrix. If so, it also returns the `k` factor. +// Otherwise, return a nullopt. +static std::optional +isScaledUnitary(const std::vector> &mat, double eps) { + typedef Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic, + Eigen::RowMajor> + RowMajorMatTy; + const int dim = std::log2(mat.size()); + Eigen::Map kMat(mat.data(), dim, dim); + if (kMat.isZero(eps)) + return 0.0; + // Check that (K_dag * K) is a scaled identity matrix + // i.e., the K matrix is a scaled unitary. + auto kdK = kMat.adjoint() * kMat; + if (!kdK.isDiagonal(eps)) + return std::nullopt; + // First element + std::complex val = kdK(0, 0); + if (val.real() > eps && std::abs(val.imag()) < eps) { + auto scaledKdK = (std::complex{1.0} / val) * kdK; + if (scaledKdK.isIdentity()) + return std::sqrt(val.real()); + } + return std::nullopt; +} + +// Helper to determine if a vector of Kraus ops are actually a unitary mixture. +// If so, it returns all the unitaries and the probabilities associated with +// each one of those unitaries. +static std::optional, + std::vector>>>> +computeUnitaryMixture( + const std::vector>> &krausOps, + double tol = 1e-6) { + std::vector probs; + std::vector>> mats; + const auto scaleMat = [](const std::vector> &mat, + double scaleFactor) { + std::vector> scaledMat = mat; + // If scaleFactor is 0, then it means the original matrix was likely all + // zeros. In that case, the probability will be 0, so the matrix doesn't + // matter, but we don't want NaNs to trickle in anywhere. + if (scaleFactor != 0.0) + for (auto &x : scaledMat) + x /= scaleFactor; + return scaledMat; + }; + for (const auto &op : krausOps) { + const auto scaledFactor = isScaledUnitary(op, tol); + if (!scaledFactor.has_value()) + return std::nullopt; + probs.emplace_back(scaledFactor.value() * scaledFactor.value()); + mats.emplace_back(scaleMat(op, scaledFactor.value())); + } + + if (std::abs(1.0 - std::reduce(probs.begin(), probs.end())) > tol) + return std::nullopt; + + return std::make_pair(probs, mats); +} + template bool isIdentity(const EigenMatTy &mat, double threshold = 1e-9) { EigenMatTy idMat = EigenMatTy::Identity(mat.rows(), mat.cols()); @@ -78,9 +143,52 @@ void validateCompletenessRelation_fp64(const std::vector &ops) { "Provided kraus_ops are not completely positive and trace preserving."); } +void generateUnitaryParameters_fp32( + const std::vector &ops, + std::vector>> &unitary_ops, + std::vector &probabilities) { + std::vector>> double_kraus_ops; + double_kraus_ops.reserve(ops.size()); + for (auto &op : ops) { + // WARNING: danger here. We are intentially treating the incoming op as fp32 + // type instead of what the compiler thinks it is (fp64). We have to do this + // because this file is compiled with cudaq::real = fp64, but the incoming + // data for this specific routine is actually fp32. + const std::complex *ptr = + reinterpret_cast *>(op.data.data()); + // Use 2 * size because pointer arithmetic is on fp32 instead of fp64 + double_kraus_ops.emplace_back( + std::vector>(ptr, ptr + 2 * op.data.size())); + } + + auto asUnitaryMixture = computeUnitaryMixture(double_kraus_ops); + if (asUnitaryMixture.has_value()) { + probabilities = std::move(asUnitaryMixture.value().first); + unitary_ops = std::move(asUnitaryMixture.value().second); + } +} + +void generateUnitaryParameters_fp64( + const std::vector &ops, + std::vector>> &unitary_ops, + std::vector &probabilities) { + std::vector>> double_kraus_ops; + double_kraus_ops.reserve(ops.size()); + for (auto &op : ops) + double_kraus_ops.emplace_back( + std::vector>(op.data.begin(), op.data.end())); + + auto asUnitaryMixture = computeUnitaryMixture(double_kraus_ops); + if (asUnitaryMixture.has_value()) { + probabilities = std::move(asUnitaryMixture.value().first); + unitary_ops = std::move(asUnitaryMixture.value().second); + } +} + kraus_channel::kraus_channel(const kraus_channel &other) : ops(other.ops), noise_type(other.noise_type), - parameters(other.parameters) {} + parameters(other.parameters), unitary_ops(other.unitary_ops), + probabilities(other.probabilities) {} std::size_t kraus_channel::size() const { return ops.size(); } @@ -94,6 +202,8 @@ kraus_channel &kraus_channel::operator=(const kraus_channel &other) { ops = other.ops; noise_type = other.noise_type; parameters = other.parameters; + unitary_ops = other.unitary_ops; + probabilities = other.probabilities; return *this; } diff --git a/runtime/common/NoiseModel.h b/runtime/common/NoiseModel.h index 6ee2bc5e3da..bdf0ff9c28a 100644 --- a/runtime/common/NoiseModel.h +++ b/runtime/common/NoiseModel.h @@ -102,6 +102,12 @@ struct kraus_op { void validateCompletenessRelation_fp32(const std::vector &ops); void validateCompletenessRelation_fp64(const std::vector &ops); +void generateUnitaryParameters_fp32( + const std::vector &ops, + std::vector>> &, std::vector &); +void generateUnitaryParameters_fp64( + const std::vector &ops, + std::vector>> &, std::vector &); /// @brief A kraus_channel represents a quantum noise channel /// on specific qubits. The action of the noise channel is @@ -143,7 +149,17 @@ class kraus_channel { // corruption. std::vector parameters; - ~kraus_channel() = default; + /// @brief If all Kraus ops are - when scaled - unitary, this holds the + /// unitary versions of those ops. These values are always "double" regardless + /// of whether cudaq::real is float or double. + std::vector>> unitary_ops; + + /// @brief If all Kraus ops are - when scaled - unitary, this holds the + /// probabilities of those ops. These values are always "double" regardless + /// of whether cudaq::real is float or double. + std::vector probabilities; + + virtual ~kraus_channel() = default; /// @brief The nullary constructor kraus_channel() = default; @@ -158,12 +174,14 @@ class kraus_channel { kraus_channel(std::initializer_list &&...inputLists) { (ops.emplace_back(std::move(inputLists)), ...); validateCompleteness(); + generateUnitaryParameters(); } /// @brief The constructor, take qubits and channel kraus_ops as lvalue /// reference kraus_channel(const std::vector &inOps) : ops(inOps) { validateCompleteness(); + generateUnitaryParameters(); } /// @brief The constructor, take qubits and channel kraus_ops as rvalue @@ -189,6 +207,23 @@ class kraus_channel { /// @brief Add a kraus_op to this channel. void push_back(kraus_op op); + + /// @brief Returns whether or not this is a unitary mixture. + bool is_unitary_mixture() const { return !unitary_ops.empty(); } + + /// @brief Checks if Kraus ops have unitary representations and saves them if + /// they do. Users should only need to call this if they have modified the + /// Kraus ops and want to recompute these values. + void generateUnitaryParameters() { + unitary_ops.clear(); + probabilities.clear(); + if constexpr (std::is_same_v) { + generateUnitaryParameters_fp32(ops, this->unitary_ops, + this->probabilities); + return; + } + generateUnitaryParameters_fp64(ops, this->unitary_ops, this->probabilities); + } }; /// @brief The noise_model type keeps track of a set of @@ -381,6 +416,7 @@ class depolarization_channel : public kraus_channel { this->parameters.push_back(probability); noise_type = noise_model_type::depolarization_channel; validateCompleteness(); + generateUnitaryParameters(); } }; @@ -396,6 +432,8 @@ class amplitude_damping_channel : public kraus_channel { this->parameters.push_back(probability); noise_type = noise_model_type::amplitude_damping_channel; validateCompleteness(); + // Note: amplitude damping is non-unitary, so there is no value in calling + // generateUnitaryParameters(). } }; @@ -412,6 +450,7 @@ class bit_flip_channel : public kraus_channel { this->parameters.push_back(probability); noise_type = noise_model_type::bit_flip_channel; validateCompleteness(); + generateUnitaryParameters(); } }; @@ -429,6 +468,7 @@ class phase_flip_channel : public kraus_channel { this->parameters.push_back(probability); noise_type = noise_model_type::phase_flip_channel; validateCompleteness(); + generateUnitaryParameters(); } }; } // namespace cudaq diff --git a/runtime/nvqir/cutensornet/simulator_cutensornet.cpp b/runtime/nvqir/cutensornet/simulator_cutensornet.cpp index a4907ccc9e9..495a6126f44 100644 --- a/runtime/nvqir/cutensornet/simulator_cutensornet.cpp +++ b/runtime/nvqir/cutensornet/simulator_cutensornet.cpp @@ -11,17 +11,6 @@ #include "cutensornet.h" #include "tensornet_spin_op.h" -namespace { -const std::vector> matPauliI = { - {1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {1.0, 0.0}}; -const std::vector> matPauliX{ - {0.0, 0.0}, {1.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}}; -const std::vector> matPauliY{ - {0.0, 0.0}, {0.0, -1.0}, {0.0, 1.0}, {0.0, 0.0}}; -const std::vector> matPauliZ{ - {1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {-1.0, 0.0}}; -} // namespace - namespace nvqir { SimulatorTensorNetBase::SimulatorTensorNetBase() @@ -148,66 +137,7 @@ void SimulatorTensorNetBase::applyGate(const GateApplicationTask &task) { } } -// Helper to check whether a matrix is a scaled unitary matrix, i.e., `k * U` -// where U is a unitary matrix. If so, it also returns the `k` factor. -// Otherwise, return a nullopt. -template -std::optional isScaledUnitary(const std::vector> &mat, - double eps) { - typedef Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic, - Eigen::RowMajor> - RowMajorMatTy; - const int dim = std::log2(mat.size()); - Eigen::Map kMat(mat.data(), dim, dim); - if (kMat.isZero()) - return std::nullopt; - // Check that (K_dag * K) is a scaled identity matrix - // i.e., the K matrix is a scaled unitary. - auto kdK = kMat.adjoint() * kMat; - if (!kdK.isDiagonal()) - return std::nullopt; - // First element - std::complex val = kdK(0, 0); - if (std::abs(val) > eps && std::abs(val.imag()) < eps) { - auto scaledKdK = (std::complex{1.0} / val) * kdK; - if (scaledKdK.isIdentity()) - return val.real(); - } - return std::nullopt; -} - -std::optional, - std::vector>>>> static computeUnitaryMixture(const std:: - vector>> - &krausOps, - double tol = 1e-6) { - std::vector probs; - std::vector>> mats; - const auto scaleMat = [](const std::vector> &mat, - double scaleFactor) { - std::vector> scaledMat = mat; - for (auto &x : scaledMat) - x /= scaleFactor; - return scaledMat; - }; - for (const auto &op : krausOps) { - const auto scaledFactor = isScaledUnitary(op, tol); - if (!scaledFactor.has_value()) - return std::nullopt; - probs.emplace_back(scaledFactor.value()); - mats.emplace_back(scaleMat(op, scaledFactor.value())); - } - - if (std::abs(1.0 - std::reduce(probs.begin(), probs.end())) > tol) - return std::nullopt; - - return std::make_pair(probs, mats); -} - -// Helper to look up a device memory pointer from a cache. +// Helper to look up a device memory pointer from a cache. // If not found, allocate a new device memory buffer and put it to the cache. static void * getOrCacheMat(const std::string &key, @@ -227,85 +157,16 @@ void SimulatorTensorNetBase::applyKrausChannel( const std::vector &qubits, const cudaq::kraus_channel &krausChannel) { LOG_API_TIME(); - switch (krausChannel.noise_type) { - case cudaq::noise_model_type::depolarization_channel: { - if (krausChannel.parameters.size() != 1) - throw std::runtime_error( - fmt::format("Invalid parameters for a depolarization channel. " - "Expecting 1 parameter, got {}.", - krausChannel.parameters.size())); - const std::vector channelMats{ - getOrCacheMat("PauliI", matPauliI, m_gateDeviceMemCache), - getOrCacheMat("PauliX", matPauliX, m_gateDeviceMemCache), - getOrCacheMat("PauliY", matPauliY, m_gateDeviceMemCache), - getOrCacheMat("PauliZ", matPauliZ, m_gateDeviceMemCache)}; - const double p = krausChannel.parameters[0]; - const std::vector probabilities = {1 - p, p / 3., p / 3., p / 3.}; - m_state->applyUnitaryChannel(qubits, channelMats, probabilities); - break; - } - case cudaq::noise_model_type::bit_flip_channel: { - if (krausChannel.parameters.size() != 1) - throw std::runtime_error( - fmt::format("Invalid parameters for a bit-flip channel. " - "Expecting 1 parameter, got {}.", - krausChannel.parameters.size())); - - const std::vector channelMats{ - getOrCacheMat("PauliI", matPauliI, m_gateDeviceMemCache), - getOrCacheMat("PauliX", matPauliX, m_gateDeviceMemCache)}; - const double p = krausChannel.parameters[0]; - const std::vector probabilities = {1 - p, p}; - m_state->applyUnitaryChannel(qubits, channelMats, probabilities); - break; - } - case cudaq::noise_model_type::phase_flip_channel: { - if (krausChannel.parameters.size() != 1) - throw std::runtime_error( - fmt::format("Invalid parameters for a phase-flip channel. " - "Expecting 1 parameter, got {}.", - krausChannel.parameters.size())); - - const std::vector channelMats{ - getOrCacheMat("PauliI", matPauliI, m_gateDeviceMemCache), - getOrCacheMat("PauliZ", matPauliZ, m_gateDeviceMemCache)}; - const double p = krausChannel.parameters[0]; - const std::vector probabilities = {1 - p, p}; - m_state->applyUnitaryChannel(qubits, channelMats, probabilities); - break; - } - case cudaq::noise_model_type::amplitude_damping_channel: { - if (krausChannel.parameters.size() != 1) - throw std::runtime_error( - fmt::format("Invalid parameters for a amplitude damping channel. " - "Expecting 1 parameter, got {}.", - krausChannel.parameters.size())); - if (krausChannel.parameters[0] != 0.0) - throw std::runtime_error("Non-unitary noise channels are not supported."); - break; - } - case cudaq::noise_model_type::unknown: { - std::vector>> mats; - for (const auto &op : krausChannel.get_ops()) - mats.emplace_back(op.data); - auto asUnitaryMixture = computeUnitaryMixture(mats); - if (asUnitaryMixture.has_value()) { - auto &[probabilities, unitaries] = asUnitaryMixture.value(); - std::vector channelMats; - for (const auto &mat : unitaries) - channelMats.emplace_back(getOrCacheMat( - "ScaledUnitary_" + std::to_string(vecComplexHash(mat)), mat, - m_gateDeviceMemCache)); - m_state->applyUnitaryChannel(qubits, channelMats, probabilities); - } else { - throw std::runtime_error("Non-unitary noise channels are not supported."); - } - break; - } - default: - throw std::runtime_error( - "Unsupported noise model type: " + - std::to_string(static_cast(krausChannel.noise_type))); + if (krausChannel.is_unitary_mixture()) { + std::vector channelMats; + for (const auto &mat : krausChannel.unitary_ops) + channelMats.emplace_back( + getOrCacheMat("ScaledUnitary_" + std::to_string(vecComplexHash(mat)), + mat, m_gateDeviceMemCache)); + m_state->applyUnitaryChannel(qubits, channelMats, + krausChannel.probabilities); + } else { + throw std::runtime_error("Non-unitary noise channels are not supported."); } } diff --git a/unittests/integration/noise_tester.cpp b/unittests/integration/noise_tester.cpp index 8368f7a6a01..409d454b616 100644 --- a/unittests/integration/noise_tester.cpp +++ b/unittests/integration/noise_tester.cpp @@ -378,71 +378,71 @@ CUDAQ_TEST(NoiseTest, checkAllQubitChannel) { // Stim does not support arbitrary cudaq::kraus_op specification. static cudaq::kraus_channel create2pNoiseChannel() { - // 1% depolarization - cudaq::kraus_op op0{cudaq::complex{0.94868329805, 0.0}, + // 20% depolarization + cudaq::kraus_op op0{cudaq::complex{0.894427191, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, - {0.94868329805, 0.0}, + {0.894427191, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, - {0.94868329805, 0.0}, + {0.894427191, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, - {0.94868329805, 0.0}}, + {0.894427191, 0.0}}, op1{cudaq::complex{0.0, 0.0}, {0.0, 0.0}, - {0.18257418583, 0.0}, + {0.25819888974, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, - {0.18257418583, 0.0}, - {0.18257418583, 0.0}, + {0.25819888974, 0.0}, + {0.25819888974, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, - {0.18257418583, 0.0}, + {0.25819888974, 0.0}, {0.0, 0.0}, {0.0, 0.0}}, op2{cudaq::complex{0.0, 0.0}, {0.0, 0.0}, - {0.0, -0.18257418583}, + {0.0, -0.25819888974}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, - {0.0, -0.18257418583}, - {0.0, 0.18257418583}, + {0.0, -0.25819888974}, + {0.0, 0.25819888974}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, - {0.0, 0.18257418583}, + {0.0, 0.25819888974}, {0.0, 0.0}, {0.0, 0.0}}, - op3{cudaq::complex{0.18257418583, 0.0}, + op3{cudaq::complex{0.25819888974, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, - {0.18257418583, 0.0}, + {0.25819888974, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, - {-0.18257418583, 0.0}, + {-0.25819888974, 0.0}, {-0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {-0.0, 0.0}, - {-0.18257418583, 0.0}}; + {-0.25819888974, 0.0}}; cudaq::kraus_channel noise2q( std::vector{op0, op1, op2, op3}); return noise2q; @@ -481,6 +481,7 @@ CUDAQ_TEST(NoiseTest, checkAllQubitChannelWithControl) { auto counts = cudaq::sample({.shots = shots, .noise = noise}, bellRandom{}, qubitIds[0], qubitIds[1]); + counts.dump(); // More than 2 entangled states due to the noise. EXPECT_GT(counts.size(), 2); } while (std::next_permutation(qubitIds.begin(), qubitIds.end())); @@ -509,6 +510,7 @@ CUDAQ_TEST(NoiseTest, checkAllQubitChannelWithControlPrefix) { auto counts = cudaq::sample({.shots = shots, .noise = noise}, bellRandom{}, qubitIds[0], qubitIds[1]); + counts.dump(); // More than 2 entangled states due to the noise. EXPECT_GT(counts.size(), 2); } while (std::next_permutation(qubitIds.begin(), qubitIds.end()));