# Why am I able to simulate such high qubit numbers on my laptop?

Recently I've been working on a VQE-related project, for which I'm using Qiskit's TwoLocal variational form (qiskit version 0.34.2). I noticed that I was able to simulate this circuit up to surprisingly high qubit numbers on my laptop, so I wanted to investigate this in a bit more detail. I ran the following piece of code for varying num_qubits:

from qiskit import Aer, execute
from qiskit.circuit.library import TwoLocal
import numpy as np
from time import time

qc = TwoLocal(num_qubits, 'ry', 'cx', reps=1, entanglement='linear')
params = np.random.rand(len(qc.parameters))
qc = qc.bind_parameters(params)
qc.measure_all()
t1 = time()
results = execute(qc, backend=Aer.get_backend('qasm_simulator'), shots=1024).result().get_counts()
t2 = time()


When I plot t2-t1 as a function of num_qubits, I get the following result: Note that linear and full refer to the corresponding entanglement modes of the TwoLocal function for reps=1. 'none' refers to the TwoLocal circuit with reps=0. In those regions without data points my code actually crashed, i.e. when running inside a jupyter notebook it says "kernel died" and when running it inside a script in the terminal it says "Segmentation fault (core dumped)". Also, there's a special mention for the case of 72 qubits, for which the code doesn't crash but the result always gives $$|0\rangle^{\otimes n}$$ with 100% probability, which obviously doesn't make sense.

My main question is the following: why am I able to simulate entangled circuits with such high qubit numbers, and quite fast at that? According to this link the only simulator which simulates more than 100 qubits is the stabilizer simulator, but I'm using gates which are not in its supported gates set (i.e. Ry gates with random angles).

I have two related bonus questions (I know each post should ideally be "laser-focused" but I think these questions are too related to warrant a separate post so please bear with me):

1. What on earth is going on in those regions where the code crashes?
2. What simulators do you think are being used in the different regions? For example, qiskit is clearly changing simulators when going from 29 qubits to 30 qubits but I'm wondering what the downside of the faster simulator is, since if there wouldn't be any then why not just always use that one?

Thanks in advance to anyone who can point me in the right direction!

• Your circuit is extremely shallow, so it can be simulated by tensor network methods. Not sure if that's actually what qsikit is doing, and also I'm not exactly sure how it's decomposing the full connectivity. But with single-layer linear it's trivial to iterate over the qubits instead of over time while using constant memory. Mar 23 at 17:35
• Hmm okay, thanks. That makes me question the practical use of VQE for combinatorial optimization problems even more, in addition to the fact that there are indications that there is no advantage in using entanglement layers in the TwoLocal variational form for these problems. Mar 25 at 14:47

The AerSimulator backend supports a variety of simulation methods, each of which supports a different set of instructions and has a different memory requirements. The method can be set manually using simulator.set_option(method=value) .

The default simulation method is automatic which will automatically select one of the other simulation methods for each circuit based on the instructions in the circuit, noise model, and available memory.

According to the source code here (C++ code), selecting the simulation method follows these steps:

1. If circuit and noise model are Clifford, then stabilizer simulator will be used.

2. For noisy simulations and number of qubits < 64 the density_matrix method is used if shots > $$2^{num\_qubits}$$.

3. Otherwise, simulation method is chosen based on the supported operations only with preference given by memory requirements: statevector > density_matrix > matrix_product_state > unitary > superop. typically any save state instructions will decide the method.

• Circuit is not Clifford -> stabilizer method is not selected.
• Number of qubits > 64 -> density_matrix method is not selected.
• Memory is not sufficient -> statevector methods is not selected.
So, matrix_product_state is chosen.
• Thanks for the clear answer and the useful links. I've been playing around with the different simulators and actually it seems that Qiskit is chosing the density_matrix method (when I fix this method I can recreate the above (strange) results). However, when I fix matrix_product_state as the method it also works and it even fixes the issue of my code crashing for specific qubit numbers and the wrong result I obtained before for 72 qubits. Mar 24 at 18:09
• I was also confused by the fact that this link says that the matrix_product_state has a maximum qubit number of 100, but I guess this limit doesn't apply when simulating on a local machine? Mar 24 at 18:11
• Edit to my first comment (not enough characters left unfortunately): Actually the density_matrix method is giving me an "insufficient memory" error from 15 qubits onwards. However, in those regions in my plot where there are no data points it doesn't give this error but instead just crashes my code (again except for the specific case of 72 qubits where it gives wrong results). matrix_product_state always works so I'll just be using that one Mar 24 at 18:23