Optimal POVM¶

Adapted from QGopt/docs/optimal-povm

In [1]:

import numpy as np
import torch

try:
    import numqi
except ImportError:
    %pip install numqi
    import numqi

import numpy as np import torch try: import numqi except ImportError: %pip install numqi import numqi

In [2]:

class OptimalPOVMModel(torch.nn.Module):
    def __init__(self):
        super().__init__()
        self.manifold = numqi.manifold.QuantumChannel(dim_in=2, dim_out=2, choi_rank=4, dtype=torch.complex128, return_kind='kraus')
        self.rho_list_T = []

    def set_rho_list(self, rho_list):
        assert (rho_list.ndim==3) and (rho_list.shape[1]==2) and (rho_list.shape[2]==2)
        self.rho_list_T = torch.tensor(rho_list.transpose(0,2,1).reshape(-1,4), dtype=torch.complex128)

    def get_povm(self, return_numpy=True):
        tmp1 = self.manifold()
        ret = torch.einsum(tmp1, [0,1,2], tmp1.conj(), [0,1,3], [0,2,3])
        if return_numpy:
            ret = ret.detach().numpy()
        return ret

    def forward(self):
        povm = self.get_povm(return_numpy=False)
        tmp0 = torch.linalg.inv(povm.reshape(4,4))
        T = torch.stack([tmp0[0].real, tmp0[3].real, tmp0[2].real, tmp0[2].imag])
        prob = torch.einsum(self.rho_list_T, [0,1], povm.reshape(4,4), [2,1], [0,2]).real
        tmp0 = torch.einsum(prob, [0,1], prob, [0,2], [0,1,2])
        tmp1 = (torch.diag_embed(prob, dim1=1, dim2=2) - tmp0).mean(axis=0)
        covariance = torch.einsum(T, [0,1], tmp1, [1,2], T, [3,2], [0,3])
        EVL = torch.linalg.eigvalsh(covariance)
        loss = torch.log(EVL[1:]).sum()
        return loss

class OptimalPOVMModel(torch.nn.Module): def __init__(self): super().__init__() self.manifold = numqi.manifold.QuantumChannel(dim_in=2, dim_out=2, choi_rank=4, dtype=torch.complex128, return_kind='kraus') self.rho_list_T = [] def set_rho_list(self, rho_list): assert (rho_list.ndim==3) and (rho_list.shape[1]==2) and (rho_list.shape[2]==2) self.rho_list_T = torch.tensor(rho_list.transpose(0,2,1).reshape(-1,4), dtype=torch.complex128) def get_povm(self, return_numpy=True): tmp1 = self.manifold() ret = torch.einsum(tmp1, [0,1,2], tmp1.conj(), [0,1,3], [0,2,3]) if return_numpy: ret = ret.detach().numpy() return ret def forward(self): povm = self.get_povm(return_numpy=False) tmp0 = torch.linalg.inv(povm.reshape(4,4)) T = torch.stack([tmp0[0].real, tmp0[3].real, tmp0[2].real, tmp0[2].imag]) prob = torch.einsum(self.rho_list_T, [0,1], povm.reshape(4,4), [2,1], [0,2]).real tmp0 = torch.einsum(prob, [0,1], prob, [0,2], [0,1,2]) tmp1 = (torch.diag_embed(prob, dim1=1, dim2=2) - tmp0).mean(axis=0) covariance = torch.einsum(T, [0,1], tmp1, [1,2], T, [3,2], [0,3]) EVL = torch.linalg.eigvalsh(covariance) loss = torch.log(EVL[1:]).sum() return loss

Generate random density matrix from Bloch sphere.

In [3]:

num_sample = 10000
tmp0 = numqi.random.rand_n_ball(3, size=num_sample).T.reshape(3,-1,1,1)/2
rho_list = numqi.gate.I/2 + numqi.gate.X*tmp0[0] + numqi.gate.Y*tmp0[1] + numqi.gate.Z*tmp0[2]

num_sample = 10000 tmp0 = numqi.random.rand_n_ball(3, size=num_sample).T.reshape(3,-1,1,1)/2 rho_list = numqi.gate.I/2 + numqi.gate.X*tmp0[0] + numqi.gate.Y*tmp0[1] + numqi.gate.Z*tmp0[2]

Build model for optimal POVM.

In [4]:

model = OptimalPOVMModel()
model.set_rho_list(rho_list)
theta_optim = numqi.optimize.minimize(model, theta0='uniform', num_repeat=1, tol=1e-7)

model = OptimalPOVMModel() model.set_rho_list(rho_list) theta_optim = numqi.optimize.minimize(model, theta0='uniform', num_repeat=1, tol=1e-7)

[round=0] min(f)=-0.3765838821647226, current(f)=-0.3765838821647226

Theoretically, the optimal POVM gives

$$ \mathrm{Tr}[M_iM_j]=\frac{2\delta_{ij}+1}{12} $$

In [5]:

povm = model.get_povm()
cross = np.einsum(povm, [0,1,2], povm, [3,2,1], [0,3], optimize=True).real
print(f'Founded POVM: Tr[Mi Mj]=\n{np.around(cross, 5)}\n')

cross_analytical = 1/12 + 1/6 * np.eye(4)
print(f'Optimal POVM (theoretical): Tr[Mi Mj]=\n{np.around(cross_analytical, 5)}')

povm = model.get_povm() cross = np.einsum(povm, [0,1,2], povm, [3,2,1], [0,3], optimize=True).real print(f'Founded POVM: Tr[Mi Mj]=\n{np.around(cross, 5)}\n') cross_analytical = 1/12 + 1/6 * np.eye(4) print(f'Optimal POVM (theoretical): Tr[Mi Mj]=\n{np.around(cross_analytical, 5)}')

Founded POVM: Tr[Mi Mj]=
[[0.25297 0.08342 0.08333 0.08324]
 [0.08342 0.24726 0.08324 0.08333]
 [0.08333 0.08324 0.25073 0.08342]
 [0.08324 0.08333 0.08342 0.24906]]

Optimal POVM (theoretical): Tr[Mi Mj]=
[[0.25    0.08333 0.08333 0.08333]
 [0.08333 0.25    0.08333 0.08333]
 [0.08333 0.08333 0.25    0.08333]
 [0.08333 0.08333 0.08333 0.25   ]]