Comments (3)
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?
from cuquantum.
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:
from cuquantum.
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.
from cuquantum.
Related Issues (20)
- Multithreaded cutn optimization issue HOT 1
- Disable slicing fails
- Releasing `qsim_mgpu` source on GitHub instead of only binaries the Docker container HOT 2
- gpu issue on qiskit - aer method , docker - cuquantum-appliance:23.10 HOT 1
- Website is down HOT 3
- Functions for arithmetic operations of two tensor networks HOT 3
- `cudaq` never giving correct result for `maxcut` QAOA problem. HOT 3
- `state_compute()` leading to kernel dying.
- Wrong sign in a single-gate-circuit statevector? HOT 3
- Sudo permission issue for cuquantum-appliance:23.10 container HOT 7
- 3XTF32 issue with the most recent cuquantum HOT 2
- Can't run pennyLane benchmarks in 23.10 cuQuantum Appliance HOT 8
- The output of apply_matrix_batched in cuQuantum Python is same as input HOT 4
- Small typo in README HOT 1
- Mid-circuit measurements cause significant slowdown HOT 3
- Unable to run cuqauntum in GH200 (grace hopper nvidia) docker arm64 HOT 3
- Support CUDA-Q kernels for CircuitToEinsum HOT 4
- Request for releasing a new version of cuQuantum Appliance HOT 4
- `CircuitToEinsum` fails for some qiskit `QuantumCircuit` HOT 3
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from cuquantum.