cuQuantum MPS Simulator vs Qiskit Aer about cuquantum

[yangcal]

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?

[millasub]

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:

[ShHsLin]

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.