In [1]:
Copied!
import numpy as np
import matplotlib.pyplot as plt
from tqdm.auto import tqdm
try:
import numqi
except ImportError:
%pip install numqi
import numqi
np_rng = np.random.default_rng()
cp_tableau = ['#006BA4', '#FF800E', '#ABABAB', '#595959', '#5F9ED1', '#C85200', '#898989', '#A2C8EC', '#FFBC79', '#CFCFCF']
import numpy as np
import matplotlib.pyplot as plt
from tqdm.auto import tqdm
try:
import numqi
except ImportError:
%pip install numqi
import numqi
np_rng = np.random.default_rng()
cp_tableau = ['#006BA4', '#FF800E', '#ABABAB', '#595959', '#5F9ED1', '#C85200', '#898989', '#A2C8EC', '#FFBC79', '#CFCFCF']
3PB, 4PB, 5PB¶
In [2]:
Copied!
dim_list = list(range(3, 9))
kwargs = dict(num_repeat=5, early_stop_threshold=1e-10, converge_tol=1e-12, dtype='float64')
# ~3 min (num_repeat=5)
udp_3pb_loss_list = []
uda_4pb_loss_list = []
uda_5pb_loss_list = []
udp_4pb_loss_list = []
udp_5pb_loss_list = []
for dim in dim_list:
print(f'd={dim}')
alpha = np.pi/dim
matB = numqi.unique_determine.get_chebshev_orthonormal(dim, alpha, with_computational_basis=False)[:(-dim)]
# matB = numqi.matrix_space.get_matrix_orthogonal_basis(matB, field='real')[0]
udp_3pb_loss_list.append(numqi.unique_determine.check_UD('udp', matB, **kwargs)[1])
matB = numqi.unique_determine.get_chebshev_orthonormal(dim, alpha, with_computational_basis=False)
udp_4pb_loss_list.append(numqi.unique_determine.check_UD('udp', matB, **kwargs)[1])
uda_4pb_loss_list.append(numqi.unique_determine.check_UD('uda', matB, **kwargs)[1])
matB = numqi.unique_determine.get_chebshev_orthonormal(dim, alpha, with_computational_basis=True)
udp_5pb_loss_list.append(numqi.unique_determine.check_UD('udp', matB, **kwargs)[1])
uda_5pb_loss_list.append(numqi.unique_determine.check_UD('uda', matB, **kwargs)[1])
dim_list = list(range(3, 9))
kwargs = dict(num_repeat=5, early_stop_threshold=1e-10, converge_tol=1e-12, dtype='float64')
# ~3 min (num_repeat=5)
udp_3pb_loss_list = []
uda_4pb_loss_list = []
uda_5pb_loss_list = []
udp_4pb_loss_list = []
udp_5pb_loss_list = []
for dim in dim_list:
print(f'd={dim}')
alpha = np.pi/dim
matB = numqi.unique_determine.get_chebshev_orthonormal(dim, alpha, with_computational_basis=False)[:(-dim)]
# matB = numqi.matrix_space.get_matrix_orthogonal_basis(matB, field='real')[0]
udp_3pb_loss_list.append(numqi.unique_determine.check_UD('udp', matB, **kwargs)[1])
matB = numqi.unique_determine.get_chebshev_orthonormal(dim, alpha, with_computational_basis=False)
udp_4pb_loss_list.append(numqi.unique_determine.check_UD('udp', matB, **kwargs)[1])
uda_4pb_loss_list.append(numqi.unique_determine.check_UD('uda', matB, **kwargs)[1])
matB = numqi.unique_determine.get_chebshev_orthonormal(dim, alpha, with_computational_basis=True)
udp_5pb_loss_list.append(numqi.unique_determine.check_UD('udp', matB, **kwargs)[1])
uda_5pb_loss_list.append(numqi.unique_determine.check_UD('uda', matB, **kwargs)[1])
d=3
d=4
d=5
d=6
d=7
d=8
In [3]:
Copied!
fig,(ax0,ax1) = plt.subplots(1,2,figsize=(8,4))
ax0.plot(dim_list, udp_4pb_loss_list, label='4PBs', marker='.')
ax0.plot(dim_list, udp_3pb_loss_list, label='3PBs', marker='x')
ax1.plot(dim_list, uda_5pb_loss_list, label='5PBs', marker='.')
ax1.plot(dim_list, uda_4pb_loss_list, label='4PBs', marker='x')
for ax in [ax0,ax1]:
ax.set_xlabel('qudit $d$')
ax.legend()
ax.set_yscale('log')
ax.grid()
ax.set_ylim(3e-14, 0.4)
ax0.set_ylabel(r'UDP loss')
ax1.set_ylabel(r'UDA loss')
ax1.yaxis.tick_right()
ax1.yaxis.set_label_position('right')
fig.tight_layout()
fig,(ax0,ax1) = plt.subplots(1,2,figsize=(8,4))
ax0.plot(dim_list, udp_4pb_loss_list, label='4PBs', marker='.')
ax0.plot(dim_list, udp_3pb_loss_list, label='3PBs', marker='x')
ax1.plot(dim_list, uda_5pb_loss_list, label='5PBs', marker='.')
ax1.plot(dim_list, uda_4pb_loss_list, label='4PBs', marker='x')
for ax in [ax0,ax1]:
ax.set_xlabel('qudit $d$')
ax.legend()
ax.set_yscale('log')
ax.grid()
ax.set_ylim(3e-14, 0.4)
ax0.set_ylabel(r'UDP loss')
ax1.set_ylabel(r'UDA loss')
ax1.yaxis.tick_right()
ax1.yaxis.set_label_position('right')
fig.tight_layout()
Stability of 5PB (UDA)¶
In [4]:
Copied!
# about 3 mins
# to save time, we choose a relative small value for 'num_repeat' and 'num_random'
dim_qudit = 6
num_random = 50
cvxpy_eps = 1e-6
noise_rate_list = np.logspace(-6, -3, 10)
matrix_subspace = numqi.unique_determine.get_chebshev_orthonormal(dim_qudit, alpha=np.pi/dim_qudit, with_computational_basis=True)
tag_ud,loss,model = numqi.unique_determine.check_UD('uda', matrix_subspace, num_repeat=5,
converge_tol=1e-10, early_stop_threshold=1e-10, dtype='float64', tag_print=2, return_model=True)
matH = model.matH.numpy().copy()
EVL,EVC = np.linalg.eigh(matH)
state_list = [EVC[:,0]] + [numqi.random.rand_haar_state(dim_qudit) for _ in range(2)]
data = []
for ind0,state_i in enumerate(state_list):
measure_no_noise = ((matrix_subspace @ state_i) @ state_i.conj()).real
for noise_rate in tqdm(noise_rate_list, desc=f'state-{ind0}'):
for _ in range(num_random):
tmp0 = np_rng.normal(size=len(matrix_subspace))
noise = tmp0 * (noise_rate/np.linalg.norm(tmp0))
tmp1,eps = numqi.unique_determine.density_matrix_recovery_SDP(matrix_subspace, measure_no_noise + noise, converge_eps=cvxpy_eps)
tmp2 = np.linalg.norm(tmp1 - state_i[:,np.newaxis]*state_i.conj(), ord='fro') #frob norm
data.append((np.linalg.norm(noise), eps, tmp2))
data = np.array(data).reshape(-1, len(noise_rate_list), num_random, 3).transpose(0,3,1,2)
# about 3 mins
# to save time, we choose a relative small value for 'num_repeat' and 'num_random'
dim_qudit = 6
num_random = 50
cvxpy_eps = 1e-6
noise_rate_list = np.logspace(-6, -3, 10)
matrix_subspace = numqi.unique_determine.get_chebshev_orthonormal(dim_qudit, alpha=np.pi/dim_qudit, with_computational_basis=True)
tag_ud,loss,model = numqi.unique_determine.check_UD('uda', matrix_subspace, num_repeat=5,
converge_tol=1e-10, early_stop_threshold=1e-10, dtype='float64', tag_print=2, return_model=True)
matH = model.matH.numpy().copy()
EVL,EVC = np.linalg.eigh(matH)
state_list = [EVC[:,0]] + [numqi.random.rand_haar_state(dim_qudit) for _ in range(2)]
data = []
for ind0,state_i in enumerate(state_list):
measure_no_noise = ((matrix_subspace @ state_i) @ state_i.conj()).real
for noise_rate in tqdm(noise_rate_list, desc=f'state-{ind0}'):
for _ in range(num_random):
tmp0 = np_rng.normal(size=len(matrix_subspace))
noise = tmp0 * (noise_rate/np.linalg.norm(tmp0))
tmp1,eps = numqi.unique_determine.density_matrix_recovery_SDP(matrix_subspace, measure_no_noise + noise, converge_eps=cvxpy_eps)
tmp2 = np.linalg.norm(tmp1 - state_i[:,np.newaxis]*state_i.conj(), ord='fro') #frob norm
data.append((np.linalg.norm(noise), eps, tmp2))
data = np.array(data).reshape(-1, len(noise_rate_list), num_random, 3).transpose(0,3,1,2)
[round=0] min(f)=3.2235452508725745e-06, current(f)=3.2235452508725745e-06
[round=1] min(f)=3.2235452508725745e-06, current(f)=4.8004268702741686e-06
[round=2] min(f)=3.2235452508725745e-06, current(f)=9.534021790661999e-05
[round=3] min(f)=2.8737687585123248e-06, current(f)=2.8737687585123248e-06
[round=4] min(f)=2.8737687585123248e-06, current(f)=9.596716561364042e-05 [13.1] 1/1/1
In [5]:
Copied!
fig,(ax0,ax1) = plt.subplots(1, 2, figsize=(8,4))
for ind0 in range(data.shape[0]):
tmp0= noise_rate_list
tmp1 = data[ind0,1]
ax0.fill_between(tmp0, tmp1.min(axis=1), tmp1.max(axis=1), alpha=0.2, color=cp_tableau[ind0])
tmp2 = r'$\psi_-$' if ind0==0 else r'random $\sigma_'+f'{ind0}$'
ax0.plot(tmp0, tmp1.mean(axis=1), color=cp_tableau[ind0], label=tmp2)
tmp1 = data[ind0,2] / (data[ind0,0] + data[ind0,1])
ax1.fill_between(tmp0, tmp1.min(axis=1), tmp1.max(axis=1), alpha=0.2, color=cp_tableau[ind0])
ax1.plot(tmp0, tmp1.mean(axis=1), color=cp_tableau[ind0])
ax0.set_ylabel(r'$\epsilon$')
ax1.set_ylabel(r'$\frac{||Y-\sigma||_F}{\epsilon+||f||_2}$')
ax1.axhline(1/np.sqrt(loss), color='r', linestyle='--', label=r'$1/\sqrt{\mathcal{L}}$')
fig.suptitle(f'5PB(d={dim_qudit})')
ax0.legend()
for ax in [ax0,ax1]:
ax.set_xlabel(r'$||f||_2$')
ax.set_xscale('log')
ax.set_yscale('log')
ax1.yaxis.tick_right()
ax1.yaxis.set_label_position('right')
fig.tight_layout()
fig,(ax0,ax1) = plt.subplots(1, 2, figsize=(8,4))
for ind0 in range(data.shape[0]):
tmp0= noise_rate_list
tmp1 = data[ind0,1]
ax0.fill_between(tmp0, tmp1.min(axis=1), tmp1.max(axis=1), alpha=0.2, color=cp_tableau[ind0])
tmp2 = r'$\psi_-$' if ind0==0 else r'random $\sigma_'+f'{ind0}$'
ax0.plot(tmp0, tmp1.mean(axis=1), color=cp_tableau[ind0], label=tmp2)
tmp1 = data[ind0,2] / (data[ind0,0] + data[ind0,1])
ax1.fill_between(tmp0, tmp1.min(axis=1), tmp1.max(axis=1), alpha=0.2, color=cp_tableau[ind0])
ax1.plot(tmp0, tmp1.mean(axis=1), color=cp_tableau[ind0])
ax0.set_ylabel(r'$\epsilon$')
ax1.set_ylabel(r'$\frac{||Y-\sigma||_F}{\epsilon+||f||_2}$')
ax1.axhline(1/np.sqrt(loss), color='r', linestyle='--', label=r'$1/\sqrt{\mathcal{L}}$')
fig.suptitle(f'5PB(d={dim_qudit})')
ax0.legend()
for ax in [ax0,ax1]:
ax.set_xlabel(r'$||f||_2$')
ax.set_xscale('log')
ax.set_yscale('log')
ax1.yaxis.tick_right()
ax1.yaxis.set_label_position('right')
fig.tight_layout()
Stability of Pauli UD Schemes¶
In [6]:
Copied!
# about 3 mins
num_qubit = 3
num_random = 50
cvxpy_eps = 1e-6
noise_rate_list = np.logspace(-6, -3, 6)
pauli_index = numqi.unique_determine.load_pauli_ud_example(num_qubit)[0] #sorted by length
pauli_group = numqi.gate.get_pauli_group(num_qubit, use_sparse=True)
matrix_subspace = numqi.unique_determine.get_matrix_list_indexing(pauli_group, pauli_index)
tag_ud,loss,model = numqi.unique_determine.check_UD('uda', matrix_subspace, num_repeat=5,
converge_tol=1e-10, early_stop_threshold=1e-10, dtype='float64', tag_print=2, return_model=True)
matH = model.matH.numpy().copy()
EVL,EVC = np.linalg.eigh(matH)
state_list = [EVC[:,0]] + [numqi.random.rand_haar_state(2**num_qubit) for _ in range(2)]
matrix_subspace_full = np.stack([x.toarray() for x in matrix_subspace])
data = []
for ind0,state_i in enumerate(state_list):
measure_no_noise = ((matrix_subspace_full @ state_i) @ state_i.conj()).real
for noise_rate in tqdm(noise_rate_list, desc=f'state-{ind0}'):
for _ in range(num_random):
tmp0 = np_rng.normal(size=len(matrix_subspace_full))
noise = tmp0 * (noise_rate/np.linalg.norm(tmp0))
tmp1,eps = numqi.unique_determine.density_matrix_recovery_SDP(matrix_subspace_full, measure_no_noise + noise, converge_eps=cvxpy_eps)
tmp2 = np.linalg.norm(tmp1 - state_i[:,np.newaxis]*state_i.conj(), ord='fro') #frob norm
data.append((np.linalg.norm(noise), eps, tmp2))
data = np.array(data).reshape(-1, len(noise_rate_list), num_random, 3).transpose(0,3,1,2)
# about 3 mins
num_qubit = 3
num_random = 50
cvxpy_eps = 1e-6
noise_rate_list = np.logspace(-6, -3, 6)
pauli_index = numqi.unique_determine.load_pauli_ud_example(num_qubit)[0] #sorted by length
pauli_group = numqi.gate.get_pauli_group(num_qubit, use_sparse=True)
matrix_subspace = numqi.unique_determine.get_matrix_list_indexing(pauli_group, pauli_index)
tag_ud,loss,model = numqi.unique_determine.check_UD('uda', matrix_subspace, num_repeat=5,
converge_tol=1e-10, early_stop_threshold=1e-10, dtype='float64', tag_print=2, return_model=True)
matH = model.matH.numpy().copy()
EVL,EVC = np.linalg.eigh(matH)
state_list = [EVC[:,0]] + [numqi.random.rand_haar_state(2**num_qubit) for _ in range(2)]
matrix_subspace_full = np.stack([x.toarray() for x in matrix_subspace])
data = []
for ind0,state_i in enumerate(state_list):
measure_no_noise = ((matrix_subspace_full @ state_i) @ state_i.conj()).real
for noise_rate in tqdm(noise_rate_list, desc=f'state-{ind0}'):
for _ in range(num_random):
tmp0 = np_rng.normal(size=len(matrix_subspace_full))
noise = tmp0 * (noise_rate/np.linalg.norm(tmp0))
tmp1,eps = numqi.unique_determine.density_matrix_recovery_SDP(matrix_subspace_full, measure_no_noise + noise, converge_eps=cvxpy_eps)
tmp2 = np.linalg.norm(tmp1 - state_i[:,np.newaxis]*state_i.conj(), ord='fro') #frob norm
data.append((np.linalg.norm(noise), eps, tmp2))
data = np.array(data).reshape(-1, len(noise_rate_list), num_random, 3).transpose(0,3,1,2)
[round=0] min(f)=0.9809962203279491, current(f)=0.9809962203279491
[round=1] min(f)=0.9809962203279491, current(f)=1.0243478958105374
[round=2] min(f)=0.9809962203279491, current(f)=1.0243463397777144
[round=3] min(f)=0.9809962203279491, current(f)=1.0000018948024187
[round=4] min(f)=0.9809962203279491, current(f)=1.0000006290249053 [7.3] 1/1/1
In [7]:
Copied!
fig,(ax0,ax1) = plt.subplots(1, 2, figsize=(8,4))
for ind0 in range(data.shape[0]):
tmp0= noise_rate_list
tmp1 = data[ind0,1]
ax0.fill_between(tmp0, tmp1.min(axis=1), tmp1.max(axis=1), alpha=0.2, color=cp_tableau[ind0])
tmp2 = r'$\psi_-$' if ind0==0 else r'random $\sigma_'+f'{ind0}$'
ax0.plot(tmp0, tmp1.mean(axis=1), color=cp_tableau[ind0], label=tmp2)
tmp1 = data[ind0,2] / (data[ind0,0] + data[ind0,1])
ax1.fill_between(tmp0, tmp1.min(axis=1), tmp1.max(axis=1), alpha=0.2, color=cp_tableau[ind0])
ax1.plot(tmp0, tmp1.mean(axis=1), color=cp_tableau[ind0])
ax0.set_ylabel(r'$\epsilon$')
ax1.set_ylabel(r'$\frac{||Y-\sigma||_F}{\epsilon+||f||_2}$')
ax1.axhline(1/np.sqrt(loss), color='r', linestyle='--', label=r'$1/\sqrt{\mathcal{L}}$')
fig.suptitle(f'Pauli UD(n={num_qubit})')
ax0.legend()
for ax in [ax0,ax1]:
ax.set_xlabel(r'$||f||_2$')
ax.set_xscale('log')
ax.set_yscale('log')
ax1.yaxis.tick_right()
ax1.yaxis.set_label_position('right')
fig.tight_layout()
fig,(ax0,ax1) = plt.subplots(1, 2, figsize=(8,4))
for ind0 in range(data.shape[0]):
tmp0= noise_rate_list
tmp1 = data[ind0,1]
ax0.fill_between(tmp0, tmp1.min(axis=1), tmp1.max(axis=1), alpha=0.2, color=cp_tableau[ind0])
tmp2 = r'$\psi_-$' if ind0==0 else r'random $\sigma_'+f'{ind0}$'
ax0.plot(tmp0, tmp1.mean(axis=1), color=cp_tableau[ind0], label=tmp2)
tmp1 = data[ind0,2] / (data[ind0,0] + data[ind0,1])
ax1.fill_between(tmp0, tmp1.min(axis=1), tmp1.max(axis=1), alpha=0.2, color=cp_tableau[ind0])
ax1.plot(tmp0, tmp1.mean(axis=1), color=cp_tableau[ind0])
ax0.set_ylabel(r'$\epsilon$')
ax1.set_ylabel(r'$\frac{||Y-\sigma||_F}{\epsilon+||f||_2}$')
ax1.axhline(1/np.sqrt(loss), color='r', linestyle='--', label=r'$1/\sqrt{\mathcal{L}}$')
fig.suptitle(f'Pauli UD(n={num_qubit})')
ax0.legend()
for ax in [ax0,ax1]:
ax.set_xlabel(r'$||f||_2$')
ax.set_xscale('log')
ax.set_yscale('log')
ax1.yaxis.tick_right()
ax1.yaxis.set_label_position('right')
fig.tight_layout()
