QCQI — Chapter 11: Quantum Shannon Theory & Channel Capacities

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

import numpy as np
from numpy.linalg import eigvalsh
def dm(psi): psi=np.asarray(psi).reshape(-1,1); return psi@psi.conj().T
def entropy(rho, tol=1e-12):
    vals = np.clip(eigvalsh((rho+rho.conj().T)/2), 0, 1)
    vals = vals[vals>tol]
    return float(-np.sum(vals*np.log2(vals)))
def kron(*Ms):
    out = np.array([[1.0+0j]])
    for M in Ms:
        out = np.kron(out, M)
    return out
def partial_trace(rho, dims, keep):
    rho = rho.reshape(*(dims+dims))
    axes = tuple(i for i in range(len(dims)) if i not in keep)
    for ax in sorted(axes, reverse=True):
        rho = np.trace(rho, axis1=ax, axis2=ax+len(dims))
    dkeep = int(np.prod([dims[i] for i in keep]))
    return rho.reshape(dkeep, dkeep)

I2 = np.eye(2, dtype=complex); X = np.array([[0,1],[1,0]], complex)
Y = np.array([[0,-1j],[1j,0]], complex); Z = np.array([[1,0],[0,-1]], complex)
def depolarizing_channel(p):
    K = [np.sqrt(1-p)*I2, np.sqrt(p/3)*X, np.sqrt(p/3)*Y, np.sqrt(p/3)*Z]
    return K
def amplitude_damping(gamma):
    E0 = np.array([[1,0],[0,np.sqrt(1-gamma)]], complex)
    E1 = np.array([[0,np.sqrt(gamma)],[0,0]], complex)
    return [E0,E1]
def apply_kraus(rho, Ks):
    out = np.zeros_like(rho, dtype=complex)
    for K in Ks:
        out += K @ rho @ K.conj().T
    return out

def holevo_of_ensemble(ensemble, Ks):
    # ensemble = list of (p, rho) with sum p = 1
    sigmas = []
    avg = np.zeros((2,2), complex)
    Savg = 0.0; Savgs = 0.0
    for p, rho in ensemble:
        sigma = apply_kraus(rho, Ks)
        sigmas.append((p, sigma))
        avg += p * sigma
    chi = entropy(avg) - sum(p*entropy(s) for p,s in sigmas)
    return float(chi)
# Example: depolarizing channel with BB84-like states
H = (1/np.sqrt(2))*np.array([[1,1],[1,-1]], complex)
states = [dm(np.array([1,0])), dm(np.array([0,1])), dm(H@np.array([1,0])), dm(H@np.array([0,1]))]
ensemble = [(0.25, s) for s in states]
for p in [0.0, 0.05, 0.2, 0.4]:
    chi = holevo_of_ensemble(ensemble, depolarizing_channel(p))
    print(f"Depolarizing p={p:.2f}: Holevo chi ≈ {chi:.3f} bits")

def coherent_information(rho, Ks):
    # Purify rho on R with |Φ> = sum_i sqrt(λ_i)|i>_R|i>_A in eigenbasis of rho
    vals, vecs = np.linalg.eigh(rho)
    vals = np.clip(vals, 0, 1)
    # Build purification in RA (dimension 2x2)
    psi = sum(np.sqrt(vals[i]) * np.kron(vecs[:,i], vecs[:,i]) for i in range(2))
    psi = psi.reshape(4,1); rhoRA = psi @ psi.conj().T
    # Apply channel on A (second qubit)
    # Kraus on A: (I⊗K) ρ (I⊗K†)
    out = np.zeros((4,4), complex)
    for K in Ks:
        Kfull = np.kron(I2, K)
        out += Kfull @ rhoRA @ Kfull.conj().T
    rhoB = partial_trace(out, [2,2], keep=(1,))
    return entropy(rhoB) - entropy(out)
# Scan over z-basis pure states for amplitude damping
def rho_on_bloch(theta):
    # pure state cos(theta/2)|0> + sin(theta/2)|1>
    v = np.array([np.cos(theta/2), np.sin(theta/2)], complex)
    return dm(v)
for gamma in [0.0, 0.2, 0.4]:
    Ks = amplitude_damping(gamma)
    thetas = np.linspace(0, np.pi, 41)
    vals = [coherent_information(rho_on_bloch(t), Ks) for t in thetas]
    print(f"Amplitude damping γ={gamma:.1f}: max I_c ≈ {max(vals):.3f} qubits")

def entanglement_assisted_MI(rho, Ks):
    # Build purification |Φ_rho> as before
    vals, vecs = np.linalg.eigh(rho)
    psi = sum(np.sqrt(max(vals[i],0)) * np.kron(vecs[:,i], vecs[:,i]) for i in range(2))
    psi = psi.reshape(4,1); rhoRA = psi @ psi.conj().T
    out = np.zeros((4,4), complex)
    for K in Ks:
        out += np.kron(I2, K) @ rhoRA @ np.kron(I2, K.conj().T)
    rhoR = partial_trace(out, [2,2], keep=(0,))
    rhoB = partial_trace(out, [2,2], keep=(1,))
    return entropy(rhoR) + entropy(rhoB) - entropy(out)
# For depolarizing, scan over Bloch radius (should peak near maximally mixed)
def rho_bloch_cart(r, nx, ny, nz):
    v = np.array([1 + 0j, 0, 0, 1], complex).reshape(2,2)
    n = np.array([nx, ny, nz]); n = n / (np.linalg.norm(n)+1e-12)
    rho = 0.5*(I2 + r*(n[0]*X + n[1]*Y + n[2]*Z))
    return rho
for p in [0.0, 0.1, 0.3]:
    Ks = depolarizing_channel(p)
    rs = np.linspace(0,1,11)
    vals = [entanglement_assisted_MI(rho_bloch_cart(r,0,0,1), Ks) for r in rs]
    print(f"Depolarizing p={p:.1f}: max I(R:B) ≈ {max(vals):.3f} bits (coarse scan)")

# We'll build the typical projector for a qubit source diagonal in Z with probs (q,1-q).
def typical_projector_Z(q, n=3):
    # keep sequences with Hamming weight close to n*(1-q)
    import itertools, math
    H = lambda x: -x*math.log2(x)-(1-x)*math.log2(1-x) if 0<x<1 else 0.0
    target = (1-q)*n
    keep = []
    for s in itertools.product([0,1], repeat=n):
        w = sum(s)
        if abs(w - target) <= 1:  # crude "typical" band
            keep.append(s)
    # Build projector on the span of |s>
    dim = 2**n
    P = np.zeros((dim,dim), complex)
    for s in keep:
        idx = sum(int(b)*(2**i) for i,b in enumerate(reversed(s)))
        e = np.zeros((dim,1), complex); e[idx,0]=1
        P += e@e.conj().T
    return P
P = typical_projector_Z(0.7, n=3); print("Typical subspace rank (n=3, q=0.7):", int(np.round(np.trace(P).real)))

if backend == "cirq":
    import cirq
    # Simple encoder for a 2-symbol source into a 1-qubit register by isometry (toy)
    q = cirq.LineQubit.range(2)
    c = cirq.Circuit()
    # Prepare |ψ0>=|0>, |ψ1>=cos α|0>+sin α|1> with prob 1/2 each; here just build the isometry stub
    alpha = 0.6
    c.append(cirq.ry(alpha).on(q[0]))  # pretend "typical" rotate
    c.append(cirq.CNOT(q[0], q[1]))
    print("Cirq stub circuit (encoder skeleton) built.")