Note

There exist several parameters for configuring and using more advanced VQE capabilities. This tutorial will cover the parameters such as initial point, expectation and gradient.

It will also cover advanced simulator use such as using Aer with the Matrix Product State method.

[1]:

from qiskit import Aer
from qiskit.opflow import X, Z, I
from qiskit.utils import QuantumInstance, algorithm_globals
from qiskit.algorithms import VQE
from qiskit.algorithms.optimizers import SLSQP
from qiskit.circuit.library import TwoLocal


Here we will use the same operator as used in the other VQE algorithms tutorials.

[2]:

H2_op = (-1.052373245772859 * I ^ I) + \
(0.39793742484318045 * I ^ Z) + \
(-0.39793742484318045 * Z ^ I) + \
(-0.01128010425623538 * Z ^ Z) + \
(0.18093119978423156 * X ^ X)


## Initial point¶

The initial_point parameter allows the optimization to begin at the given point, where the point is a list of parameters that will configure the ansatz. By default the initial point is None which means that VQE will choose one. The choice in in this case is if the supplied ansatz has a preferred point, based on the initial state provided to it, then this will be chosen, otherwise a random initial point that fits with any bounds the ansatz has will be chosen. If an initial point is supplied it will take priority though and be used - note though it must match in length to the number of parameters in the ansatz circuit.

Why to use a initial point? One reason would be if you have guess a reasonable starting point for the problem or perhaps know have information from a prior experiment.

To demonstrate the use let’s first simply repeat the first working example from the algorithms introduction tutorial to get a solution’s optimal point.

[3]:

seed = 50
algorithm_globals.random_seed = seed
qi = QuantumInstance(Aer.get_backend('statevector_simulator'), seed_transpiler=seed, seed_simulator=seed)

ansatz = TwoLocal(rotation_blocks='ry', entanglement_blocks='cz')
slsqp = SLSQP(maxiter=1000)
vqe = VQE(ansatz, optimizer=slsqp, quantum_instance=qi)
result = vqe.compute_minimum_eigenvalue(operator=H2_op)
print(result)
optimizer_evals = result.optimizer_evals

{   'aux_operator_eigenvalues': None,
'cost_function_evals': 65,
'eigenstate': array([ 9.55387084e-05+0.j, -9.93766273e-01+0.j,  1.11483569e-01+0.j,
1.76709597e-05+0.j]),
'eigenvalue': (-1.8572750175665305+0j),
'optimal_parameters': {   ParameterVectorElement(θ[3]): 6.092947700274438,
ParameterVectorElement(θ[4]): -2.59832583298058,
ParameterVectorElement(θ[5]): 1.5683259187297158,
ParameterVectorElement(θ[6]): -4.717618240434922,
ParameterVectorElement(θ[7]): 0.3602071714990988,
ParameterVectorElement(θ[0]): 4.296520433306488,
ParameterVectorElement(θ[2]): 0.5470752974469076,
ParameterVectorElement(θ[1]): 4.426962161670812},
'optimal_point': array([ 4.29652043,  4.42696216,  0.5470753 ,  6.0929477 , -2.59832583,
1.56832592, -4.71761824,  0.36020717]),
'optimal_value': -1.8572750175665305,
'optimizer_evals': None,
'optimizer_time': 0.1315295696258545}


Now we can take the optimal_point from the above result and use it as the initial_point here.

[4]:

initial_pt = result.optimal_point

algorithm_globals.random_seed = seed
qi = QuantumInstance(Aer.get_backend('statevector_simulator'), seed_transpiler=seed, seed_simulator=seed)

ansatz = TwoLocal(rotation_blocks='ry', entanglement_blocks='cz')
slsqp = SLSQP(maxiter=1000)
vqe = VQE(ansatz, optimizer=slsqp, initial_point=initial_pt, quantum_instance=qi)
result1 = vqe.compute_minimum_eigenvalue(operator=H2_op)
print(result1)
optimizer_evals1 = result1.optimizer_evals
print()
print(f'optimizer_evals is {optimizer_evals1} with initial point versus {optimizer_evals} without it.')

{   'aux_operator_eigenvalues': None,
'cost_function_evals': 9,
'eigenstate': array([ 9.55387084e-05+0.j, -9.93766273e-01+0.j,  1.11483569e-01+0.j,
1.76709597e-05+0.j]),
'eigenvalue': (-1.8572750175665305+0j),
'optimal_parameters': {   ParameterVectorElement(θ[1]): 4.426962161670812,
ParameterVectorElement(θ[0]): 4.296520433306488,
ParameterVectorElement(θ[2]): 0.5470752974469076,
ParameterVectorElement(θ[3]): 6.092947700274438,
ParameterVectorElement(θ[4]): -2.59832583298058,
ParameterVectorElement(θ[5]): 1.5683259187297158,
ParameterVectorElement(θ[6]): -4.717618240434922,
ParameterVectorElement(θ[7]): 0.3602071714990988},
'optimal_point': array([ 4.29652043,  4.42696216,  0.5470753 ,  6.0929477 , -2.59832583,
1.56832592, -4.71761824,  0.36020717]),
'optimal_value': -1.8572750175665305,
'optimizer_evals': None,
'optimizer_time': 0.02120184898376465}

optimizer_evals is None with initial point versus None without it.


Here we see that result was arrived at much more quickly with optimizer_evals when it started from a random value when the initial point was not supplied (default of None).

Where this becomes useful for examples where we the solution to one problem can be used to for a guess for the solution to a very close similar problem. Chemistry is very good example where we change the inter-atomic distance(s) of molecule to plot a dissociation profile. When the distance changes are small we expect the solution to still be nearby the prior one. One technique is to simply use the optimal point from one solution as the starting point for the next step. Now more complex techniques are possible that do some extrapolation to compute an initial position based on prior solution(s) rather than directly use the prior solution.

## Expectation¶

The energy of the Hamiltonian operator that VQE is working on is the expectation value when evaluated with the parameterized ansatz. To compute the expectation value VQE uses an instance of an expectation object. Such an instance may be supplied via the expectation parameter, or in the default case, where it has a value of None, VQE will use the ExpectationFactory to create itself a suitable instance based on the supplied backend.

For most cases letting VQE create a suitable instance is sufficient. However the Qiskit Aer aer_simulator supports a snapshot instruction that can be used in conjunction with the operator expectation computation. If used then the outcome is ideal, i.e. like the statevector simulator, and has no shot noise. Since people normally choose the aer_simulator to have shot noise (sampling noise), and be more like a real-device outcome, VQE has an include_custom flag that is passed on to the ExpectationFactory. When using Aer qasm simulator, and this is set True, the factory will return AerPauliExpectation which uses the snapshot instruction, when False, default, then the regular PauliExpectation is returned.

The following example shows include_custom=True where the outcome matches the statevector simulator. In fact it can be better/faster to do this than use the statevector_simulator directly. This is because in the latter case when the Hamiltonian is a sum of Paulis it must be converted to matrix form, and this is avoided when using the snapshot instruction done when include_custom is True.

[5]:

from qiskit import Aer

algorithm_globals.random_seed = seed
qi = QuantumInstance(Aer.get_backend('aer_simulator'), seed_transpiler=seed, seed_simulator=seed)

ansatz = TwoLocal(rotation_blocks='ry', entanglement_blocks='cz')
slsqp = SLSQP(maxiter=1000)
vqe = VQE(ansatz, optimizer=slsqp, quantum_instance=qi, include_custom=True)
result = vqe.compute_minimum_eigenvalue(operator=H2_op)
optimal_value1 = result.optimal_value
print(result)

{   'aux_operator_eigenvalues': None,
'cost_function_evals': 65,
'eigenstate': {'01': 0.9921567416492215, '10': 0.125},
'eigenvalue': (-1.8572750175571595+0j),
'optimal_parameters': {   ParameterVectorElement(θ[0]): 4.2965205340503685,
ParameterVectorElement(θ[1]): 4.426962242132452,
ParameterVectorElement(θ[2]): 0.5470754193210003,
ParameterVectorElement(θ[7]): 0.36020735708081203,
ParameterVectorElement(θ[6]): -4.717618177195121,
ParameterVectorElement(θ[5]): 1.5683259454122547,
ParameterVectorElement(θ[4]): -2.5983258639687397,
ParameterVectorElement(θ[3]): 6.092947857528147},
'optimal_point': array([ 4.29652053,  4.42696224,  0.54707542,  6.09294786, -2.59832586,
1.56832595, -4.71761818,  0.36020736]),
'optimal_value': -1.8572750175571595,
'optimizer_evals': None,
'optimizer_time': 0.10388422012329102}


In case you have doubts here is the aer_simulator again but include_custom has been left to default to False. The optimization ended abruptly, presumably due to the shot noise confusing the SLSQP optimizer.

[6]:

algorithm_globals.random_seed = seed
qi = QuantumInstance(Aer.get_backend('aer_simulator'), seed_transpiler=seed, seed_simulator=seed)

ansatz = TwoLocal(rotation_blocks='ry', entanglement_blocks='cz')
slsqp = SLSQP(maxiter=1000)
vqe = VQE(ansatz, optimizer=slsqp, quantum_instance=qi)
result = vqe.compute_minimum_eigenvalue(operator=H2_op)
optimal_value = result.optimal_value
print('The optimal value can be seen to be wrong too, i.e. '
f'{optimal_value:.3f} versus the correct {optimal_value1:.3f}.')
print()
print(result)

The optimal value can be seen to be wrong too, i.e. -1.099 versus the correct -1.857.

{   'aux_operator_eigenvalues': None,
'cost_function_evals': 9,
'eigenstate': {   '00': 0.7781187248742958,
'01': 0.4881406047441659,
'10': 0.39404750665370286,
'11': 0.03125},
'eigenvalue': (-1.0987888676631705+0j),
'optimal_parameters': {   ParameterVectorElement(θ[2]): 0.6019852007557844,
ParameterVectorElement(θ[4]): -3.3070470445355764,
ParameterVectorElement(θ[5]): 1.8462931831829383,
ParameterVectorElement(θ[6]): -5.466043598406607,
ParameterVectorElement(θ[7]): 0.6984088030463615,
ParameterVectorElement(θ[0]): 3.611860069224077,
ParameterVectorElement(θ[1]): 4.19301252102391,
ParameterVectorElement(θ[3]): 5.949536809130025},
'optimal_point': array([ 3.61186007,  4.19301252,  0.6019852 ,  5.94953681, -3.30704704,
1.84629318, -5.4660436 ,  0.6984088 ]),
'optimal_value': -1.0987888676631705,
'optimizer_evals': None,
'optimizer_time': 0.34603023529052734}


Changing the optimizer to SPSA, which is designed to work in noisy environments, gets us a better result. Though the noise has affected the outcome so it’s not as accurate.

[7]:

from qiskit.algorithms.optimizers import SPSA

algorithm_globals.random_seed = seed
qi = QuantumInstance(Aer.get_backend('aer_simulator'), seed_transpiler=seed, seed_simulator=seed)

ansatz = TwoLocal(rotation_blocks='ry', entanglement_blocks='cz')
slsqp = SPSA(maxiter=100)
vqe = VQE(ansatz, optimizer=slsqp, quantum_instance=qi)
result = vqe.compute_minimum_eigenvalue(operator=H2_op)
print(result)

{   'aux_operator_eigenvalues': None,
'cost_function_evals': 200,
'eigenstate': {'01': 0.9916644782889019, '10': 0.1288470508005519},
'eigenvalue': (-1.8627691667457411+0j),
'optimal_parameters': {   ParameterVectorElement(θ[6]): -4.750881193409766,
ParameterVectorElement(θ[4]): -2.365851315862917,
ParameterVectorElement(θ[3]): 6.659266627719733,
ParameterVectorElement(θ[5]): 2.9099466356599697,
ParameterVectorElement(θ[7]): 0.17033837870153928,
ParameterVectorElement(θ[2]): -0.274076481794564,
ParameterVectorElement(θ[1]): 3.413112968624178,
ParameterVectorElement(θ[0]): 3.702453458416716},
'optimal_point': array([ 3.70245346,  3.41311297, -0.27407648,  6.65926663, -2.36585132,
2.90994664, -4.75088119,  0.17033838]),
'optimal_value': -1.8627691667457411,
'optimizer_evals': None,
'optimizer_time': 7.184031963348389}


As mentioned above, an expectation object can be explicitly given (so the internal ExpectationFactory and include_custom are never used/needed. Below we create an AerPauliExpectation and pass this to VQE. We can see the result matches that above where we set include_custom to True and let VQE create its own expectation object.

[8]:

from qiskit.opflow import AerPauliExpectation

algorithm_globals.random_seed = seed
qi = QuantumInstance(Aer.get_backend('aer_simulator'), seed_transpiler=seed, seed_simulator=seed)

ansatz = TwoLocal(rotation_blocks='ry', entanglement_blocks='cz')
slsqp = SLSQP(maxiter=1000)
vqe = VQE(ansatz, optimizer=slsqp, quantum_instance=qi,
expectation=AerPauliExpectation())
result = vqe.compute_minimum_eigenvalue(operator=H2_op)
print(result)

{   'aux_operator_eigenvalues': None,
'cost_function_evals': 65,
'eigenstate': {'01': 0.9921567416492215, '10': 0.125},
'eigenvalue': (-1.8572750175571595+0j),
'optimal_parameters': {   ParameterVectorElement(θ[5]): 1.5683259454122547,
ParameterVectorElement(θ[3]): 6.092947857528147,
ParameterVectorElement(θ[7]): 0.36020735708081203,
ParameterVectorElement(θ[6]): -4.717618177195121,
ParameterVectorElement(θ[4]): -2.5983258639687397,
ParameterVectorElement(θ[2]): 0.5470754193210003,
ParameterVectorElement(θ[0]): 4.2965205340503685,
ParameterVectorElement(θ[1]): 4.426962242132452},
'optimal_point': array([ 4.29652053,  4.42696224,  0.54707542,  6.09294786, -2.59832586,
1.56832595, -4.71761818,  0.36020736]),
'optimal_value': -1.8572750175571595,
'optimizer_evals': None,
'optimizer_time': 0.10378384590148926}


By default, the PauliExpectation object, that would have be chosen when include_custom is False (or when using Aer aer_simulator, or a real device) groups Paulis into commuting sets. This is efficient as it runs less circuits to compute the expectation. However, if for some reason you wanted to run a circuit for each Pauli then then grouping can be turned off when constructing the PauliExpectation. You need to explicitly pass in such an expectation instance to VQE to have it work this way though as shown below.

[9]:

from qiskit.opflow import PauliExpectation

algorithm_globals.random_seed = seed
qi = QuantumInstance(Aer.get_backend('aer_simulator'), seed_transpiler=seed, seed_simulator=seed)

ansatz = TwoLocal(rotation_blocks='ry', entanglement_blocks='cz')
slsqp = SPSA(maxiter=100)
vqe = VQE(ansatz, optimizer=slsqp, quantum_instance=qi,
expectation=PauliExpectation(group_paulis=False))
result = vqe.compute_minimum_eigenvalue(operator=H2_op)
print(result)

{   'aux_operator_eigenvalues': None,
'cost_function_evals': 200,
'eigenstate': {'01': 0.9916644782889019, '10': 0.1288470508005519},
'eigenvalue': (-1.8649830308114583+0j),
'optimal_parameters': {   ParameterVectorElement(θ[7]): 0.18495131387681987,
ParameterVectorElement(θ[1]): 3.3827104248942383,
ParameterVectorElement(θ[4]): -2.282559635034206,
ParameterVectorElement(θ[2]): -0.26516964794294456,
ParameterVectorElement(θ[3]): 6.699789960055847,
ParameterVectorElement(θ[0]): 3.683395848045204,
ParameterVectorElement(θ[6]): -4.766514920799031,
ParameterVectorElement(θ[5]): 2.9887235542230135},
'optimal_point': array([ 3.68339585,  3.38271042, -0.26516965,  6.69978996, -2.28255964,
2.98872355, -4.76651492,  0.18495131]),
'optimal_value': -1.8649830308114583,
'optimizer_evals': None,
'optimizer_time': 10.103214740753174}


Optimizers that use a gradient-based technique can be supplied with a user defined gradient that will be used instead of their default gradient computation which is usually done by simple finite difference. Gradients are passed indirectly via to the optimizer via its gradient parameter.

As the use of a user supplied gradient was shown in the Monitoring VQE Convergence tutorial I will simply refer you there. Also the Gradients framework tutorial has much more about the gradients themselves.

## Quantum Instance and advanced simulation¶

While you may be familiar with passing a QuantumInstancen created from a aer_simulator_statevector a aer_simulator or real device backend, it is possible to use the advanced simulation modes of Aer too when applicable. For instance we can easily use the Aer Matrix Product State method, that has the potential to scale to larger numbers of qubits.

[10]:

algorithm_globals.random_seed = seed

from qiskit.providers.aer import QasmSimulator
quantum_instance = QuantumInstance(QasmSimulator(method='matrix_product_state'), shots=1)

ansatz = TwoLocal(rotation_blocks='ry', entanglement_blocks='cz')
slsqp = SLSQP(maxiter=1000)
vqe = VQE(ansatz, optimizer=slsqp, quantum_instance=qi, include_custom=True)
result = vqe.compute_minimum_eigenvalue(operator=H2_op)
print(result)

{   'aux_operator_eigenvalues': None,
'cost_function_evals': 65,
'eigenstate': {'01': 0.9921567416492215, '10': 0.125},
'eigenvalue': (-1.8572750175571595+0j),
'optimal_parameters': {   ParameterVectorElement(θ[2]): 0.5470754193210003,
ParameterVectorElement(θ[4]): -2.5983258639687397,
ParameterVectorElement(θ[7]): 0.36020735708081203,
ParameterVectorElement(θ[5]): 1.5683259454122547,
ParameterVectorElement(θ[6]): -4.717618177195121,
ParameterVectorElement(θ[3]): 6.092947857528147,
ParameterVectorElement(θ[1]): 4.426962242132452,
ParameterVectorElement(θ[0]): 4.2965205340503685},
'optimal_point': array([ 4.29652053,  4.42696224,  0.54707542,  6.09294786, -2.59832586,
1.56832595, -4.71761818,  0.36020736]),
'optimal_value': -1.8572750175571595,
'optimizer_evals': None,
'optimizer_time': 0.10611629486083984}

[11]:

import qiskit.tools.jupyter
%qiskit_version_table


### Version Information

Qiskit SoftwareVersion
qiskit-terra0.21.0
qiskit-aer0.10.4
qiskit-ibmq-provider0.19.2
qiskit0.37.0
qiskit-nature0.4.1
qiskit-finance0.3.3
qiskit-optimization0.4.0
qiskit-machine-learning0.4.0
System information
Python version3.8.13
Python compilerGCC 9.4.0
Python builddefault, Jun 20 2022 14:28:56
OSLinux
CPUs2
Memory (Gb)6.783603668212891
Thu Jun 30 18:27:32 2022 UTC

### This code is a part of Qiskit

[ ]: