maximum entropy method¶

see paper "Maximum entropy methods for quantum state compatibility problems" arxiv-link for more details.

In [1]:

import functools
import numpy as np
import matplotlib.pyplot as plt

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

hf_kron = lambda *x: functools.reduce(np.kron, x)
hf_trace0 = lambda x: x-(np.trace(x)/x.shape[0])*np.eye(x.shape[0])
np_rng = np.random.default_rng()

import functools import numpy as np import matplotlib.pyplot as plt try: import numqi except ImportError: %pip install numqi import numqi hf_kron = lambda *x: functools.reduce(np.kron, x) hf_trace0 = lambda x: x-(np.trace(x)/x.shape[0])*np.eye(x.shape[0]) np_rng = np.random.default_rng()

2 qubits, Pauli Operator¶

In [2]:

op_list = [
    hf_kron(numqi.gate.X, numqi.gate.X),
    hf_kron(numqi.gate.Z, numqi.gate.I),
]
model = numqi.maximum_entropy.MaximumEntropyModel(op_list)

term_value = np.array([0.5, 0.5])
model.set_target(term_value)
theta_optim0 = numqi.optimize.minimize(model, theta0='uniform', num_repeat=1, tol=1e-10)
model.theta
model.dm_torch
model.term_value

term_value = np.array([2, 2])
model.set_target(term_value)
theta_optim0 = numqi.optimize.minimize(model, theta0='uniform', num_repeat=1, tol=1e-10)

term_value_target, term_value_list, EVL_list = numqi.maximum_entropy.get_maximum_entropy_model_boundary(model, radius=1.5)
term_value = np.array([-1.3,-1.5])
coeffA, coeffC = model.get_witness(term_value)
fig,ax = numqi.maximum_entropy.draw_maximum_entropy_model_boundary(term_value_target, term_value_list,
                EVL_list, witnessA=coeffA, witnessC=coeffC)

op_list = [ hf_kron(numqi.gate.X, numqi.gate.X), hf_kron(numqi.gate.Z, numqi.gate.I), ] model = numqi.maximum_entropy.MaximumEntropyModel(op_list) term_value = np.array([0.5, 0.5]) model.set_target(term_value) theta_optim0 = numqi.optimize.minimize(model, theta0='uniform', num_repeat=1, tol=1e-10) model.theta model.dm_torch model.term_value term_value = np.array([2, 2]) model.set_target(term_value) theta_optim0 = numqi.optimize.minimize(model, theta0='uniform', num_repeat=1, tol=1e-10) term_value_target, term_value_list, EVL_list = numqi.maximum_entropy.get_maximum_entropy_model_boundary(model, radius=1.5) term_value = np.array([-1.3,-1.5]) coeffA, coeffC = model.get_witness(term_value) fig,ax = numqi.maximum_entropy.draw_maximum_entropy_model_boundary(term_value_target, term_value_list, EVL_list, witnessA=coeffA, witnessC=coeffC)

[round=0] min(f)=3.518634948575557e-17, current(f)=3.518634948575557e-17
[round=0] min(f)=3.343145758020391, current(f)=3.343145758020391

[round=0] min(f)=0.970113352635388, current(f)=0.970113352635388

3 qubits, random local operators¶

In [3]:

np_rng = np.random.default_rng(233) #233
tmp0 = [
    hf_kron(numqi.random.rand_hermitian_matrix(4, seed=np_rng), numqi.gate.I),
    hf_kron(numqi.gate.I, numqi.random.rand_hermitian_matrix(4, seed=np_rng)),
]
tmp0 = [hf_trace0(x) for x in tmp0]
op_list = [x*(2/np.linalg.norm(x.reshape(-1))) for x in tmp0] #make it norm-2, better in plotting
model = numqi.maximum_entropy.MaximumEntropyModel(op_list)

term_value_target, term_value_list, EVL_list = numqi.maximum_entropy.get_maximum_entropy_model_boundary(model, radius=1.5, index=(0,1))
term_value = np.array([-1.4,-0.5])
coeffA, coeffC = model.get_witness(term_value)
assert coeffA is not None
fig,ax = numqi.maximum_entropy.draw_maximum_entropy_model_boundary(term_value_target, term_value_list,
                EVL_list, witnessA=coeffA, witnessC=coeffC)

np_rng = np.random.default_rng(233) #233 tmp0 = [ hf_kron(numqi.random.rand_hermitian_matrix(4, seed=np_rng), numqi.gate.I), hf_kron(numqi.gate.I, numqi.random.rand_hermitian_matrix(4, seed=np_rng)), ] tmp0 = [hf_trace0(x) for x in tmp0] op_list = [x*(2/np.linalg.norm(x.reshape(-1))) for x in tmp0] #make it norm-2, better in plotting model = numqi.maximum_entropy.MaximumEntropyModel(op_list) term_value_target, term_value_list, EVL_list = numqi.maximum_entropy.get_maximum_entropy_model_boundary(model, radius=1.5, index=(0,1)) term_value = np.array([-1.4,-0.5]) coeffA, coeffC = model.get_witness(term_value) assert coeffA is not None fig,ax = numqi.maximum_entropy.draw_maximum_entropy_model_boundary(term_value_target, term_value_list, EVL_list, witnessA=coeffA, witnessC=coeffC)

[round=0] min(f)=0.13754543304124442, current(f)=0.13754543304124442

Maximum entropy method can produce witness.

In [4]:

num_qubit = 3
op_list = numqi.maximum_entropy.get_1dchain_2local_pauli_basis(num_qubit, with_I=False)
model = numqi.maximum_entropy.MaximumEntropyModel(op_list)

state = np.array([0,1,1,0,1,0,0,0])/np.sqrt(3) #W-state
term_value_target = ((op_list @ state) @ state.conj()).real
model.set_target(term_value_target)
theta_optim0 = numqi.optimize.minimize(model, theta0='uniform', num_repeat=3, tol=1e-12)
if theta_optim0.fun < 1e-9: #fail to converge sometime, just re-run it
    rho = model.dm_torch.detach().numpy()
    rank = np.sum(np.linalg.eigvalsh(rho)>1e-4)
    assert rank==1

state = np.array([1,0,0,0,0,0,0,1])/np.sqrt(2) #GHZ-state
term_value_target = ((op_list @ state) @ state.conj()).real
model.set_target(term_value_target)
theta_optim0 = numqi.optimize.minimize(model, theta0='uniform', num_repeat=3, tol=1e-12)
if theta_optim0.fun < 1e-9: #fail to converge sometime, just re-run it
    rho = model.dm_torch.detach().numpy()
    rank = np.sum(np.linalg.eigvalsh(rho)>1e-4)
    assert rank==2

    coeffA,coeffC = model.get_witness(term_value_target*1.1)
    assert coeffA is not None
    for _ in range(1000):
        rho = numqi.random.rand_density_matrix(2**num_qubit)
        z0 = np.trace(op_list @ rho, axis1=1, axis2=2).real
        assert np.dot(z0 - coeffA, coeffC) >= 0

num_qubit = 3 op_list = numqi.maximum_entropy.get_1dchain_2local_pauli_basis(num_qubit, with_I=False) model = numqi.maximum_entropy.MaximumEntropyModel(op_list) state = np.array([0,1,1,0,1,0,0,0])/np.sqrt(3) #W-state term_value_target = ((op_list @ state) @ state.conj()).real model.set_target(term_value_target) theta_optim0 = numqi.optimize.minimize(model, theta0='uniform', num_repeat=3, tol=1e-12) if theta_optim0.fun < 1e-9: #fail to converge sometime, just re-run it rho = model.dm_torch.detach().numpy() rank = np.sum(np.linalg.eigvalsh(rho)>1e-4) assert rank==1 state = np.array([1,0,0,0,0,0,0,1])/np.sqrt(2) #GHZ-state term_value_target = ((op_list @ state) @ state.conj()).real model.set_target(term_value_target) theta_optim0 = numqi.optimize.minimize(model, theta0='uniform', num_repeat=3, tol=1e-12) if theta_optim0.fun < 1e-9: #fail to converge sometime, just re-run it rho = model.dm_torch.detach().numpy() rank = np.sum(np.linalg.eigvalsh(rho)>1e-4) assert rank==2 coeffA,coeffC = model.get_witness(term_value_target*1.1) assert coeffA is not None for _ in range(1000): rho = numqi.random.rand_density_matrix(2**num_qubit) z0 = np.trace(op_list @ rho, axis1=1, axis2=2).real assert np.dot(z0 - coeffA, coeffC) >= 0

[round=0] min(f)=1.3311598878066474e-10, current(f)=1.3311598878066474e-10

[round=1] min(f)=1.3311598878066474e-10, current(f)=4.612732138078387e-08

[round=2] min(f)=1.3311598878066474e-10, current(f)=1.4847696843901684e-08
[round=0] min(f)=8.13419532911629e-12, current(f)=8.13419532911629e-12

[round=1] min(f)=8.13419532911629e-12, current(f)=2.454520574161648e-11
[round=2] min(f)=3.5732396639786645e-12, current(f)=3.5732396639786645e-12

[round=0] min(f)=0.020000012798044944, current(f)=0.020000012798044944

4 qubits, entanglement monogamy¶

TODO more details

In [5]:

num_qubit = 4
op_2qubit_list = numqi.maximum_entropy.get_1dchain_2local_pauli_basis(2, with_I=False)
op_list = numqi.maximum_entropy.get_1dchain_2local_pauli_basis(num_qubit, with_I=False)
model = numqi.maximum_entropy.MaximumEntropyModel(op_list)

state = np.array([1,0,0,1])/np.sqrt(2) #Bell state
tmp0 = ((op_2qubit_list @ state) @ state.conj()).real
term_value_target = np.tile(tmp0, num_qubit-1)

# due to entanglement monogamy, this must be a witness
coeffA,coeffC = model.get_witness(term_value_target*1.1)
assert coeffA is not None
for _ in range(1000):
    rho = numqi.random.rand_density_matrix(2**num_qubit)
    z0 = np.trace(op_list @ rho, axis1=1, axis2=2).real
    assert np.dot(z0 - coeffA, coeffC) >= 0

num_qubit = 4 op_2qubit_list = numqi.maximum_entropy.get_1dchain_2local_pauli_basis(2, with_I=False) op_list = numqi.maximum_entropy.get_1dchain_2local_pauli_basis(num_qubit, with_I=False) model = numqi.maximum_entropy.MaximumEntropyModel(op_list) state = np.array([1,0,0,1])/np.sqrt(2) #Bell state tmp0 = ((op_2qubit_list @ state) @ state.conj()).real term_value_target = np.tile(tmp0, num_qubit-1) # due to entanglement monogamy, this must be a witness coeffA,coeffC = model.get_witness(term_value_target*1.1) assert coeffA is not None for _ in range(1000): rho = numqi.random.rand_density_matrix(2**num_qubit) z0 = np.trace(op_list @ rho, axis1=1, axis2=2).real assert np.dot(z0 - coeffA, coeffC) >= 0

[round=0] min(f)=1.5563379819561582, current(f)=1.5563379819561582

[round=0] min(f)=4.175245018373061e-09, current(f)=4.175245018373061e-09

4 qubits, random local operators¶

3 minutes on my laptop, hard to converge (the plotting is not good)

In [6]:

# np_rng = np.random.default_rng(233)
# num_qubit = 4
# op_2qubit_list = numqi.maximum_entropy.get_1dchain_2local_pauli_basis(2, with_I=False)
# op_list = numqi.maximum_entropy.get_1dchain_2local_pauli_basis(num_qubit, with_I=False)
# model = numqi.maximum_entropy.MaximumEntropyModel(op_list)

# tmp0 = [numqi.random.rand_density_matrix(4, seed=np_rng) for _ in range(num_qubit-1)]
# term_value = np.concatenate([np.trace(op_2qubit_list@x, axis1=1, axis2=2).real for x in tmp0], axis=0)

# index = np.sort(np_rng.permutation(len(term_value))[:2])
# term_value_target, term_value_list, EVL_list = numqi.maximum_entropy.get_maximum_entropy_model_boundary(
#         model, radius=1.2, index=index, term_value_target=term_value, tol=1e-7, num_repeat=3, num_point=50)
# fig,ax = numqi.maximum_entropy.draw_maximum_entropy_model_boundary(term_value_target, term_value_list, EVL_list, index=index, rank_radius=-1)

# coeffA,coeffC = model.get_witness(term_value_target[20])
# if coeffA is not None:
#     fig,ax = numqi.maximum_entropy.draw_maximum_entropy_model_boundary(term_value_target, term_value_list, EVL_list,
#                 index=index, witnessA=coeffA, witnessC=coeffC, rank_radius=-1)
#     for _ in range(1000):
#         rho = numqi.random.rand_density_matrix(2**num_qubit)
#         z0 = np.trace(op_list @ rho, axis1=1, axis2=2).real
#         assert np.dot(z0 - coeffA, coeffC) >= 0

# np_rng = np.random.default_rng(233) # num_qubit = 4 # op_2qubit_list = numqi.maximum_entropy.get_1dchain_2local_pauli_basis(2, with_I=False) # op_list = numqi.maximum_entropy.get_1dchain_2local_pauli_basis(num_qubit, with_I=False) # model = numqi.maximum_entropy.MaximumEntropyModel(op_list) # tmp0 = [numqi.random.rand_density_matrix(4, seed=np_rng) for _ in range(num_qubit-1)] # term_value = np.concatenate([np.trace(op_2qubit_list@x, axis1=1, axis2=2).real for x in tmp0], axis=0) # index = np.sort(np_rng.permutation(len(term_value))[:2]) # term_value_target, term_value_list, EVL_list = numqi.maximum_entropy.get_maximum_entropy_model_boundary( # model, radius=1.2, index=index, term_value_target=term_value, tol=1e-7, num_repeat=3, num_point=50) # fig,ax = numqi.maximum_entropy.draw_maximum_entropy_model_boundary(term_value_target, term_value_list, EVL_list, index=index, rank_radius=-1) # coeffA,coeffC = model.get_witness(term_value_target[20]) # if coeffA is not None: # fig,ax = numqi.maximum_entropy.draw_maximum_entropy_model_boundary(term_value_target, term_value_list, EVL_list, # index=index, witnessA=coeffA, witnessC=coeffC, rank_radius=-1) # for _ in range(1000): # rho = numqi.random.rand_density_matrix(2**num_qubit) # z0 = np.trace(op_list @ rho, axis1=1, axis2=2).real # assert np.dot(z0 - coeffA, coeffC) >= 0