cuQuantum MPS Simulator vs Qiskit Aer

[millasub]

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 the apply_gate function (modifying the exact_gate_algorithm as exact_gate_algorithm = {'qr_method': False, 'svd_method':{'partition': 'UV', 'abs_cutoff':1e-12, 'max_extent': bondDim}}, with bondDim being the parameter that controls the bond dimension), and finally compute some observable with mps_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 as AerEstimator(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, with bondDim=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 using cuquantum-appliance:23.10)
Thank you!

[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.

[yangcal]

Hello,

Thanks for sharing the script and thanks @ShHsLin for helpful inputs on the topic.

I do want to first clarify that the notebooks in our sample directory are only meant as basic tutorials to help users understand how to integrate cuquantum functionalities into their own work stack. They're not meant as research quality showcases (they can be made to for sure).

Back to the accuracy observation, like @ShHsLin pointed out, there are different algorithms to actually perform MPS gating, for instance, applying canonicalization before contract/decompose generally improves the accuracy. We're not sure what qiskit does in their implementation, but the demo notebook performs a basic direct contract-decompose operation. Therefore the lower accuracy is not surprising. One simple way to improve the MPS accuracy without full canonicalization is to maintain gauges (basically singular values) for the MPS (some people call this method simple update) and the script for doing that is attached at the end of this script (pay attention to the gauges and use_gauge variable).

On my end with t/ntrot = 6/10, I'm seeing improved numbers using gauges, and they're also about the same as qiskit's numbers (so maybe they did this too instead of full canonicalization).

In any case, while we're looking into providing more MPS features directly in our public API, you're always welcome to play around with these demos or build on top of it and let us know if you have more questions.

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_tensor_absorb_gauge(t, s, direction='left'):
    if s is None:
        return t
    assert direction in {'left', 'right'}
    subscripts = {'left': 'i', 'right': 'k'}[direction] + ',ijk->ijk'
    return cp.einsum(subscripts, s, t)

def mps_site_right_swap(mps_tensors,i, gauges=None, algorithm=None, options=None):
    if gauges is not None:
        s = gauges.get(i+1, None)
        mps_tensors[i] = mps_tensor_absorb_gauge(mps_tensors[i], s, direction='right')
    a, s, b = contract_decompose('ipj,jqk->iqj,jpk', *mps_tensors[i:i+2], algorithm=algorithm, options=options)
    mps_tensors[i:i+2] = (a, b)
    if s is not None:
        gauges[i+1] = s
    return mps_tensors

def apply_gate(mps_tensors, gate, qubits, gauges=None, use_gauges=False, algorithm=None, options=None):
    n_qubits = len(qubits)
    if gauges is None:
        gauges = dict()
    if use_gauges:
        assert algorithm['svd_method']['partition'] is None, "For MPS with gauges, SVD partition must be set to None"
    else:
        assert algorithm['svd_method']['partition'] == 'UV', "For MPS without gauges, SVD partition must be set to UV"
    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), gauges=gauges, use_gauges=use_gauges, algorithm=algorithm, options=options)
        elif i+1 == j:
            # two adjacent qubits
            sa = gauges.get(i, None) # left gauge of i   
            s = gauges.get(i+1, None) # shared gauge of i, j
            sb = gauges.get(j+1, None) # right gauge of j

            # absorb all gauges before contract & decompose
            mps_tensors[i] = mps_tensor_absorb_gauge(mps_tensors[i], sa, direction="left")
            mps_tensors[i] = mps_tensor_absorb_gauge(mps_tensors[i], s, direction="right")
            mps_tensors[j] = mps_tensor_absorb_gauge(mps_tensors[j], sb, direction="right")
            
            a, s, b = contract_decompose('ipj,jqk,rspq->irj,jsk', *mps_tensors[i:i+2], gate, algorithm=algorithm, options=options)
            # update shared gauge
            if s is not None:
                gauges[i+1] = s
            # remove gauge effect back
            mps_tensors[i] = mps_tensor_absorb_gauge(a, None if sa is None else 1/sa, direction="left")
            mps_tensors[j] = mps_tensor_absorb_gauge(b, None if sb is None else 1/sb, direction="right")
        else:
            # non-adjacent two-qubit gate
            # step 1: swap i with i+1
            mps_site_right_swap(mps_tensors, i, gauges=gauges, 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), gauges=gauges, use_gauges=use_gauges, algorithm=algorithm, options=options)
            # step 3: swap back i and i+1
            mps_site_right_swap(mps_tensors, i, gauges=gauges, algorithm=algorithm, options=options)
    else:
        raise NotImplementedError("Only one- and two-qubit gates supported")
    if use_gauges:
        return mps_tensors, gauges
    else:
        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 absorb_gauges(self, mps_tensors, gauges=None):
        if gauges is None:
            return mps_tensors
        for i, o in enumerate(mps_tensors):
            s = gauges.get(i, None)
            mps_tensors[i] = mps_tensor_absorb_gauge(o, s, direction='left')
            
    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, gauges=None, options=None):
        self.absorb_gauges(mps_tensors, gauges=gauges)
        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 = 10
ntrot = 10
bondDim = 10
circ = SecondTrotter(num_qubits,ntrot,t)

# cuQuantum MPS simulator

handle = cutn.create()
options = {'handle': handle}

# directly perform contract and decompose to perform the gate operation
basic_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=basic_gate_algorithm, 
               use_gauges=False, 
               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("Basic MPS without gauging")
print(cp.asnumpy(cp.array(zvalC).flatten()))

# use gauged gate method to perform another MPS simulation
mps_tensors = get_initial_mps(num_qubits, dtype=dtype)
gauges = dict()
gauged_gate_algorithm = {'qr_method': False, 
                         'svd_method':{'partition': None, 'abs_cutoff':1e-12, 'max_extent': bondDim}}
for (gate, qubits) in gates:
    qubits = [gate_map[q] for q in qubits]
    mps_tensors, gauges = apply_gate(mps_tensors, 
                                     gate, 
                                     qubits, 
                                     algorithm=gauged_gate_algorithm, 
                                     gauges=gauges, 
                                     use_gauges=True, 
                                     options=options)
# absorb gauges to the MPS state
mps_helper.absorb_gauges(mps_tensors, gauges=gauges)

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("MPS with gauging")
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)

[sunggonkim]

We are running the exactly the same notebook.

The error occurs when execute this part of the code

for (gate, qubits) in gates:
    qubits = [gate_map[q] for q in qubits]
    apply_gate(mps_tensors, 
               gate, 
               qubits, 
               algorithm=basic_gate_algorithm, 
               use_gauges=False, 
               options=options)

In here, the apply_gate is called multiple times through for loop.
While first few iterations works without any error, the floating point error occurs after the few iteration.

Following is the apply_gate code (identical to the notebook but with just some print)

   n_qubits = len(qubits)
    if n_qubits == 1:
        print("SKIM: MPS single qubit gate")
        # 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:
        print("SKIM: MPS 2 qubit gate")
        # two-qubit gate
        i, j = qubits
        if i > j:
            # swap qubits order
            print("SKIM: 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
            print("two adjacent qubits")
            a, _, b = contract_decompose('ipj,jqk,rspq->irj,jsk', *mps_tensors[i:i+2], gate, algorithm=algorithm, options=options)
            #floting point error after few iterations
            print("after contract_decompose")
            mps_tensors[i:i+2] = (a, b) # in-place update
        else:
            print("non-adjacent")
            # 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)

And here is the output running the notebook

SKIM: MPS single qubit gate
SKIM: MPS single qubit gate
…
SKIM: MPS single qubit gate
SKIM: MPS single qubit gate
two adjacent qubits
after contract_decompose
SKIM: MPS 2 qubit gate
two adjacent qubits
after contract_decompose
SKIM: MPS 2 qubit gate
two adjacent qubits
after contract_decompose
SKIM: MPS 2 qubit gate
two adjacent qubits
after contract_decompose
SKIM: MPS single qubit gate
SKIM: MPS single qubit gate
…
SKIM: MPS single qubit gate
SKIM: MPS 2 qubit gate
two adjacent qubits
Floating point exception

We are running this on perlmutter system in Berkeley lab.

[yangcal]

Hello,

I'm not able to reproduce this issue locally with 24.03. Can you make the following changes to the code right before contract_decompose call and share the data.npz file here? BTW, what is the cuda version in your environment?

print("two adjacent qubits")
import numpy as np
a, _, b = contract_decompose('ipj,jqk,rspq->irj,jsk', *mps_tensors[i:i+2], gate, algorithm=algorithm, options=options)
np.savez('data.npz', mps_left=mps_tensors[i].get(), mps_right=mps_tensors[i+1].get(), gate=gate.get(), algorithm=algorithm)

[sunggonkim]

Sure,

data.npz.txt

I added .txt at the end of file name due to the github not allowing the file with .npz extension.
Please remove the .txt.

We are using cuda 12.2.

Here is the output from nvcc --version

nvcc --version
nvcc: NVIDIA (R) Cuda compiler driver
Copyright (c) 2005-2023 NVIDIA Corporation
Built on Tue_Jun_13_19:16:58_PDT_2023
Cuda compilation tools, release 12.2, V12.2.91
Build cuda_12.2.r12.2/compiler.32965470_0

and part of echo $Path

/global/common/software/nersc/bin:/opt/nvidia/hpc_sdk/Linux_x86_64/23.9/cuda/12.2/compute-sanitizer:/opt/nvidia/hpc_sdk/Linux_x86_64/23.9/cuda/12.2/bin:/opt/nvidia/hpc_sdk/Linux_x86_64/23.9/cuda/12.2/libnvvp

[yangcal]

@sunggonkim Using the data you provided, I'm still not able to reproduce the Floating point exception.

Can you put up a new post in our dicussion or issue page where we can continue to follow up on this. I wanted to keep this post clean such that previous commenters no longer receive these notifications.

In the new post, can you share your results running the following standalone script?

import numpy as np
import cupy as cp
from cuquantum.cutensornet.experimental import contract_decompose
data = np.load('data.npz', allow_pickle=True)

a = cp.asarray(data['mps_left'])
b = cp.asarray(data['mps_right'])
c = cp.asarray(data['gate'])
algorithm = data['algorithm'].item()

u, s, v = contract_decompose('ipj,jqk,rspq->irj,jsk', a, b, c, algorithm=algorithm)

Also, what is the gpu in your system? You may check this via nvidia-smi in the command line.

[sunggonkim]

Sorry about this. I thought this was just minor bug fix.
I will post the result with a new issue.