QCQI — Chapter 5: Density Operators & Quantum Operations

# Backend selector
backend = "backend"  # change: "pennylane", "braket", "pyquil", "qiskit"
print("Selected backend:", backend)
# Optional installs (if needed)
# !pip install cirq pennylane amazon-braket-sdk pyquil qiskit qiskit-aer

import numpy as np
I2 = np.eye(2, dtype=complex)
X = np.array([[0,1],[1,0]], dtype=complex)
Y = np.array([[0,-1j],[1j,0]], dtype=complex)
Z = np.array([[1,0],[0,-1]], dtype=complex)
def dm(psi): return psi @ psi.conj().T
def bloch(rho):
    return np.array([np.real(np.trace(rho@X)),
                     np.real(np.trace(rho@Y)),
                     np.real(np.trace(rho@Z))])
def trace_distance(rho, sigma):
    w = np.linalg.eigvalsh(rho - sigma)
    return 0.5*float(np.sum(np.abs(w)))
def fidelity(rho, sigma):
    A = rho @ (sigma @ rho)
    w = np.linalg.eigvalsh( (rho**0.5) @ sigma @ (rho**0.5) )  # fallback for qubits
    # for qubits you can also use explicit formula
    evals = np.linalg.eigvalsh((np.linalg.matrix_power(rho,1) @ sigma))
    # safer generic:
    from scipy.linalg import sqrtm
    M = sqrtm(rho) @ sigma @ sqrtm(rho)
    from numpy.linalg import eigvalsh
    vals = eigvalsh((M + M.conj().T)/2)
    vals = np.clip(vals, 0, None)
    return float(np.real(np.square(np.sum(np.sqrt(vals)))))

def apply_kraus(rho, Ks):
    # Check TP
    tp = sum([K.conj().T @ K for K in Ks])
    if not np.allclose(tp, I2, atol=1e-7):
        print("Warning: not TP to tolerance:", np.linalg.norm(tp - I2))
    out = sum([K @ rho @ K.conj().T for K in Ks])
    return out

def depolarizing(p):
    K0 = np.sqrt(1-p) * I2
    Kx = np.sqrt(p/3) * X
    Ky = np.sqrt(p/3) * Y
    Kz = np.sqrt(p/3) * Z
    return [K0, Kx, Ky, Kz]

def dephasing(p):
    K0 = np.sqrt(1-p) * I2
    K1 = np.sqrt(p) * Z
    return [K0/np.sqrt(2), K0/np.sqrt(2), K1/np.sqrt(2), -K1/np.sqrt(2)]  # equivalent expansion

def amplitude_damping(gamma):
    K0 = np.array([[1,0],[0,np.sqrt(1-gamma)]], dtype=complex)
    K1 = np.array([[0, np.sqrt(gamma)],[0,0]], dtype=complex)
    return [K0, K1]

def choi_from_kraus(Ks):
    # |Phi> = |00> + |11|  (unnormalized)
    Phi = np.zeros((4,4), dtype=complex)
    # Build |Phi><Phi|
    v = np.zeros((4,1), dtype=complex); v[0,0]=1; v[3,0]=1
    Phi = v @ v.conj().T
    # (E \otimes I) Phi
    J = sum([(np.kron(K, np.eye(2)) @ Phi @ np.kron(K.conj().T, np.eye(2))) for K in Ks])
    return J

def is_psd(M, tol=1e-9):
    w = np.linalg.eigvalsh((M + M.conj().T)/2)
    return bool(np.all(w > -tol)), w

def tp_check_from_choi(J):
    # trace over output (first subsystem if we took K acting on input); here dims are 2x2
    # reshape: indices (out,in; out',in') ordering if using kron(K, I)
    # We'll do explicit partial trace over "out" using basis
    Jmat = J.reshape(2,2,2,2)
    red = np.einsum('oi,oj->ij', Jmat[:, :, :, :].transpose(0,2,1,3).reshape(2,2,2,2), np.eye(2))
    # Simpler: compute partial trace via blocks
    J_blocks = [J[0:2,0:2], J[2:4,2:4]]
    return np.allclose(J_blocks[0]+J_blocks[1], np.eye(2), atol=1e-7)

Ks = amplitude_damping(0.3)
J = choi_from_kraus(Ks)
psd, eigs = is_psd(J)
print("Choi PSD?", psd, " min eig ~", float(np.min(eigs)))
# TP via sum K^\dagger K
print("TP check sum K^dag K ~ I:", np.allclose(sum([K.conj().T@K for K in Ks]), I2, atol=1e-7))

def bloch_affine(Ks):
    # Map rho = (I + r·σ)/2 -> (I + (A r + b) · σ)/2
    # Determine A (3x3) and b (3)
    B = []
    # Basis states to probe
    basis = [np.array([[1,0],[0,0]], complex), np.array([[0,0],[0,1]], complex),
             np.array([[0.5,0.5],[0.5,0.5]], complex), np.array([[0.5,-0.5j],[0.5j,0.5]], complex)]
    rs = [np.array([0,0,1]), np.array([0,0,-1]), np.array([1,0,0]), np.array([0,1,0])]
    outs = [apply_kraus(r, Ks) for r in basis]
    routs = [bloch(R) for R in outs]
    # Solve for A,b using linear system on these four probes
    # r_out = A r + b
    R = np.stack(rs, axis=1)  # 3x4
    Y = np.stack(routs, axis=1)  # 3x4
    # append row of ones to R for affine; solve Y = [A | b] @ [R; ones]
    ones = np.ones((1, R.shape[1]))
    M = np.vstack([R, ones])  # 4x4
    AB = Y @ np.linalg.inv(M)  # 3x4
    A = AB[:, :3]; b = AB[:, 3]
    return A, b
A,b = bloch_affine(amplitude_damping(0.3))
print("A ~\n", np.round(A,3)); print("b ~", np.round(b,3))

def rand_pure():
    v = np.random.default_rng(2).normal(size=(2,)) + 1j*np.random.default_rng(3).normal(size=(2,))
    v = v/np.linalg.norm(v)
    return v.reshape(2,1)
psi = np.array([[1],[0]], complex)
phi = (1/np.sqrt(2))*np.array([[1],[1]], complex)
rho = dm(psi); sigma = dm(phi)
D0 = trace_distance(rho, sigma)
# simple fidelity for qubits:
from numpy.linalg import eigvalsh
def Uhlmann(r, s):
    from scipy.linalg import sqrtm
    M = sqrtm(r) @ s @ sqrtm(r)
    ev = eigvalsh((M + M.conj().T)/2)
    ev = np.clip(ev,0,None)
    return float((np.sum(np.sqrt(ev)))**2)
F0 = Uhlmann(rho, sigma)
Ks = amplitude_damping(0.4)
r1 = apply_kraus(rho, Ks); s1 = apply_kraus(sigma, Ks)
D1 = trace_distance(r1, s1); F1 = Uhlmann(r1, s1)
print("D before, after =", round(D0,4), round(D1,4))
print("F before, after =", round(F0,4), round(F1,4))

# POVM: E0 = |0><0| * 0.8 + |1><1| * 0.2  (diagonal's fine for demo), E1=I-E0
E0 = np.array([[0.8,0],[0,0.2]], complex); E1 = I2 - E0
# Kraus via spectral decomposition
w0, V0 = np.linalg.eigh(E0)
K0 = V0 @ np.diag(np.sqrt(np.clip(w0,0,None))) @ V0.conj().T
w1, V1 = np.linalg.eigh(E1)
K1 = V1 @ np.diag(np.sqrt(np.clip(w1,0,None))) @ V1.conj().T
# Isometry V = K0 ⊗ |0> + K1 ⊗ |1>  (map S -> S⊗E)
V = np.block([[K0, K1]])
# Check probabilities on |+>
plus = (1/np.sqrt(2))*np.array([[1],[1]], complex)
rho = dm(plus)
p0 = np.real(np.trace(E0 @ rho)); p1 = np.real(np.trace(E1 @ rho))
print("POVM probs (direct) :", round(p0,4), round(p1,4))
# Direct Kraus as channel with classical register omitted
rho0 = K0 @ rho @ K0.conj().T; rho1 = K1 @ rho @ K1.conj().T
print("Sum prob via Kraus:", round(np.trace(rho0).real,4), round(np.trace(rho1).real,4))

if backend == "cirq":
    import cirq
    q = cirq.LineQubit(0)
    gamma = 0.4
    c = cirq.Circuit(cirq.H(q),
                     cirq.amplitude_damp(gamma).on(q),
                     cirq.measure(q, key='m'))
    print(c)
    sim = cirq.DensityMatrixSimulator()
    res = sim.run(c, repetitions=2000)
    print("Counts:", dict(res.histogram(key='m')))

if backend == "pennylane":
    import pennylane as qml, numpy as np
    dev = qml.device("default.mixed", wires=1, shots=2000)
    p = 0.3
    @qml.qnode(dev)
    def circ():
        qml.Hadamard(0)
        qml.PhaseDamping(p, wires=0)
        return qml.sample(qml.PauliZ(0))
    s = circ()
    import numpy as np
    zeros = int((s==0).sum()); ones = int((s==1).sum())
    print({"0": zeros, "1": ones})

if backend == "braket":
    from braket.circuits import Circuit, noises
    from braket.devices import LocalSimulator
    circ = Circuit().h(0).apply_gate_noise(noises.AmplitudeDamping(0.4), target_qubits=[0]).measure(0)
    dev = LocalSimulator("braket_dm")
    res = dev.run(circ, shots=2000).result().measurement_counts
    print(res)

if backend == "pyquil":
    from pyquil import Program
    from pyquil.gates import H, MEASURE
    from pyquil.api import get_qc
    # QVM doesn't natively support Kraus, but we can sample H then post-process probabilities
    p = Program(); ro = p.declare('ro','BIT',1); p += H(0); p += MEASURE(0, ro[0])
    qc = get_qc('1q-qvm'); res = qc.run(p.wrap_in_numshots_loop(2000))
    bits = res.get_register_map()['ro']
    print("Raw counts ~", {0:int((bits==0).sum()), 1:int((bits==1).sum())})

if backend == "qiskit":
    from qiskit import QuantumCircuit
    from qiskit_aer import AerSimulator
    from qiskit_aer.noise import amplitude_damping_error
    qc = QuantumCircuit(1,1); qc.h(0); qc.measure(0,0)
    noise_model = None
    try:
        from qiskit_aer.noise import NoiseModel
        noise_model = NoiseModel()
        amp = amplitude_damping_error(0.4)
        noise_model.add_quantum_error(amp, ['measure'], [0])
    except Exception as e:
        pass
    sim = AerSimulator(method='density_matrix')
    res = sim.run(qc, noise_model=noise_model, shots=2000).result().get_counts()
    print(res)