Average error#
import quantum_simulation_recipe as qsr
from quantum_simulation_recipe import spin_ham
H = spin_ham.Nearest_Neighbour_1d(4)
H.ham
SparsePauliOp(['IIXX', 'IXXI', 'XXII', 'IIYY', 'IYYI', 'YYII', 'IIZZ', 'IZZI', 'ZZII', 'IIIX', 'IIXI', 'IXII', 'XIII'],
coeffs=[1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j,
1. +0.j, 0.2+0.j, 0.2+0.j, 0.2+0.j, 0.2+0.j])
State: Worst-case VS average-case Bound#
Worst initial state for Trotter#
from qiskit.quantum_info import random_unitary, random_hermitian
from quantum_simulation_recipe.spin import Nearest_Neighbour_1d
from quantum_simulation_recipe.trotter import pf, expH
import numpy as np
def error_operator(H_list, t=2.453, r=1):
appro_U, exact_U = pf(H_list, t, r, return_exact=True)
return exact_U - appro_U
MFI = Nearest_Neighbour_1d(2, hx=0.8090, hy=0.9045, Jx=1, pbc=False)
TFI = Nearest_Neighbour_1d(2, Jz=2, hx=1, pbc=False)
TFI_H_list = [h.to_matrix() for h in TFI.ham_par]
# H_list = [h.to_matrix() for h in TFI.ham_xyz]
MFI_H_list = [h.to_matrix() for h in MFI.ham_par]
RND_H_list = [random_hermitian(4), random_hermitian(4)]
TFI_E = error_operator(TFI_H_list)
MFI_E = error_operator(MFI_H_list)
RND_E = error_operator(RND_H_list)
# D = expH(H1+H2, 1) - expH(H1, 1)@expH(H2, 1)
print('E of MFI is normal: ', np.allclose(MFI_E@MFI_E.conj().T, MFI_E.conj().T@MFI_E))
print('E of TFI is normal: ', np.allclose(TFI_E@TFI_E.conj().T, TFI_E.conj().T@TFI_E))
print('E of TFI is Hermitian: ', np.allclose(TFI_E, TFI_E.conj().T))
print('E of TFI is anti-Hermitian: ', np.allclose(TFI_E, -TFI_E.conj().T))
print('E of TFI is unitary: ', np.allclose(np.dot(TFI_E, TFI_E.conj().T), np.eye(TFI_E.shape[0])))
E of MFI is normal: False
E of TFI is normal: True
E of TFI is Hermitian: False
E of TFI is anti-Hermitian: False
E of TFI is unitary: False
D = MFI_E
# D = TFI_E
# D = RND_E
print('DD* is Hermitian (trivial): ', np.allclose(D@D.conj().T, (D@D.conj().T).conj().T))
print('eigenvals of D (by eigh): ', np.linalg.eigh(D)[0])
print('eigenvals of D*D (by eigh): \n', np.linalg.eigh(D@D.conj().T)[0])
# print(np.linalg.eigh(D.conj().T@D)[0])
print('eigenvals of D*D (by eig, not sorted): \n', np.linalg.eig(D@D.conj().T)[0])
DD* is Hermitian (trivial): True
eigenvals of D (by eigh): [-0.79155 -0.37126 0.02682 1.54871]
eigenvals of D*D (by eigh):
[0.01403 0.11221 2.4417 2.85321]
eigenvals of D*D (by eig, not sorted):
[2.85321-0.j 2.4417 -0.j 0.11221+0.j 0.01403-0.j]
print('D is normal: ', np.allclose(D@D.conj().T, D.conj().T@D))
# print(D@D.conj().T-D.conj().T@D)
print('---------eigen of D D*')
print(np.linalg.eigh(D @ D.conj().T)[0])
# print(np.linalg.eigh(D @ D.conj().T)[1])
print('---------eigen of D* D')
print(np.linalg.eigh(D.conj().T @ D)[0])
# print(np.linalg.eigh(D.conj().T @ D)[1])
print('DD*!=D*D, but they have the same eigenvalues (similar??).')
print(np.dot(np.linalg.eigh(D.conj().T @ D)[1][:, 0], np.linalg.eigh(D @ D.conj().T)[1][:, 0].conj().T))
D is normal: False
---------eigen of D D*
[0.01403 0.11221 2.4417 2.85321]
---------eigen of D* D
[0.01403 0.11221 2.4417 2.85321]
DD*!=D*D, but they have the same eigenvalues (similar??).
(0.5497812019484642+0.050036165227177194j)
Worst-case: singular value#
np.set_printoptions(precision=5, suppress=True)
print('spetral norm: ', np.linalg.norm(D, ord=2))
U, S, V = np.linalg.svd(D)
print("S: ", S)
print('square of singular values: \n', sorted(S**2))
print("U: \n", U)
print("V: \n", V)
# init_st = V[0].T
init_st_svd = V[0].conj().T
print("init_st_svd: ", init_st_svd)
spetral norm: 1.6891439100095893
S: [1.68914 1.56259 0.33498 0.11845]
square of singular values:
[0.014029287629003814, 0.11220908840624161, 2.441700759298198, 2.8532071487224835]
U:
[[-0.20572+0.57631j -0.4436 +0.04446j -0.40134-0.40178j 0.31204+0.08321j]
[ 0.05161-0.00535j -0.09772-0.78246j -0.3361 +0.47573j 0.18116+0.05845j]
[ 0.09959-0.40388j -0.05016+0.27609j -0.06563+0.08953j 0.55808+0.65149j]
[ 0.05193+0.66867j 0.31645-0.02386j 0.4988 +0.27766j 0.16947+0.308j ]]
V:
[[-0.67068+0.j 0.39496+0.13057j 0.00134+0.05187j -0.55865-0.24973j]
[-0.31735+0.j 0.07078-0.27154j 0.03861+0.78759j 0.44569-0.01098j]
[-0.57087+0.j 0.0138 -0.11015j 0.06228-0.57914j 0.54609+0.15585j]
[-0.35154+0.j -0.83982+0.17489j -0.13854+0.13054j -0.22333+0.23328j]]
init_st_svd: [-0.67068-0.j 0.39496-0.13057j 0.00134-0.05187j -0.55865+0.24973j]
the worst-input is the U, V or eigenstate of largest singular value?
If E is normal e.g. TFI, they are the same
Otherwise e.g. MFI, the worst input (largest singular value) is V[0].conj().T
print('eigen values: ', np.linalg.eigh(D @ D.conj().T)[0])
print("squre of singualr values: ", sorted(S**2))
print(np.linalg.eigh(D @ D.conj().T)[1])
print(D@np.linalg.eigh(D @ D.conj().T)[1][:,-1])
print('eigenvector with largest eigval: \n', np.divide(D @ D.conj().T@np.linalg.eigh(D @ D.conj().T)[1][:,-1], np.linalg.eigh(D @ D.conj().T)[0][-1]))
# print(np.linalg.norm(np.linalg.eigh(D @ D.conj().T)[1][:,-1])) verify normalized vector
print(np.linalg.norm(D@np.linalg.eigh(D @ D.conj().T)[1][:,-1])) # .conj().T
print('achieve spectral norm (square root of largest eigenval): ', np.linalg.norm(D@init_st_svd))
print(np.linalg.norm(D@V[0]))
print(np.linalg.norm(D@U[:, 0].conj().T))
eigen values: [0.01403 0.11221 2.4417 2.85321]
squre of singualr values: [0.014029287629003814, 0.11220908840624161, 2.441700759298198, 2.8532071487224835]
[[-0.32294+0.j 0.56789+0.j 0.44582-0.j -0.61192+0.j ]
[-0.1901 -0.0098j -0.09905-0.574j 0.0192 +0.7883j 0.02239+0.0468j ]
[-0.7071 -0.4857j -0.01695-0.10971j 0.07744-0.26971j 0.41386-0.04199j]
[-0.2431 -0.25393j -0.54895+0.15667j -0.31725-0.00782j -0.61229+0.27371j]]
[-0.36482+0.50027j 0.66319-0.04171j 0.00876-0.58929j 0.05976+1.13047j]
eigenvector with largest eigval:
[-0.61192-0.j 0.02239+0.0468j 0.41386-0.04199j -0.61229+0.27371j]
1.5664627453249478
achieve spectral norm (square root of largest eigenval): 1.689143910009589
1.5873420245042686
1.350723692946984