Closed millasub closed 5 months ago
Hello,
Can you provide a minimal reproducer to get the results from both our MPSHelper
and qiskit? I'm surprised that bondDim=50
is suffice to converge an MPS with 100 qubits and depth 100-700 unless the state is very sparse. And how are the "true" values obtained in your case?
For sure. Here is an example (with 20 qubits):
import cupy as cp
import numpy as np
import scipy as sp
from qiskit import execute
from qiskit.circuit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit_aer.primitives import Estimator as AerEstimator
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.providers.aer import *
from cuquantum import contract, contract_path, CircuitToEinsum, tensor
from cuquantum import cutensornet as cutn
from cuquantum.cutensornet.experimental import contract_decompose
np.random.seed(0)
cp.random.seed(0)
def Hpauli(dict,N):
tableop = ""
tableind = []
for index, Pauli in sorted(dict.items(),reverse=True):
tableop+=Pauli
tableind.append(index)
operator = SparsePauliOp.from_sparse_list([(tableop, tableind, 1)], num_qubits = N)
return operator.simplify()
def Hm(N):
op = SparsePauliOp.from_sparse_list([("I", [0], 0)], num_qubits = N)
for n in range(N):
op = op + (-1)**(n)/2 * Hpauli({n: "Z"},N)
return op.simplify()
def Hkin(N):
op = SparsePauliOp.from_sparse_list([("I", [0], 0)], num_qubits = N)
for n in range(N-1):
op = op + 1/2 * (Hpauli({n: "X", n+1: "X"},N) + Hpauli({n: "Y", n+1: "Y"},N))
return (1/2 * op).simplify()
def Hint(N):
op = SparsePauliOp.from_sparse_list([("I", [0], 0)], num_qubits = N)
for n in range(int(N/2)):
op = op + (1/4+1/8*(-1)**n)*Hpauli({n: "Z"},N)
for n in range(N-1,int(N/2)-1,-1):
op = op - (1/4+1/8*(-1)**(n+1))*Hpauli({n: "Z"},N)
for n in range(N-1):
op = op + abs(int(N/2)-1-n)/8*Hpauli({n: "Z", n+1: "Z"},N)
return op.simplify()
def SecondTrotter(N,ntrot,t):
circ = QuantumCircuit(N)
for i in range(int(N/2)):
circ.x(2*i)
hmass = Hm(N)
hkin = Hkin(N)
hint = Hint(N)
for _ in range(ntrot):
Ukin = PauliEvolutionGate(hkin,time=(t/2/ntrot))
circ.append(Ukin, range(N))
Umass = PauliEvolutionGate(hmass,time=0.5*(t/ntrot))
circ.append(Umass, range(N))
Uint = PauliEvolutionGate(hint,time=0.3**2*(t/ntrot))
circ.append(Uint, range(N))
Ukin = PauliEvolutionGate(hkin[::-1],time=(t/2/ntrot))
circ.append(Ukin, range(N))
return circ.decompose().decompose().decompose()
def get_initial_mps(num_qubits, dtype='complex128'):
state_tensor = cp.asarray([1, 0], dtype=dtype).reshape(1,2,1)
mps_tensors = [state_tensor] * num_qubits
return mps_tensors
def mps_site_right_swap(mps_tensors,i,algorithm=None,options=None):
a, _, b = contract_decompose('ipj,jqk->iqj,jpk', *mps_tensors[i:i+2], algorithm=algorithm, options=options)
mps_tensors[i:i+2] = (a, b)
return mps_tensors
def apply_gate(mps_tensors,gate,qubits,algorithm=None,options=None):
n_qubits = len(qubits)
if n_qubits == 1:
# single-qubit gate
i = qubits[0]
mps_tensors[i] = contract('ipj,qp->iqj', mps_tensors[i], gate, options=options) # in-place update
elif n_qubits == 2:
# two-qubit gate
i, j = qubits
if i > j:
# swap qubits order
return apply_gate(mps_tensors, gate.transpose(1,0,3,2), (j, i), algorithm=algorithm, options=options)
elif i+1 == j:
# two adjacent qubits
a, _, b = contract_decompose('ipj,jqk,rspq->irj,jsk', *mps_tensors[i:i+2], gate, algorithm=algorithm, options=options)
mps_tensors[i:i+2] = (a, b) # in-place update
else:
# non-adjacent two-qubit gate
# step 1: swap i with i+1
mps_site_right_swap(mps_tensors, i, algorithm=algorithm, options=options)
# step 2: apply gate to (i+1, j) pair. This amounts to a recursive swap until the two qubits are adjacent
apply_gate(mps_tensors, gate, (i+1, j), algorithm=algorithm, options=options)
# step 3: swap back i and i+1
mps_site_right_swap(mps_tensors, i, algorithm=algorithm, options=options)
else:
raise NotImplementedError("Only one- and two-qubit gates supported")
return mps_tensors
class MPSContractionHelper:
def __init__(self, num_qubits):
self.num_qubits = num_qubits
self.path_cache = dict()
self.bra_modes = [(2*i, 2*i+1, 2*i+2) for i in range(num_qubits)]
offset = 2*num_qubits+1
self.ket_modes = [(i+offset, 2*i+1, i+1+offset) for i in range(num_qubits)]
def contract_norm(self, mps_tensors, options=None):
interleaved_inputs = []
for i, o in enumerate(mps_tensors):
interleaved_inputs.extend([o, self.bra_modes[i], o.conj(), self.ket_modes[i]])
interleaved_inputs.append([]) # output
return self._contract('norm', interleaved_inputs, options=options).real
def contract_state_vector(self, mps_tensors, options=None):
interleaved_inputs = []
for i, o in enumerate(mps_tensors):
interleaved_inputs.extend([o, self.bra_modes[i]])
output_modes = tuple([bra_modes[1] for bra_modes in self.bra_modes])
interleaved_inputs.append(output_modes) # output
return self._contract('sv', interleaved_inputs, options=options)
def contract_expectation(self, mps_tensors, operator, qubits, normalize=False, options=None):
interleaved_inputs = []
extra_mode = 3 * self.num_qubits + 2
operator_modes = [None] * len(qubits) + [self.bra_modes[q][1] for q in qubits]
qubits = list(qubits)
for i, o in enumerate(mps_tensors):
interleaved_inputs.extend([o, self.bra_modes[i]])
k_modes = self.ket_modes[i]
if i in qubits:
k_modes = (k_modes[0], extra_mode, k_modes[2])
q = qubits.index(i)
operator_modes[q] = extra_mode # output modes
extra_mode += 1
interleaved_inputs.extend([o.conj(), k_modes])
interleaved_inputs.extend([operator, tuple(operator_modes)])
interleaved_inputs.append([]) # output
if normalize:
norm = self.contract_norm(mps_tensors, options=options)
else:
norm = 1
return self._contract(f'exp{qubits}', interleaved_inputs, options=options) / norm
def contract_mps_mpo_to_state_vector(self, mps_tensors, mpo_tensors, options=None):
interleaved_inputs = []
for i, o in enumerate(mps_tensors):
interleaved_inputs.extend([o, self.bra_modes[i]])
output_modes = []
offset = 2 * self.num_qubits + 1
for i, o in enumerate(mpo_tensors):
mpo_modes = (2*i+offset, 2*i+offset+1, 2*i+offset+2, 2*i+1)
output_modes.append(2*i+offset+1)
interleaved_inputs.extend([o, mpo_modes])
interleaved_inputs.append(output_modes)
return self._contract('mps_mpo', interleaved_inputs, options=options)
def _contract(self, key, interleaved_inputs, options=None):
if key not in self.path_cache:
self.path_cache[key] = contract_path(*interleaved_inputs, options=options)[0]
path = self.path_cache[key]
return contract(*interleaved_inputs, options=options, optimize={'path':path})
num_qubits = 20
t = 6
ntrot = 6
bondDim = 10
circ = SecondTrotter(num_qubits,ntrot,t)
# cuQuantum MPS simulator
handle = cutn.create()
options = {'handle': handle}
exact_gate_algorithm = {'qr_method': False, 'svd_method':{'partition': 'UV', 'abs_cutoff':1e-12, 'max_extent': bondDim}}
dtype = 'complex128'
mps_helper = MPSContractionHelper(num_qubits)
mps_tensors = get_initial_mps(num_qubits, dtype=dtype)
myconverter = CircuitToEinsum(circ, dtype=dtype, backend=cp)
gates = myconverter.gates
gate_map = dict(zip(myconverter.qubits, range(num_qubits)))
for (gate, qubits) in gates:
qubits = [gate_map[q] for q in qubits]
apply_gate(mps_tensors, gate, qubits, algorithm=exact_gate_algorithm, options=options)
zvalC = []
for i in range(num_qubits):
target_qubits = (i,)
gate_z = cp.asarray([[1, 0], [0, -1]]).astype(dtype)
expec_mps_reference = mps_helper.contract_expectation(mps_tensors, gate_z, target_qubits, options=options, normalize=True)
zvalC.append(expec_mps_reference.real)
print(cp.asnumpy(cp.array(zvalC).flatten()))
cutn.destroy(handle)
# Qiskit-Aer MPS simulator
estimator = AerEstimator(run_options= {"method": "matrix_product_state", "shots": None, "matrix_product_state_max_bond_dimension": bondDim}, approximation=True)
zvalQ = []
for i in range(num_qubits):
op = Hpauli({i: "Z"},num_qubits)
expectation_value = estimator.run(circ, op).result().values
zvalQ.append(expectation_value[0])
print(zvalQ)
# Qiskit-Aer exact statevector
backend = AerSimulator(method='statevector')
circ.save_statevector()
job = execute(circ, backend)
svcirc = job.result().get_statevector(circ)
exact = []
for i in range(num_qubits):
exact.append(np.real(svcirc.expectation_value(Hpauli({i: "Z"},num_qubits))))
print(exact)
And here are the results for different circuit depths (trotter steps) and bond dimensions:
A beginner starting to look into cuQuantum here. From my naive glance, I suspect that the "exact_gate_algorithm" is the part that goes wrong. It is the example code to check when MPS can exactly represent the states when no truncation occurs.
In your case, truncations are necessary and to have good truncation, TEBD-like of apply_gate is required, i.e. utilizing canonical form to make optimal local truncation. I believe the qiskit simulator must have done this correctly.
I was playing with the MPS simulator from https://github.com/NVIDIA/cuQuantum/blob/main/python/samples/cutensornet/tn_algorithms/mps_algorithms.ipynb, specifically trying to test its performance when reducing the bond dimension when simulating some Trotterized time evolution circuits. Following that notebook, I transform my qiskit circuit to tensor form with
CircuitToEinsum
, and then apply the gates with theapply_gate
function (modifying theexact_gate_algorithm
asexact_gate_algorithm = {'qr_method': False, 'svd_method':{'partition': 'UV', 'abs_cutoff':1e-12, 'max_extent': bondDim}}
, withbondDim
being the parameter that controls the bond dimension), and finally compute some observable withmps_helper.contract_expectation
. What I'm seeing is that if I compare with the MPS simulator from Qiskit Aer, for the same value of bond dimension (controlled asAerEstimator(run_options= {"method": "matrix_product_state", "shots": None, "matrix_product_state_max_bond_dimension": bondDim}, approximation=True)
), the cuQuantum result is further away from the "true" value than the qiskit-aer one (for circuits with around 100 qubits and depth between 100 and 700, withbondDim=50
, qiskit results already look they have converged, but not the cuQuantum ones). Am I doing something wrong on the cuQuantum side, or I cannot compare the results from the two methods? (I am usingcuquantum-appliance:23.10
) Thank you!