QuantumChannel.h

Go to the documentation of this file.

 
#pragma once
 
#ifndef _SIMULATORS_QUANTUM_CHANNEL_H_
#define _SIMULATORS_QUANTUM_CHANNEL_H_
 
#include <Eigen/Eigen>
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <cstddef>
#include <limits>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
 
namespace Simulators
{

class QuantumChannel
{
  public:
    using Matrix = Eigen::MatrixXcd;
    using KrausOperators = std::vector<Matrix>;

    explicit QuantumChannel(KrausOperators krausOperators) : krausOperators_(std::move(krausOperators)), numQubits_(0)
    {
        Validate();
    }
 
    const KrausOperators &GetKrausOperators() const
    {
        return krausOperators_;
    }
 
    size_t GetNumberOfQubits() const
    {
        return numQubits_;
    }

    const std::vector<double> &GetPauliProbabilities() const
    {
        return pauliProbabilities_;
    }

    bool IsApprox(const QuantumChannel &other, double tolerance = 1e-12) const
    {
        if (tolerance < 0.0 || !std::isfinite(tolerance))
            throw std::invalid_argument("Quantum-channel comparison tolerance must be finite and nonnegative");
        if (numQubits_ != other.numQubits_)
            return false;
 
        const Matrix left = ChoiMatrix();
        const Matrix right = other.ChoiMatrix();
        if (left.rows() != right.rows() || left.cols() != right.cols())
            return false;
        return (left - right).cwiseAbs().maxCoeff() <= tolerance;
    }

    Matrix ChoiMatrix() const
    {
        const Eigen::Index dimension = krausOperators_.front().rows();
        const Eigen::Index choiDimension = dimension * dimension;
        Matrix choi = Matrix::Zero(choiDimension, choiDimension);
        for (const Matrix &krausOperator : krausOperators_)
        {
            const Eigen::Map<const Eigen::VectorXcd> vectorized(krausOperator.data(), choiDimension);
            choi.noalias() += vectorized * vectorized.adjoint();
        }
        return choi;
    }

    static QuantumChannel BitFlip(double probability)
    {
        return Pauli({1.0 - CheckedProbability(probability), probability, 0.0, 0.0});
    }

    static QuantumChannel BitPhaseFlip(double probability)
    {
        return Pauli({1.0 - CheckedProbability(probability), 0.0, probability, 0.0});
    }

    static QuantumChannel PhaseFlip(double probability)
    {
        return Pauli({1.0 - CheckedProbability(probability), 0.0, 0.0, probability});
    }

    static QuantumChannel Pauli(double px, double py, double pz)
    {
        CheckedProbability(px);
        CheckedProbability(py);
        CheckedProbability(pz);
        const double total = px + py + pz;
        if (!std::isfinite(total) || total > 1.0 + kProbabilityTolerance)
            throw std::invalid_argument("Pauli error probabilities must sum to at most one");
        return Pauli({std::max(0.0, 1.0 - total), px, py, pz});
    }

    static QuantumChannel Pauli(const std::vector<double> &probabilities)
    {
        if (probabilities.empty())
            throw std::invalid_argument("A Pauli channel must contain probabilities");
 
        size_t count = probabilities.size();
        size_t numQubits = 0;
        while (count > 1 && count % 4 == 0)
        {
            count /= 4;
            ++numQubits;
        }
        if (count != 1 || numQubits == 0)
            throw std::invalid_argument("A Pauli channel needs exactly 4^n probabilities for n >= 1");
 
        double total = 0.0;
        for (const double probability : probabilities)
        {
            CheckedProbability(probability);
            total += probability;
        }
        if (!std::isfinite(total) || std::abs(total - 1.0) > kProbabilityTolerance)
            throw std::invalid_argument("Pauli channel probabilities must sum to one");
 
        KrausOperators kraus;
        kraus.reserve(static_cast<size_t>(std::count_if(probabilities.begin(), probabilities.end(), [](double probability) { return probability > 0.0; })));
        for (size_t index = 0; index < probabilities.size(); ++index)
        {
            if (probabilities[index] == 0.0)
                continue;
 
            Matrix pauli(1, 1);
            pauli(0, 0) = 1.0;
            size_t divisor = 1;
            for (size_t q = 0; q < numQubits; ++q)
                divisor *= 4;
            for (size_t q = numQubits; q-- > 0;)
            {
                divisor /= 4;
                const size_t pauliId = (index / divisor) % 4;
                pauli = Kronecker(pauli, PauliMatrix(pauliId));
            }
            kraus.emplace_back(std::sqrt(probabilities[index]) * pauli);
        }
        QuantumChannel channel(std::move(kraus));
        channel.pauliProbabilities_ = probabilities;
        return channel;
    }

    static QuantumChannel Depolarizing(double errorProbability)
    {
        CheckedProbability(errorProbability);
        return Pauli({1.0 - errorProbability, errorProbability / 3.0, errorProbability / 3.0, errorProbability / 3.0});
    }

    static QuantumChannel DepolarizingMixing(double mixingProbability)
    {
        CheckedProbability(mixingProbability);
        return Depolarizing(0.75 * mixingProbability);
    }

    static QuantumChannel AmplitudeDamping(double gamma)
    {
        CheckedProbability(gamma);
        Matrix k0 = Matrix::Zero(2, 2);
        k0(0, 0) = 1.0;
        k0(1, 1) = std::sqrt(1.0 - gamma);
        Matrix k1 = Matrix::Zero(2, 2);
        k1(0, 1) = std::sqrt(gamma);
        return QuantumChannel({std::move(k0), std::move(k1)});
    }

    static QuantumChannel PhaseDamping(double gamma)
    {
        CheckedProbability(gamma);
        Matrix k0 = Matrix::Zero(2, 2);
        k0(0, 0) = 1.0;
        k0(1, 1) = std::sqrt(1.0 - gamma);
        Matrix k1 = Matrix::Zero(2, 2);
        k1(1, 1) = std::sqrt(gamma);
        return QuantumChannel({std::move(k0), std::move(k1)});
    }

    static QuantumChannel GeneralizedAmplitudeDamping(double gamma, double excitedStatePopulation)
    {
        CheckedProbability(gamma);
        CheckedProbability(excitedStatePopulation);
        const double groundStatePopulation = 1.0 - excitedStatePopulation;
 
        Matrix k0 = Matrix::Zero(2, 2);
        k0(0, 0) = std::sqrt(groundStatePopulation);
        k0(1, 1) = std::sqrt(groundStatePopulation * (1.0 - gamma));
        Matrix k1 = Matrix::Zero(2, 2);
        k1(0, 1) = std::sqrt(groundStatePopulation * gamma);
        Matrix k2 = Matrix::Zero(2, 2);
        k2(0, 0) = std::sqrt(excitedStatePopulation * (1.0 - gamma));
        k2(1, 1) = std::sqrt(excitedStatePopulation);
        Matrix k3 = Matrix::Zero(2, 2);
        k3(1, 0) = std::sqrt(excitedStatePopulation * gamma);
        return QuantumChannel({std::move(k0), std::move(k1), std::move(k2), std::move(k3)});
    }

    static QuantumChannel ThermalRelaxation(double duration, double t1, double t2, double excitedStatePopulation = 0.0)
    {
        if (!std::isfinite(duration) || duration < 0.0)
            throw std::invalid_argument("Thermal-relaxation duration must be finite and nonnegative");
        ValidatePositiveTime(t1, "T1");
        ValidatePositiveTime(t2, "T2");
        CheckedProbability(excitedStatePopulation);
 
        if (std::isinf(t2) && !std::isinf(t1))
            throw std::invalid_argument("Thermal relaxation requires T2 <= 2*T1");
        if (std::isfinite(t1) && std::isfinite(t2) && t2 > 2.0 * t1 * (1.0 + kProbabilityTolerance))
            throw std::invalid_argument("Thermal relaxation requires T2 <= 2*T1");
 
        const double coherenceMultiplier = std::isinf(t2) ? 1.0 : std::exp(-duration / t2);
        const double gamma = std::isinf(t1) ? 0.0 : -std::expm1(-duration / t1);
        const double up = gamma * excitedStatePopulation;
        const double down = gamma * (1.0 - excitedStatePopulation);
 
        // The Choi matrix has a 2x2 block
        // [[1-up, coherence], [coherence, 1-down]]. Decomposing that block gives
        // at most two diagonal Kraus operators; excitation and decay add one each.
        const double a = 1.0 - up;
        const double c = 1.0 - down;
        if (coherenceMultiplier * coherenceMultiplier > a * c + kProbabilityTolerance)
            throw std::invalid_argument("The supplied T1 and T2 do not define a completely-positive map");
 
        const double discriminant = std::hypot(a - c, 2.0 * coherenceMultiplier);
        const double lambdaPlus = 0.5 * (a + c + discriminant);
        const double lambdaMinus = std::max(0.0, 0.5 * (a + c - discriminant));
        const double angle = 0.5 * std::atan2(2.0 * coherenceMultiplier, a - c);
        const double cosine = std::cos(angle);
        const double sine = std::sin(angle);
 
        KrausOperators kraus;
        if (lambdaPlus > kZeroTolerance)
        {
            Matrix diagonal = Matrix::Zero(2, 2);
            diagonal(0, 0) = std::sqrt(lambdaPlus) * cosine;
            diagonal(1, 1) = std::sqrt(lambdaPlus) * sine;
            kraus.emplace_back(std::move(diagonal));
        }
        if (lambdaMinus > kZeroTolerance)
        {
            Matrix diagonal = Matrix::Zero(2, 2);
            diagonal(0, 0) = -std::sqrt(lambdaMinus) * sine;
            diagonal(1, 1) = std::sqrt(lambdaMinus) * cosine;
            kraus.emplace_back(std::move(diagonal));
        }
        if (up > kZeroTolerance)
        {
            Matrix excitation = Matrix::Zero(2, 2);
            excitation(1, 0) = std::sqrt(up);
            kraus.emplace_back(std::move(excitation));
        }
        if (down > kZeroTolerance)
        {
            Matrix decay = Matrix::Zero(2, 2);
            decay(0, 1) = std::sqrt(down);
            kraus.emplace_back(std::move(decay));
        }
        return QuantumChannel(std::move(kraus));
    }

    static QuantumChannel CorrelatedPhaseFlip(double probability)
    {
        return CorrelatedPhaseFlip(probability, 1.0);
    }

    static QuantumChannel CorrelatedPhaseFlip(double probability, double correlation)
    {
        CheckedProbability(probability);
        CheckedProbability(correlation);
        std::vector<double> probabilities(16, 0.0);
        const double independent = 1.0 - correlation;
        probabilities[0] = independent * (1.0 - probability) * (1.0 - probability) + correlation * (1.0 - probability);
        probabilities[3] = independent * probability * (1.0 - probability);
        probabilities[12] = probabilities[3];
        probabilities[15] = independent * probability * probability + correlation * probability;
        return Pauli(probabilities);
    }

    static QuantumChannel TwoQubitDepolarizing(double errorProbability)
    {
        CheckedProbability(errorProbability);
        std::vector<double> probabilities(16, errorProbability / 15.0);
        probabilities[0] = 1.0 - errorProbability;
        return Pauli(probabilities);
    }

    static QuantumChannel TwoQubitDepolarizingMixing(double mixingProbability)
    {
        CheckedProbability(mixingProbability);
        return TwoQubitDepolarizing((15.0 / 16.0) * mixingProbability);
    }
 
  private:
    static constexpr double kProbabilityTolerance = 1e-12;
    static constexpr double kZeroTolerance = 1e-15;
 
    static double CheckedProbability(double probability)
    {
        if (!std::isfinite(probability) || probability < 0.0 || probability > 1.0)
            throw std::invalid_argument("Quantum-channel probabilities must be finite and in [0, 1]");
        return probability;
    }
 
    static void ValidatePositiveTime(double value, const char *name)
    {
        if (std::isnan(value) || value <= 0.0)
            throw std::invalid_argument(std::string(name) + " must be positive (infinity is allowed)");
    }
 
    static Matrix PauliMatrix(size_t pauliId)
    {
        Matrix matrix = Matrix::Zero(2, 2);
        switch (pauliId)
        {
        case 0:
            matrix.setIdentity();
            break;
        case 1:
            matrix(0, 1) = 1.0;
            matrix(1, 0) = 1.0;
            break;
        case 2:
            matrix(0, 1) = std::complex<double>(0.0, -1.0);
            matrix(1, 0) = std::complex<double>(0.0, 1.0);
            break;
        case 3:
            matrix(0, 0) = 1.0;
            matrix(1, 1) = -1.0;
            break;
        default:
            throw std::invalid_argument("Invalid Pauli index");
        }
        return matrix;
    }
 
    static Matrix Kronecker(const Matrix &left, const Matrix &right)
    {
        Matrix result(left.rows() * right.rows(), left.cols() * right.cols());
        for (Eigen::Index row = 0; row < left.rows(); ++row)
            for (Eigen::Index column = 0; column < left.cols(); ++column)
                result.block(row * right.rows(), column * right.cols(), right.rows(), right.cols()) = left(row, column) * right;
        return result;
    }
 
    void Validate()
    {
        if (krausOperators_.empty())
            throw std::invalid_argument("A quantum channel must contain at least one Kraus operator");
 
        const Eigen::Index dimension = krausOperators_.front().rows();
        if (dimension < 2 || krausOperators_.front().cols() != dimension)
            throw std::invalid_argument("Kraus operators must be nonempty square qubit matrices");
 
        size_t remainingDimension = static_cast<size_t>(dimension);
        while (remainingDimension > 1 && remainingDimension % 2 == 0)
        {
            remainingDimension /= 2;
            ++numQubits_;
        }
        if (remainingDimension != 1)
            throw std::invalid_argument("Kraus-operator dimensions must be powers of two");
 
        Matrix completeness = Matrix::Zero(dimension, dimension);
        for (const Matrix &krausOperator : krausOperators_)
        {
            if (krausOperator.rows() != dimension || krausOperator.cols() != dimension)
                throw std::invalid_argument("All Kraus operators must have the same dimensions");
            if (!krausOperator.allFinite())
                throw std::invalid_argument("Kraus operators must contain only finite values");
            completeness.noalias() += krausOperator.adjoint() * krausOperator;
        }
 
        const Matrix identity = Matrix::Identity(dimension, dimension);
        const double tolerance = 1e-10 * std::max<Eigen::Index>(1, dimension);
        if ((completeness - identity).norm() > tolerance)
            throw std::invalid_argument("Kraus operators do not define a trace-preserving channel");
    }
 
    KrausOperators krausOperators_;
    size_t numQubits_;
    std::vector<double> pauliProbabilities_;
};
 
} // namespace Simulators
 
#endif // _SIMULATORS_QUANTUM_CHANNEL_H_