Skip to content

Commit 4c63f25

Browse files
SexyCarrotsSexyCarrots
committed
2 parents 970b35c + 80ecf2c commit 4c63f25

File tree

6 files changed

+474
-6
lines changed

6 files changed

+474
-6
lines changed

‎CHANGELOG.md

+11-3
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,14 @@
22

33
## Unreleased
44

5+
### Added
6+
7+
- add QuOperator convert tools which can convert MPO in the form of TensorNetwork and Quimb into MPO in the form of QuOperator
8+
9+
### Changed
10+
11+
- quantum Hamiltonian generation now support the direct return of numpy form matrix
12+
513
## 0.0.220318
614

715
### Added
@@ -16,7 +24,7 @@
1624

1725
- futher refactor VQNHE code in applications
1826

19-
- add alias ``sample`` for ``perfect_sampling`` method
27+
- add alias `sample` for `perfect_sampling` method
2028

2129
- optimize VQNHE pipeline for a more stable training loop (breaking changes in some APIs)
2230

@@ -32,9 +40,9 @@
3240

3341
- add MPO expectation template function for MPO evaluation on circuit
3442

35-
- add ``operator_expectation`` in templates.measurements for a unified expectation interface
43+
- add `operator_expectation` in templates.measurements for a unified expectation interface
3644

37-
- add ``templates.chems`` module for interface between tc and openfermion on quantum chemistry related tasks
45+
- add `templates.chems` module for interface between tc and openfermion on quantum chemistry related tasks
3846

3947
- add tc.Circuit to Qiskit QuantumCircuit transformation
4048

‎examples/hamiltonian_building.py

+45
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
1+
"""
2+
benchmark sparse hamiltonian building
3+
"""
4+
import time
5+
import numpy as np
6+
import quimb
7+
import tensorcircuit as tc
8+
9+
10+
nwires = 20
11+
# tc outperforms quimb in large nwires
12+
13+
14+
# 1. tc approach for TFIM
15+
16+
g = tc.templates.graphs.Line1D(nwires, pbc=False)
17+
time0 = time.time()
18+
h1 = tc.quantum.heisenberg_hamiltonian(
19+
g, hzz=1, hxx=0, hyy=0, hz=0, hx=-1, hy=0, sparse=True, numpy=True
20+
)
21+
time1 = time.time()
22+
23+
print("tc time: ", time1 - time0)
24+
25+
# 2. quimb approach for TFIM
26+
27+
builder = quimb.tensor.tensor_gen.SpinHam()
28+
# spin operator instead of Pauli matrix
29+
builder += 4, "Z", "Z"
30+
builder += -2, "X"
31+
time0 = time.time()
32+
h2 = builder.build_sparse(nwires)
33+
h2 = h2.tocoo()
34+
time1 = time.time()
35+
36+
print("quimb time: ", time1 - time0)
37+
38+
39+
def assert_equal(h1, h2):
40+
np.testing.assert_allclose(h1.row, h2.row, atol=1e-5)
41+
np.testing.assert_allclose(h1.col, h2.col, atol=1e-5)
42+
np.testing.assert_allclose(h1.data, h2.data, atol=1e-5)
43+
44+
45+
assert_equal(h1, h2)

‎examples/vqe_extra_mpo.py

+146
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,146 @@
1+
"""
2+
Demonstration of TFIM VQE with extra size in MPO formulation
3+
"""
4+
5+
import time
6+
import logging
7+
import sys
8+
import numpy as np
9+
10+
logger = logging.getLogger("tensorcircuit")
11+
logger.setLevel(logging.INFO)
12+
ch = logging.StreamHandler()
13+
ch.setLevel(logging.DEBUG)
14+
logger.addHandler(ch)
15+
16+
sys.setrecursionlimit(10000)
17+
18+
import tensorflow as tf
19+
import tensornetwork as tn
20+
import cotengra as ctg
21+
22+
import optax
23+
import tensorcircuit as tc
24+
25+
opt = ctg.ReusableHyperOptimizer(
26+
methods=["greedy", "kahypar"],
27+
parallel="ray",
28+
minimize="combo",
29+
max_time=360,
30+
max_repeats=4096,
31+
progbar=True,
32+
)
33+
34+
35+
def opt_reconf(inputs, output, size, **kws):
36+
tree = opt.search(inputs, output, size)
37+
tree_r = tree.subtree_reconfigure_forest(
38+
parallel="ray",
39+
progbar=True,
40+
num_trees=20,
41+
num_restarts=20,
42+
subtree_weight_what=("size",),
43+
)
44+
return tree_r.path()
45+
46+
47+
tc.set_contractor("custom", optimizer=opt_reconf, preprocessing=True)
48+
tc.set_dtype("complex64")
49+
tc.set_backend("tensorflow")
50+
# jax backend is incompatible with keras.save
51+
52+
dtype = np.complex64
53+
54+
nwires, nlayers = 150, 7
55+
56+
57+
Jx = np.array([1.0 for _ in range(nwires - 1)]) # strength of xx interaction (OBC)
58+
Bz = np.array([-1.0 for _ in range(nwires)]) # strength of transverse field
59+
hamiltonian_mpo = tn.matrixproductstates.mpo.FiniteTFI(
60+
Jx, Bz, dtype=dtype
61+
) # matrix product operator
62+
hamiltonian_mpo = tc.quantum.tn2qop(hamiltonian_mpo)
63+
64+
65+
def vqe_forward(param):
66+
print("compiling")
67+
split_conf = {
68+
"max_singular_values": 2,
69+
"fixed_choice": 1,
70+
}
71+
c = tc.Circuit(nwires, split=split_conf)
72+
for i in range(nwires):
73+
c.H(i)
74+
for j in range(nlayers):
75+
for i in range(0, nwires - 1):
76+
c.exp1(
77+
i,
78+
(i + 1) % nwires,
79+
theta=param[4 * j, i],
80+
unitary=tc.gates._xx_matrix,
81+
)
82+
83+
for i in range(nwires):
84+
c.rz(i, theta=param[4 * j + 1, i])
85+
for i in range(nwires):
86+
c.ry(i, theta=param[4 * j + 2, i])
87+
for i in range(nwires):
88+
c.rz(i, theta=param[4 * j + 3, i])
89+
return tc.templates.measurements.mpo_expectation(c, hamiltonian_mpo)
90+
91+
92+
if __name__ == "__main__":
93+
refresh = False
94+
95+
time0 = time.time()
96+
if refresh:
97+
tc_vag = tf.function(
98+
tc.backend.value_and_grad(vqe_forward),
99+
input_signature=[tf.TensorSpec([4 * nlayers, nwires], tf.float32)],
100+
)
101+
tc.keras.save_func(tc_vag, "./funcs/%s_%s_tfim_mpo" % (nwires, nlayers))
102+
time1 = time.time()
103+
print("staging time: ", time1 - time0)
104+
105+
tc_vag_loaded = tc.keras.load_func("./funcs/%s_%s_tfim_mpo" % (nwires, nlayers))
106+
107+
lr1 = 0.008
108+
lr2 = 0.06
109+
steps = 2000
110+
switch = 400
111+
debug_steps = 20
112+
113+
if tc.backend.name == "jax":
114+
opt = tc.backend.optimizer(optax.adam(lr1))
115+
opt2 = tc.backend.optimizer(optax.sgd(lr2))
116+
else:
117+
opt = tc.backend.optimizer(tf.keras.optimizers.Adam(lr1))
118+
opt2 = tc.backend.optimizer(tf.keras.optimizers.SGD(lr2))
119+
120+
times = []
121+
param = tc.backend.implicit_randn(stddev=0.1, shape=[4 * nlayers, nwires])
122+
123+
for j in range(steps):
124+
loss, gr = tc_vag_loaded(param)
125+
if j < switch:
126+
param = opt.update(gr, param)
127+
else:
128+
if j == switch:
129+
print("switching the optimizer")
130+
param = opt2.update(gr, param)
131+
if j % debug_steps == 0 or j == steps - 1:
132+
times.append(time.time())
133+
print("loss", tc.backend.numpy(loss))
134+
if j > 0:
135+
print("running time:", (times[-1] - times[0]) / j)
136+
137+
138+
"""
139+
# Baseline code: obtained from DMRG using quimb
140+
141+
import quimb
142+
143+
h = quimb.tensor.tensor_gen.MPO_ham_ising(nwires, 4, 2, cyclic=False)
144+
dmrg = quimb.tensor.tensor_dmrg.DMRG2(m, bond_dims=[10, 20, 100, 100, 200], cutoffs=1e-10)
145+
dmrg.solve(tol=1e-9, verbosity=1) # may require repetition of this API
146+
"""

‎examples/vqeh2o_benchmark.py

+153
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,153 @@
1+
"""
2+
Time comparison for different evaluation approach on molecule VQE
3+
"""
4+
import time
5+
import numpy as np
6+
from openfermion.chem import MolecularData
7+
from openfermion.transforms import (
8+
get_fermion_operator,
9+
binary_code_transform,
10+
checksum_code,
11+
reorder,
12+
)
13+
from openfermion.chem import geometry_from_pubchem
14+
from openfermion.utils import up_then_down
15+
from openfermionpyscf import run_pyscf
16+
17+
import tensorcircuit as tc
18+
19+
K = tc.set_backend("tensorflow")
20+
21+
22+
n = 12
23+
nlayers = 2
24+
multiplicity = 1
25+
basis = "sto-3g"
26+
# 14 spin orbitals for H2O
27+
geometry = geometry_from_pubchem("h2o")
28+
description = "h2o"
29+
molecule = MolecularData(geometry, basis, multiplicity, description=description)
30+
molecule = run_pyscf(molecule, run_mp2=True, run_cisd=True, run_ccsd=True, run_fci=True)
31+
mh = molecule.get_molecular_hamiltonian()
32+
fh = get_fermion_operator(mh)
33+
b = binary_code_transform(reorder(fh, up_then_down), 2 * checksum_code(7, 1))
34+
lsb, wb = tc.templates.chems.get_ps(b, 12)
35+
print("%s terms in H2O qubit Hamiltonian" % len(wb))
36+
mb = tc.quantum.PauliStringSum2COO_numpy(lsb, wb)
37+
mbd = mb.todense()
38+
mb = K.coo_sparse_matrix(
39+
np.transpose(np.stack([mb.row, mb.col])), mb.data, shape=(2 ** n, 2 ** n)
40+
)
41+
mbd = tc.array_to_tensor(mbd)
42+
43+
44+
def ansatz(param):
45+
c = tc.Circuit(n)
46+
for i in [0, 1, 2, 3, 4, 6, 7, 8, 9, 10]:
47+
c.X(i)
48+
for j in range(nlayers):
49+
for i in range(n - 1):
50+
c.cz(i, i + 1)
51+
for i in range(n):
52+
c.rx(i, theta=param[j, i])
53+
return c
54+
55+
56+
def benchmark(vqef, tries=3):
57+
if tries < 0:
58+
vagf = K.value_and_grad(vqef)
59+
else:
60+
vagf = K.jit(K.value_and_grad(vqef))
61+
time0 = time.time()
62+
v, g = vagf(K.zeros([nlayers, n]))
63+
time1 = time.time()
64+
for _ in range(tries):
65+
v, g = vagf(K.zeros([nlayers, n]))
66+
time2 = time.time()
67+
print("staging time: ", time1 - time0)
68+
if tries > 0:
69+
print("running time: ", (time2 - time1) / tries)
70+
return (v, g), (time1 - time0, (time2 - time1) / tries)
71+
72+
73+
# 1. plain Pauli sum
74+
75+
76+
def vqe1(param):
77+
c = ansatz(param)
78+
loss = 0.0
79+
for ps, w in zip(lsb, wb):
80+
obs = []
81+
for i, p in enumerate(ps):
82+
if p == 1:
83+
obs.append([tc.gates.x(), [i]])
84+
elif p == 2:
85+
obs.append([tc.gates.y(), [i]])
86+
elif p == 3:
87+
obs.append([tc.gates.z(), [i]])
88+
89+
loss += w * c.expectation(*obs)
90+
return K.real(loss)
91+
92+
93+
# 2. vmap the Pauli sum
94+
95+
96+
def measurement(s, structure):
97+
c = tc.Circuit(n, inputs=s)
98+
return tc.templates.measurements.parameterized_measurements(
99+
c, structure, onehot=True
100+
)
101+
102+
103+
measurement = K.jit(K.vmap(measurement, vectorized_argnums=1))
104+
structures = tc.array_to_tensor(lsb)
105+
weights = tc.array_to_tensor(wb, dtype="float32")
106+
107+
108+
def vqe2(param):
109+
c = ansatz(param)
110+
s = c.state()
111+
ms = measurement(s, structures)
112+
return K.sum(ms * weights)
113+
114+
115+
# 3. dense matrix
116+
117+
118+
def vqe3(param):
119+
c = ansatz(param)
120+
return tc.templates.measurements.operator_expectation(c, mbd)
121+
122+
123+
# 4. sparse matrix
124+
125+
126+
def vqe4(param):
127+
c = ansatz(param)
128+
return tc.templates.measurements.operator_expectation(c, mb)
129+
130+
131+
# 5. mpo (ommited, since it is not that applicable for molecule/long range Hamiltonian
132+
# due to large bond dimension)
133+
134+
135+
if __name__ == "__main__":
136+
r0 = None
137+
des = [
138+
"plain Pauli sum",
139+
"vmap Pauli sum",
140+
"dense Hamiltonian matrix",
141+
"sparse Hamiltonian matrix",
142+
]
143+
vqef = [vqe1, vqe2, vqe3, vqe4]
144+
tries = [-1, 2, 5, 5]
145+
for i in range(4):
146+
print(des[i])
147+
r1, _ = benchmark(vqef[i], tries=tries[i])
148+
# plain approach takes too long to jit
149+
if r0 is not None:
150+
np.testing.assert_allclose(r0[0], r1[0], atol=1e-5)
151+
np.testing.assert_allclose(r0[1], r1[1], atol=1e-5)
152+
r0 = r1
153+
print("------------------")

0 commit comments

Comments
 (0)