What explains my anomalously scaled up VQE?

Question

I am trying to implement VQE from the Qiskit to obtain the ground state of a very specific Hamiltonian that has been generated via a docplex minimized quadratic model. The model has been converted to an Ising Hamiltonian using Qiskit's Optimization module. The resultant Hamiltonian denoted by H is as follows:

from qiskit.providers.aer import AerSimulator, QasmSimulator
from qiskit.algorithms import VQE
from qiskit.algorithms.optimizers import COBYLA 
from qiskit.circuit.library import TwoLocal
from qiskit import *
from qiskit.opflow import OperatorBase
from qiskit.opflow import Z, X, I  # Pauli Z, X matrices and identity
import pylab
import matplotlib.pyplot as plt
import numpy as np

H =   504.0 * I^I^I^I^I^I^I^Z+1008.0 * I^I^I^I^I^I^Z^I+2016.0 * I^I^I^I^I^Z^I^I+504.0 * I^I^I^I^Z^I^I^I+1143.7999999999997 * I^I^I^Z^I^I^I^I+2287.6 * I^I^Z^I^I^I^I^I+4575.200000000001 * I^Z^I^I^I^I^I^I+1143.7999999999997 * Z^I^I^I^I^I^I^I+98.0 * I^I^I^I^I^I^Z^Z+196.0 * I^I^I^I^I^Z^I^Z+392.0 * I^I^I^I^I^Z^Z^I+49.0 * I^I^I^I^Z^I^I^Z+98.0 * I^I^I^I^Z^I^Z^I+196.0 * I^I^I^I^Z^Z^I^I+93.1 * I^I^Z^Z^I^I^I^I+186.2 * I^Z^I^Z^I^I^I^I+372.4 * I^Z^Z^I^I^I^I^I+46.55 * Z^I^I^Z^I^I^I^I+93.1 * Z^I^Z^I^I^I^I^I+186.2 * Z^Z^I^I^I^I^I^I

backend = QasmSimulator()
optimizer = COBYLA(maxiter=2000)
ansatz = TwoLocal(num_qubits=8, rotation_blocks='ry', entanglement_blocks=None, entanglement='full', reps=1, skip_unentangled_qubits=False, skip_final_rotation_layer=False)
# set the algorithm
vqe = VQE(ansatz, optimizer, quantum_instance=backend)

#run it with the Hamiltonian we defined above
result = vqe.compute_minimum_eigenvalue(H) 

This however yields the error:

'Circuit execution failed: ERROR:  [Experiment 0] QasmSimulator: Insufficient memory for 141-qubit circuit using "statevector" method. You could try using the "matrix_product_state" or "extended_stabilizer" method instead.'

My questions are:

1. How and why does my circuit yield 141 qubits when there are only 8 Pauli Operators in each term of my Hamiltonian? What am I missing conceptually?
2. How do we calculate the number of qubits required when solving this sort of problem?

Answer 1

Easy Fix:

It seems like it is because of the way you define $H$. You need the parenthesis around each of the term!

So something like:

H =    (504.0 * I^I^I^I^I^I^I^Z) + (1008.0 * I^I^I^I^I^I^Z^I) + ( 2016.0 *  I^I^I^I^I^Z^I^I)

Just replace this in your code then it will work!

---

Alternative (longer) way:

Here I will offer another way to define the Hamiltonian in case you are curious. Instead of doing the above, if you replace it with something like:

from qiskit.aqua.operators import *
pauli_terms = ['IIIIIIIZ', 'IIIIIIZI', 'IIIIIZII' ]
pauli_weights = [504.0, 1008.0, 2016.0]
pauli_dict = {'paulis': [{"coeff": {"imag": 0., "real": pauli_weights[i] }, "label": pauli_terms[i]} \
                         for i in range(len(pauli_terms))]}
H = WeightedPauliOperator.from_dict(pauli_dict)

This should work. I just grabbed the first 3 entries of your Hamiltonian. Here is the full script for you to reproduce the result:

import numpy as np
from qiskit.providers.aer import AerSimulator, QasmSimulator
from qiskit.algorithms.optimizers import COBYLA 
from qiskit.circuit.library import TwoLocal
from qiskit.aqua.operators import *
from qiskit.aqua import set_qiskit_aqua_logging, QuantumInstance
from qiskit.aqua.algorithms import NumPyMinimumEigensolver, VQE, NumPyEigensolver
from qiskit.circuit import QuantumCircuit, ParameterVector

pauli_terms = ['IIIIIIIZ', 'IIIIIIZI', 'IIIIIZII' ]
pauli_weights = [504.0, 1008.0, 2016.0]
pauli_dict = {'paulis': [{"coeff": {"imag": 0., "real": pauli_weights[i] }, "label": pauli_terms[i]} \
                         for i in range(len(pauli_terms))]}
Hamiltonian = WeightedPauliOperator.from_dict(pauli_dict)
ansatz = TwoLocal(num_qubits=8, rotation_blocks='ry', entanglement_blocks=None, entanglement='full', reps=1, skip_unentangled_qubits=False, skip_final_rotation_layer=False)
print(ansatz)
backend = QasmSimulator()
quantum_instance = QuantumInstance(backend,
                                   shots = 8192,
                                   initial_layout = None,
                                   optimization_level = 3)  

optimizer = COBYLA(maxiter= 100, tol=0.000000001)   
vqe = VQE(Hamiltonian, ansatz, optimizer, initial_point= None , include_custom = False)
print('We are using:', quantum_instance.backend)
vqe_result = vqe.run(quantum_instance)
vqe_result['eigenvalue']

output:

     ┌──────────┐ ┌──────────┐
q_0: ┤ RY(θ[0]) ├─┤ RY(θ[8]) ├
     ├──────────┤ ├──────────┤
q_1: ┤ RY(θ[1]) ├─┤ RY(θ[9]) ├
     ├──────────┤┌┴──────────┤
q_2: ┤ RY(θ[2]) ├┤ RY(θ[10]) ├
     ├──────────┤├───────────┤
q_3: ┤ RY(θ[3]) ├┤ RY(θ[11]) ├
     ├──────────┤├───────────┤
q_4: ┤ RY(θ[4]) ├┤ RY(θ[12]) ├
     ├──────────┤├───────────┤
q_5: ┤ RY(θ[5]) ├┤ RY(θ[13]) ├
     ├──────────┤├───────────┤
q_6: ┤ RY(θ[6]) ├┤ RY(θ[14]) ├
     ├──────────┤├───────────┤
q_7: ┤ RY(θ[7]) ├┤ RY(θ[15]) ├
     └──────────┘└───────────┘
We are using: qasm_simulator

(-3528+0j)

Although I am not sure why you would have two consecutive layers of $RY$ rotation though.