Skip to content

Commit c5d9ff3

Browse files
rroodllrroodll
authored
Merge pull request #41 from isi-usc-edu/31-one-norm-does-not-account-for-double-factorization-truncation-threshold
truncate alpha for doublefactorized encoding
2 parents 7e44c44 + dc36f32 commit c5d9ff3

File tree

5 files changed

+216
-25
lines changed

5 files changed

+216
-25
lines changed

‎src/pyLIQTR/BlockEncodings/DoubleFactorized.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -103,7 +103,7 @@ def __init__(self,ProblemInstance,df_error_threshold:float=1e-3,sf_error_thresho
103103
@cached_property
104104
def alpha(self):
105105
"""returns the double factorized Hamiltonian norm."""
106-
return self.PI.get_alpha(encoding='DF',df_error_threshold=self.df_error_threshold,sf_error_threshold=self.sf_error_threshold)
106+
return self.PI.get_alpha(encoding='DF',df_cutoffs=self.Xi_l_data)
107107

108108
@cached_property
109109
def control_registers(self) -> Tuple[Register]:

‎src/pyLIQTR/BlockEncodings/tests/test_DoubleFactorized.py‎

Lines changed: 73 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -5,22 +5,88 @@
55
import pytest
66
import platform
77
import cirq
8+
import numpy as np
89
from pyLIQTR.ProblemInstances.getInstance import *
910
from pyLIQTR.BlockEncodings.getEncoding import *
1011
from pyLIQTR.utils.printing import openqasm
1112
from pyLIQTR.utils.resource_analysis import estimate_resources
12-
from openfermion.chem import MolecularData
13+
from openfermion import InteractionOperator
1314

14-
@pytest.mark.skipif(platform.system() == 'Windows', reason = "cannot run this test on Windows")
15+
@pytest.mark.skipif(platform.system() == 'Windows', reason = "pyscf not supported on Windows")
1516
class TestDoubleFactorizedEncoding:
1617

1718
@pytest.fixture(scope="class")
1819
def h2_instance(self):
19-
from openfermionpyscf import run_pyscf
20-
mol_data = MolecularData([('H', (0.0, 0.0, 0.63164)), ('H', (0.0, 0.0, 1.76836))],\
21-
'sto-3g', 1.0, 0, 'H2')
22-
mol = run_pyscf(mol_data, run_scf=1, run_mp2=0, run_cisd=0, run_ccsd=0, run_fci=0, verbose=0)
23-
mol_instance = getInstance('ChemicalHamiltonian',mol_ham=mol.get_molecular_hamiltonian(),mol_name='H2')
20+
obt = np.array([[-1.04689337, 0.0, 0.0, 0.0],
21+
[ 0.0, -1.04689337, 0.0, 0.0],
22+
[ 0.0, 0.0, -0.62248428, 0.0],
23+
[ 0.0, 0.0, 0.0, -0.62248428]])
24+
tbt = np.array([[[[0.30156089, 0.0, 0.0, 0.0],
25+
[0.0, 0.0, 0.0, 0.0],
26+
[0.0, 0.0, 0.10281136, 0.0],
27+
[0.0, 0.0, 0.0, 0.0]],
28+
[[0.0, 0.0, 0.0, 0.0],
29+
[0.30156089, 0.0, 0.0, 0.0],
30+
[0.0, 0.0, 0.0, 0.0],
31+
[0.0, 0.0, 0.10281136, 0.0]],
32+
[[0.0, 0.0, 0.10281136, 0.0],
33+
[0.0, 0.0, 0.0, 0.0],
34+
[0.3011551 , 0.0, 0.0, 0.0],
35+
[0.0, 0.0, 0.0, 0.0]],
36+
[[0.0, 0.0, 0.0, 0.0],
37+
[0.0, 0.0, 0.10281136, 0.0],
38+
[0.0, 0.0, 0.0, 0.0],
39+
[0.3011551 , 0.0, 0.0, 0.0]]],
40+
[[[0.0, 0.30156089, 0.0, 0.0],
41+
[0.0, 0.0, 0.0, 0.0],
42+
[0.0, 0.0, 0.0, 0.10281136],
43+
[0.0, 0.0, 0.0, 0.0]],
44+
[[0.0, 0.0, 0.0, 0.0],
45+
[0.0, 0.30156089, 0.0, 0.0],
46+
[0.0, 0.0, 0.0, 0.0],
47+
[0.0, 0.0, 0.0, 0.10281136]],
48+
[[0.0, 0.0, 0.0, 0.10281136],
49+
[0.0, 0.0, 0.0, 0.0],
50+
[0.0, 0.3011551 , 0.0, 0.0],
51+
[0.0, 0.0, 0.0, 0.0]],
52+
[[0.0, 0.0, 0.0, 0.0],
53+
[0.0, 0.0, 0.0, 0.10281136],
54+
[0.0, 0.0, 0.0, 0.0],
55+
[0.0, 0.3011551 , 0.0, 0.0]]],
56+
[[[0.0, 0.0, 0.3011551 , 0.0],
57+
[0.0, 0.0, 0.0, 0.0],
58+
[0.10281136, 0.0, 0.0, 0.0],
59+
[0.0, 0.0, 0.0, 0.0]],
60+
[[0.0, 0.0, 0.0, 0.0],
61+
[0.0, 0.0, 0.3011551 , 0.0],
62+
[0.0, 0.0, 0.0, 0.0],
63+
[0.10281136, 0.0, 0.0, 0.0]],
64+
[[0.10281136, 0.0, 0.0, 0.0],
65+
[0.0, 0.0, 0.0, 0.0],
66+
[0.0, 0.0, 0.31597877, 0.0],
67+
[0.0, 0.0, 0.0, 0.0]],
68+
[[0.0, 0.0, 0.0, 0.0],
69+
[0.10281136, 0.0, 0.0, 0.0],
70+
[0.0, 0.0, 0.0, 0.0],
71+
[0.0, 0.0, 0.31597877, 0.0]]],
72+
[[[0.0, 0.0, 0.0, 0.3011551 ],
73+
[0.0, 0.0, 0.0, 0.0],
74+
[0.0, 0.10281136, 0.0, 0.0],
75+
[0.0, 0.0, 0.0, 0.0]],
76+
[[0.0, 0.0, 0.0, 0.0],
77+
[0.0, 0.0, 0.0, 0.3011551 ],
78+
[0.0, 0.0, 0.0, 0.0],
79+
[0.0, 0.10281136, 0.0, 0.0]],
80+
[[0.0, 0.10281136, 0.0, 0.0],
81+
[0.0, 0.0, 0.0, 0.0],
82+
[0.0, 0.0, 0.0, 0.31597877],
83+
[0.0, 0.0, 0.0, 0.0]],
84+
[[0.0, 0.0, 0.0, 0.0],
85+
[0.0, 0.10281136, 0.0, 0.0],
86+
[0.0, 0.0, 0.0, 0.0],
87+
[0.0, 0.0, 0.0, 0.31597877]]]])
88+
h2_op = InteractionOperator(constant=0.4655299554155818,one_body_tensor=obt,two_body_tensor=tbt)
89+
mol_instance = getInstance('ChemicalHamiltonian',mol_ham=h2_op,mol_name='H2')
2490
return mol_instance
2591

2692
@pytest.fixture(scope="class")

‎src/pyLIQTR/ProblemInstances/ChemicalHamiltonian.py‎

Lines changed: 36 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -36,8 +36,9 @@ def __init__(self, mol_ham: InteractionOperator, mol_name=None, **kwargs):
3636
juliapkg.resolve()
3737

3838
# Now start the real initialization
39-
self._mol_ham = mol_ham
40-
self._mol_name = mol_name
39+
self._mol_ham = mol_ham
40+
self._mol_name = mol_name
41+
self._df_cutoffs = dict()
4142

4243
super(ProblemInstance, self).__init__(**kwargs)
4344

@@ -53,11 +54,6 @@ def __str__(self):
5354
def n_qubits(self):
5455
return self._mol_ham.n_qubits
5556

56-
57-
# # @property
58-
# def n_terms(self,**kwargs):
59-
# return len(list(self._terms_jw))
60-
6157
@cached_property
6258
def terms_jw(self):
6359
return jordan_wigner(self._mol_ham).terms
@@ -110,14 +106,29 @@ def hamiltonian_tensors(self):
110106

111107
return {'h0':h0, 'one_body_tensor':one_body_tensor, 'two_body_tensor':two_body_tensor}
112108

109+
@cached_property
110+
def obt_fragment(self):
111+
from pyLIQTR.utils.df_utils import to_OBF
112+
_, obt, _ = self.hamiltonian_tensors.values()
113+
obt_frag = to_OBF(obt)
114+
return obt_frag
115+
113116
def DF_fragments(self,sf_error_threshold):
114117
from pyLIQTR.utils.df_utils import DF_decomposition
115118
h0, one_body_tensor, two_body_tensor = self.hamiltonian_tensors.values()
116119
return DF_decomposition(h0, one_body_tensor, two_body_tensor,tol=sf_error_threshold)
117120

118-
def get_alpha(self,encoding:str='PauliLCU',df_error_threshold=None,sf_error_threshold=None):
121+
def df_cutoffs(self,df_error_threshold:float=1e-3,sf_error_threshold:float=1e-10):
122+
error_pair = (df_error_threshold,sf_error_threshold)
123+
if error_pair in self._df_cutoffs:
124+
return self._df_cutoffs[error_pair]
125+
else:
126+
__,__,__,__ = self.yield_DF_Info(df_error_threshold=df_error_threshold,sf_error_threshold=sf_error_threshold)
127+
return self._df_cutoffs[error_pair]
128+
129+
def get_alpha(self,encoding:str='PauliLCU',sf_error_threshold:float=1e-10,df_error_threshold:float=None,df_cutoffs:list=None):
119130
if encoding == 'PauliLCU':
120-
return(self._ops.get_alpha())
131+
return(self._ops.get_coeff_norm())
121132
elif encoding == 'DF':
122133
from pyLIQTR.utils.df_utils import to_OBF
123134
if 'sphinx' not in sys.modules:
@@ -126,15 +137,24 @@ def get_alpha(self,encoding:str='PauliLCU',df_error_threshold=None,sf_error_thre
126137
jl.seval('Pkg.add("QuantumMAMBO")')
127138
jl.seval("using QuantumMAMBO")
128139
mambo = jl.QuantumMAMBO
140+
141+
if df_cutoffs is None:
142+
if df_error_threshold is None:
143+
df_cutoffs = self.df_cutoffs(df_error_threshold=1e-3,sf_error_threshold=sf_error_threshold)
144+
else:
145+
df_cutoffs = self.df_cutoffs(df_error_threshold=df_error_threshold,sf_error_threshold=sf_error_threshold)
146+
else:
147+
if df_error_threshold is not None:
148+
raise ValueError("provide only df_error_threshold or df_cutoffs")
149+
129150
_, one_body_tensor, two_body_tensor = self.hamiltonian_tensors.values()
130151
DF_frags = self.DF_fragments(sf_error_threshold)
131152
one_body_correction = 2*sum([two_body_tensor[:,:,r,r] for r in range(two_body_tensor.shape[0])])
132153
one_body_fragment = to_OBF(one_body_tensor + one_body_correction)
133154
lambdaTprime = mambo.OBF_L1(one_body_fragment)
134155
lambdaDF = 0.0
135-
for frag in DF_frags:
136-
lambdaDF += mambo.DF_L1(frag)
137-
# TODO: lambdaDF should exclude terms based on error threshold
156+
for l,frag in enumerate(DF_frags):
157+
lambdaDF += 0.5 * abs(frag.coeff) * ((sum(np.abs(frag.C.λ[:df_cutoffs[l]])))**2)
138158
return lambdaTprime + lambdaDF
139159

140160

@@ -176,7 +196,7 @@ def yield_PauliLCU_Info(self,return_as='arrays',do_pad=0,pad_value=1.0):
176196

177197

178198
def yield_DF_Info(self, df_error_threshold: float,sf_error_threshold:float=1e-10):
179-
from pyLIQTR.utils.df_utils import to_OBF, U_to_Givens, calc_xi
199+
from pyLIQTR.utils.df_utils import U_to_Givens, calc_xi
180200

181201
_, obt, tbt = self.hamiltonian_tensors.values()
182202

@@ -187,7 +207,7 @@ def yield_DF_Info(self, df_error_threshold: float,sf_error_threshold:float=1e-10
187207
mus_mat = np.zeros(shape = (num_frags + 1, num_orbs))
188208
thetas_tsr = np.zeros(shape = (num_frags + 1, num_orbs, num_orbs - 1))
189209

190-
obt_frag = to_OBF(obt)
210+
obt_frag = self.obt_fragment
191211
for i in range(len(obt_frag.C.λ)):
192212
mus_mat[0][i] = obt_frag.C.λ[i]
193213

@@ -228,9 +248,9 @@ def yield_DF_Info(self, df_error_threshold: float,sf_error_threshold:float=1e-10
228248
f_p_vals.append(abs(k))
229249
f_p_full.append([f_p_signs, f_p_vals])
230250

231-
xi = calc_xi(f_p_abs, df_error_threshold)
251+
self._df_cutoffs[(df_error_threshold,sf_error_threshold)] = calc_xi(f_p_abs, df_error_threshold)
232252

233-
return T_prime_full, f_p_full, xi, thetas_tsr
253+
return T_prime_full, f_p_full, self._df_cutoffs[(df_error_threshold,sf_error_threshold)], thetas_tsr
234254

235255

236256

‎src/pyLIQTR/ProblemInstances/tests/test_ChemicalHamiltonian.py‎

Lines changed: 105 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
"""
2+
Copyright (c) 2024 Massachusetts Institute of Technology
3+
SPDX-License-Identifier: BSD-2-Clause
4+
"""
5+
import pytest
6+
import platform
7+
import numpy as np
8+
from openfermion import InteractionOperator
9+
from pyLIQTR.ProblemInstances.getInstance import *
10+
11+
@pytest.mark.skipif(platform.system() == 'Windows', reason = "pyscf not supported on Windows")
12+
class TestChemicalHamiltonian:
13+
14+
@pytest.fixture(scope="class")
15+
def h2_instance(self):
16+
obt = np.array([[-1.04689337, 0.0, 0.0, 0.0],
17+
[ 0.0, -1.04689337, 0.0, 0.0],
18+
[ 0.0, 0.0, -0.62248428, 0.0],
19+
[ 0.0, 0.0, 0.0, -0.62248428]])
20+
tbt = np.array([[[[0.30156089, 0.0, 0.0, 0.0],
21+
[0.0, 0.0, 0.0, 0.0],
22+
[0.0, 0.0, 0.10281136, 0.0],
23+
[0.0, 0.0, 0.0, 0.0]],
24+
[[0.0, 0.0, 0.0, 0.0],
25+
[0.30156089, 0.0, 0.0, 0.0],
26+
[0.0, 0.0, 0.0, 0.0],
27+
[0.0, 0.0, 0.10281136, 0.0]],
28+
[[0.0, 0.0, 0.10281136, 0.0],
29+
[0.0, 0.0, 0.0, 0.0],
30+
[0.3011551 , 0.0, 0.0, 0.0],
31+
[0.0, 0.0, 0.0, 0.0]],
32+
[[0.0, 0.0, 0.0, 0.0],
33+
[0.0, 0.0, 0.10281136, 0.0],
34+
[0.0, 0.0, 0.0, 0.0],
35+
[0.3011551 , 0.0, 0.0, 0.0]]],
36+
[[[0.0, 0.30156089, 0.0, 0.0],
37+
[0.0, 0.0, 0.0, 0.0],
38+
[0.0, 0.0, 0.0, 0.10281136],
39+
[0.0, 0.0, 0.0, 0.0]],
40+
[[0.0, 0.0, 0.0, 0.0],
41+
[0.0, 0.30156089, 0.0, 0.0],
42+
[0.0, 0.0, 0.0, 0.0],
43+
[0.0, 0.0, 0.0, 0.10281136]],
44+
[[0.0, 0.0, 0.0, 0.10281136],
45+
[0.0, 0.0, 0.0, 0.0],
46+
[0.0, 0.3011551 , 0.0, 0.0],
47+
[0.0, 0.0, 0.0, 0.0]],
48+
[[0.0, 0.0, 0.0, 0.0],
49+
[0.0, 0.0, 0.0, 0.10281136],
50+
[0.0, 0.0, 0.0, 0.0],
51+
[0.0, 0.3011551 , 0.0, 0.0]]],
52+
[[[0.0, 0.0, 0.3011551 , 0.0],
53+
[0.0, 0.0, 0.0, 0.0],
54+
[0.10281136, 0.0, 0.0, 0.0],
55+
[0.0, 0.0, 0.0, 0.0]],
56+
[[0.0, 0.0, 0.0, 0.0],
57+
[0.0, 0.0, 0.3011551 , 0.0],
58+
[0.0, 0.0, 0.0, 0.0],
59+
[0.10281136, 0.0, 0.0, 0.0]],
60+
[[0.10281136, 0.0, 0.0, 0.0],
61+
[0.0, 0.0, 0.0, 0.0],
62+
[0.0, 0.0, 0.31597877, 0.0],
63+
[0.0, 0.0, 0.0, 0.0]],
64+
[[0.0, 0.0, 0.0, 0.0],
65+
[0.10281136, 0.0, 0.0, 0.0],
66+
[0.0, 0.0, 0.0, 0.0],
67+
[0.0, 0.0, 0.31597877, 0.0]]],
68+
[[[0.0, 0.0, 0.0, 0.3011551 ],
69+
[0.0, 0.0, 0.0, 0.0],
70+
[0.0, 0.10281136, 0.0, 0.0],
71+
[0.0, 0.0, 0.0, 0.0]],
72+
[[0.0, 0.0, 0.0, 0.0],
73+
[0.0, 0.0, 0.0, 0.3011551 ],
74+
[0.0, 0.0, 0.0, 0.0],
75+
[0.0, 0.10281136, 0.0, 0.0]],
76+
[[0.0, 0.10281136, 0.0, 0.0],
77+
[0.0, 0.0, 0.0, 0.0],
78+
[0.0, 0.0, 0.0, 0.31597877],
79+
[0.0, 0.0, 0.0, 0.0]],
80+
[[0.0, 0.0, 0.0, 0.0],
81+
[0.0, 0.10281136, 0.0, 0.0],
82+
[0.0, 0.0, 0.0, 0.0],
83+
[0.0, 0.0, 0.0, 0.31597877]]]])
84+
h2_op = InteractionOperator(constant=0.4655299554155818,one_body_tensor=obt,two_body_tensor=tbt)
85+
mol_instance = getInstance('ChemicalHamiltonian',mol_ham=h2_op,mol_name='H2')
86+
return mol_instance
87+
88+
def test_ChemicalHamiltonian_alpha(self, h2_instance):
89+
lcu_alpha1 = h2_instance.get_alpha()
90+
assert lcu_alpha1 == 1.4527183600000002
91+
92+
lcu_alpha2 = h2_instance.get_alpha(encoding='PauliLCU')
93+
assert lcu_alpha2 == 1.4527183600000002
94+
95+
df_alpha = h2_instance.get_alpha(encoding='DF')
96+
assert df_alpha == 1.2543735419921236
97+
98+
df_alpha_pass_err = h2_instance.get_alpha(encoding='DF',df_error_threshold=0.99)
99+
assert df_alpha_pass_err == 0.8006020005033776
100+
101+
df_alpha_pass_xi = h2_instance.get_alpha(encoding='DF',df_cutoffs=[1,1])
102+
assert df_alpha_pass_xi == 0.6463849605033776
103+
104+
with pytest.raises(ValueError,match='provide only df_error_threshold or df_cutoffs'):
105+
alpha = h2_instance.get_alpha(encoding='DF',df_cutoffs=[1,1],df_error_threshold=1e-2)

‎src/pyLIQTR/_version.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1 @@
1-
__version__ = "1.3.4"
1+
__version__ = "1.3.5"

0 commit comments

Comments
 (0)