Instancias y extensiones

En este capítulo se presentan varios algoritmos variacionales cuánticos, entre ellos:

• Variational Quantum Eigensolver (VQE)
• Subspace Search VQE (SSVQE)
• Variational Quantum Deflation (VQD)
• Quantum Sampling Regression (QSR)

A través de estos algoritmos aprenderás varias ideas de diseño que pueden incorporarse en algoritmos variacionales personalizados, como pesos, penalizaciones, sobremuestreo y submuestreo. Te animamos a experimentar con estos conceptos y a compartir tus hallazgos con la comunidad.

El marco de patrones de Qiskit aplica a todos estos algoritmos, pero solo señalaremos los pasos explícitamente en el primer ejemplo.

Variational Quantum Eigensolver (VQE)

VQE es uno de los algoritmos variacionales cuánticos más utilizados, y establece una plantilla sobre la que se construyen otros algoritmos.

Paso 1: Mapear las entradas clásicas a un problema cuántico

Esquema teórico

El esquema de VQE es sencillo:

• Preparar los operadores de referencia $U_R$
  • Partimos del estado $|0\rangle$ y llegamos al estado de referencia $|\rho\rangle$
• Aplicar la forma variacional $U_V(\vec\theta_{i,j})$ para crear un ansatz $U_A(\vec\theta_{i,j})$
  • Pasamos del estado $|\rho\rangle$ a $U_V(\vec\theta_{i,j})|\rho\rangle = |\psi(\vec\theta_{i,j})\rangle$
• Arranque en $i=0$ si tenemos un problema similar (normalmente obtenido mediante simulación clásica o muestreo)
  • Cada optimizador se inicializará de forma diferente, generando un conjunto inicial de vectores de parámetros $\Theta_0 := \\{ {\vec\theta_{0,j} | j \in \mathcal{J}_\text{opt}^0} \\}$ (por ejemplo, a partir de un punto inicial $\vec\theta_0$).
• Evaluar la función de costo $C(\vec\theta_{i,j}) := \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ para todos los estados preparados en un ordenador cuántico.
• Usar un optimizador clásico para seleccionar el siguiente conjunto de parámetros $\Theta_{i+1}$.
• Repetir el proceso hasta alcanzar la convergencia.

Este es un ciclo de optimización clásico sencillo en el que evaluamos la función de costo. Algunos optimizadores pueden requerir varias evaluaciones para calcular un gradiente, determinar la siguiente iteración o verificar la convergencia.

A continuación se muestra el ejemplo para el siguiente observable:

$$\hat{O}_1 = 2 II - 2 XX + 3 YY - 3 ZZ,$$

Implementación

# Added by doQumentation — required packages for this notebook
!pip install -q numpy qiskit qiskit-ibm-runtime scipy

from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit.circuit.library import TwoLocal
import numpy as np

theta_list = (2 * np.pi * np.random.rand(1, 8)).tolist()
observable = SparsePauliOp.from_list([("II", 2), ("XX", -2), ("YY", 3), ("ZZ", -3)])

reference_circuit = QuantumCircuit(2)
reference_circuit.x(0)

variational_form = TwoLocal(
    2,
    rotation_blocks=["rz", "ry"],
    entanglement_blocks="cx",
    entanglement="linear",
    reps=1,
)

ansatz = reference_circuit.compose(variational_form)

ansatz.decompose().draw("mpl")

def cost_func_vqe(parameters, ansatz, hamiltonian, estimator):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance

    Returns:
        float: Energy estimate
    """

    estimator_job = estimator.run([(ansatz, hamiltonian, [parameters])])
    estimator_result = estimator_job.result()[0]

    cost = estimator_result.data.evs[0]
    return cost

from qiskit.primitives import StatevectorEstimator

estimator = StatevectorEstimator()

Podemos usar esta función de costo para calcular los parámetros óptimos:

# SciPy minimizer routine
from scipy.optimize import minimize

x0 = np.ones(8)

result = minimize(
    cost_func_vqe, x0, args=(ansatz, observable, estimator), method="COBYLA"
)

result

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -5.999999982445723
       x: [ 1.741e+00  9.606e-01  1.571e+00  2.115e-05  1.899e+00
            1.243e+00  6.063e-01  6.063e-01]
    nfev: 136
   maxcv: 0.0

Paso 2: Optimizar el problema para la ejecución cuántica

Seleccionaremos el backend menos ocupado e importaremos los componentes necesarios de qiskit_ibm_runtime.

from qiskit_ibm_runtime import SamplerV2 as Sampler
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit_ibm_runtime import Session, EstimatorOptions
from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)
print(backend)

<IBMBackend('ibm_brisbane')>

Transpilaremos el circuito usando el gestor de pasadas predefinido con nivel de optimización 3, y aplicaremos el diseño correspondiente al observable.

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
isa_ansatz = pm.run(ansatz)
isa_observable = observable.apply_layout(layout=isa_ansatz.layout)

Paso 3: Ejecutar usando las primitivas de Qiskit Runtime

Ya estamos listos para ejecutar nuestro cálculo en hardware de IBM Quantum®. Dado que la minimización de la función de costo es muy iterativa, iniciaremos una sesión de Runtime. De este modo, solo tendremos que esperar en la cola una vez. Una vez que el trabajo comience a ejecutarse, cada iteración con actualizaciones de parámetros se ejecutará de inmediato.

x0 = np.ones(8)

estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    result = minimize(
        cost_func_vqe,
        x0,
        args=(isa_ansatz, isa_observable, estimator),
        method="COBYLA",
        options={"maxiter": 200, "disp": True},
    )
session.close()
print(result)

Paso 4: Posprocesar y devolver el resultado en formato clásico

Podemos ver que la rutina de minimización terminó con éxito, lo que significa que se alcanzó la tolerancia predeterminada del optimizador clásico COBYLA. Si necesitamos un resultado más preciso, podemos especificar una tolerancia menor. Esto podría ser necesario, ya que el resultado se alejó varios puntos porcentuales del obtenido por el simulador anterior.

El valor de x obtenido es la mejor estimación actual de los parámetros que minimizan la función de costo. Si se itera para obtener mayor precisión, esos valores deben usarse en lugar del x0 inicial (un vector de unos).

Por último, cabe señalar que la función fue evaluada 96 veces durante el proceso de optimización. Este número puede diferir de la cantidad de pasos de optimización, ya que algunos optimizadores requieren múltiples evaluaciones de la función en un solo paso, como ocurre al estimar un gradiente.

Subspace Search VQE (SSVQE)

SSVQE es una variante de VQE que permite obtener los primeros $k$ valores propios de un observable $\hat{H}$ con valores propios $\{\lambda_0, \lambda_1,...,\lambda_{N-1}\}$, donde $N\geq k$. Sin pérdida de generalidad, suponemos que $\lambda_0<\lambda_1<...<\lambda_{N-1}$. SSVQE introduce una nueva idea al añadir pesos para priorizar la optimización del término con mayor peso.

Para implementar este algoritmo, necesitamos $k$ estados de referencia mutuamente ortogonales $\{ |\rho_j\rangle \}_{j=0}^{k-1}$, lo que significa que $\langle \rho_j | \rho_l \rangle = \delta_{jl}$ para $j,l<k$. Estos estados pueden construirse usando operadores de Pauli. La función de costo de este algoritmo es:

$$\begin{aligned} C(\vec{\theta}) & := \sum_{j=0}^{k-1} w_j \langle \rho_j | U_{V}^{\dagger}(\vec{\theta})\hat{H} U_{V}(\vec{\theta})|\rho_j \rangle \\[1mm] & := \sum_{j=0}^{k-1} w_j \langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle \\[1mm] \end{aligned}$$

donde $w_j$ es un número positivo arbitrario tal que si $j<l<k$ entonces $w_j>w_l$, y $U_V(\vec{\theta})$ es la forma variacional definida por el usuario.

El algoritmo SSVQE se basa en el hecho de que los estados propios correspondientes a distintos valores propios son mutuamente ortogonales. En particular, el producto interno de $U_V(\vec{\theta})|\rho_j\rangle$ y $U_V(\vec{\theta})|\rho_l\rangle$ puede expresarse como:

$$\begin{aligned} \langle \rho_j | U_{V}^{\dagger}(\vec{\theta})U_{V}(\vec{\theta})|\rho_l \rangle & = \langle \rho_j | I |\rho_l \rangle \\[1mm] & = \langle \rho_j | \rho_l \rangle \\[1mm] & = \delta_{jl} \end{aligned}$$

La primera igualdad es válida porque $U_{V}(\vec{\theta})$ es un operador cuántico y, por tanto, unitario. La última igualdad se cumple por la ortogonalidad de los estados de referencia $|\rho_j\rangle$. El hecho de que la ortogonalidad se conserve a través de transformaciones unitarias está profundamente relacionado con el principio de conservación de la información, tal como se expresa en la ciencia de la información cuántica. Desde esta perspectiva, las transformaciones no unitarias representan procesos en los que se pierde o se inyecta información.

Los pesos $w_j$ ayudan a garantizar que todos los estados sean estados propios. Si los pesos son suficientemente distintos, el término con mayor peso (es decir, $w_0$) recibirá prioridad durante la optimización sobre los demás. Como resultado, el estado resultante $U_{V}(\vec{\theta})|\rho_0 \rangle$ se convertirá en el estado propio correspondiente a $\lambda_0$. Dado que $\{ U_{V}(\vec{\theta})|\rho_j\rangle \}_{j=0}^{k-1}$ son mutuamente ortogonales, los estados restantes serán ortogonales a este y, por tanto, estarán contenidos en el subespacio correspondiente a los valores propios $\{\lambda_1,...,\lambda_{N-1}\}$.

Aplicando el mismo argumento al resto de los términos, la siguiente prioridad sería el término con peso $w_1$, de modo que $U_{V}(\vec{\theta})|\rho_1 \rangle$ sería el estado propio correspondiente a $\lambda_1$, y los demás términos estarían contenidos en el espacio propio de $\{\lambda_2,...,\lambda_{N-1}\}$.

Razonando inductivamente, deducimos que $U_{V}(\vec{\theta})|\rho_j \rangle$ será un estado propio aproximado de $\lambda_j$ para $0\leq j < k.$

Esquema teórico

El algoritmo SSVQE puede resumirse de la siguiente manera:

• Preparar varios estados de referencia aplicando un unitario U_R a $k$ estados distintos de la base computacional
  • Este algoritmo requiere el uso de $k$ estados de referencia mutuamente ortogonales $\{ |\rho_j\rangle \}_{j=0}^{k-1}$, de modo que $\langle \rho_j | \rho_l \rangle = \delta_{jl}$ para $j,l<k$.
• Aplicar la forma variacional $U_V(\vec\theta_{i,j})$ a cada estado de referencia, obteniendo el siguiente ansatz $U_A(\vec\theta_{i,j})$.
• Arranque en $i=0$ si se dispone de un problema similar (normalmente obtenido mediante simulación clásica o muestreo).
• Evaluar la función de costo $C(\vec\theta_{i,j}) := \sum_{j=0}^{k-1} w_j \langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ para todos los estados preparados en un ordenador cuántico.
  • Esto puede separarse en el cálculo del valor esperado de un observable $\langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ y la multiplicación del resultado por $w_j$.
  • A continuación, la función de costo devuelve la suma de todos los valores esperados ponderados.
• Usar un optimizador clásico para determinar el siguiente conjunto de parámetros $\Theta_{i+1}$.
• Repetir los pasos anteriores hasta alcanzar la convergencia.

Reconstruirás la función de costo de SSVQE en la evaluación, pero aquí tienes el siguiente fragmento para orientar tu solución:

import numpy as np

def cost_func_ssvqe(
    params, initialized_anastz_list, weights, ansatz, hamiltonian, estimator
):
    # """Return estimate of energy from estimator

    # Parameters:
    #     params (ndarray): Array of ansatz parameters
    #     initialized_anastz_list (list QuantumCircuit): Array of initialised ansatz with reference
    #     weights (list): List of weights
    #     ansatz (QuantumCircuit): Parameterized ansatz circuit
    #     hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
    #     estimator (Estimator): Estimator primitive instance

    # Returns:
    #     float: Weighted energy estimate
    # """

    energies = []

    # Define SSVQE

    weighted_energy_sum = np.dot(energies, weights)
    return weighted_energy_sum

Deflación cuántica variacional (VQD)

VQD es un método iterativo que extiende VQE para obtener los primeros $k$ valores propios de un observable $\hat{H}$ con valores propios $\{\lambda_0, \lambda_1,...,\lambda_{N-1}\}$, donde $N\geq k$, en lugar de solo el primero. Para el resto de esta sección, asumiremos, sin pérdida de generalidad, que $\lambda_0\leq\lambda_1\leq...\leq\lambda_{N-1}$. VQD introduce la noción de un coste de penalización para guiar el proceso de optimización.

VQD introduce un término de penalización, denotado como $\beta$, para equilibrar la contribución de cada término de solapamiento al coste. Este término de penalización sirve para sancionar el proceso de optimización si no se logra la ortogonalidad. Imponemos esta restricción porque los estados propios de un observable, o de un operador hermítico, correspondientes a distintos valores propios son siempre mutuamente ortogonales, o pueden hacerse así en el caso de degeneración o valores propios repetidos. Así, al imponer la ortogonalidad con el estado propio correspondiente a $\lambda_0$, estamos optimizando efectivamente sobre el subespacio que corresponde al resto de los valores propios $\{\lambda_1, \lambda_2,..., \lambda_{N-1}\}$. Aquí, $\lambda_1$ es el valor propio más bajo de entre los restantes y, por tanto, la solución óptima del nuevo problema puede obtenerse usando el teorema variacional.

La idea general detrás de VQD es usar VQE de la forma habitual para obtener el valor propio más bajo $\lambda_0 := C_0(\vec\theta^0) \equiv C_\text{VQE}(\vec\theta^0)$ junto con el estado propio (aproximado) correspondiente $|\psi(\vec{\theta^0})\rangle$ para algún vector de parámetros óptimo $\vec{\theta^0}$. Luego, para obtener el siguiente valor propio $\lambda_1 > \lambda_0$, en lugar de minimizar la función de coste $C_0(\vec{\theta}) := \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$, optimizamos:

$$C_1(\vec{\theta}) := C_0(\vec{\theta})+ \beta_0 |\langle \psi(\vec{\theta})| \psi(\vec{\theta^0})\rangle |^2$$

El valor positivo $\beta_0$ debería ser idealmente mayor que $\lambda_1-\lambda_0$.

Esto introduce una nueva función de coste que puede verse como un problema con restricciones, donde minimizamos $C_\text{VQE}(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ sujeto a la restricción de que el estado debe ser ortogonal al $|\psi(\vec{\theta^0})\rangle$ previamente obtenido, con $\beta_0$ actuando como término de penalización si la restricción no se satisface.

Alternativamente, este nuevo problema puede interpretarse como ejecutar VQE sobre el observable:

$$\hat{H_1} := \hat{H} + \beta_0 |\psi(\vec{\theta^0})\rangle \langle \psi(\vec{\theta^0})| \quad \Rightarrow \quad C_1(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H_1} | \psi(\vec{\theta})\rangle,$$

Asumiendo que la solución al nuevo problema es $|\psi(\vec{\theta^1})\rangle$, el valor esperado de $\hat{H}$ (no de $\hat{H_1}$) debería ser $\langle \psi(\vec{\theta^1}) | \hat{H} | \psi(\vec{\theta^1})\rangle = \lambda_1$. Para obtener el tercer valor propio $\lambda_2$, la función de coste a optimizar es:

$$C_2(\vec{\theta}) := C_1(\vec{\theta}) + \beta_1 |\langle \psi(\vec{\theta})| \psi(\vec{\theta^1})\rangle |^2$$

donde $\beta_1$ es una constante positiva suficientemente grande para imponer la ortogonalidad del estado solución tanto con $|\psi(\vec{\theta^0})\rangle$ como con $|\psi(\vec{\theta^1})\rangle$. Esto penaliza los estados del espacio de búsqueda que no cumplen este requisito, restringiendo efectivamente dicho espacio. Así, la solución óptima del nuevo problema debería ser el estado propio correspondiente a $\lambda_2$.

Al igual que en el caso anterior, este nuevo problema también puede interpretarse como VQE con el observable:

$$\hat{H_2} := \hat{H_1} + \beta_1 |\psi(\vec{\theta^1})\rangle \langle \psi(\vec{\theta^1})| \quad \Rightarrow \quad C_2(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H_2} | \psi(\vec{\theta})\rangle.$$

Si la solución a este nuevo problema es $|\psi(\vec{\theta^2})\rangle$, el valor esperado de $\hat{H}$ (no de $\hat{H_2}$) debería ser $\langle \psi(\vec{\theta^2}) | \hat{H} | \psi(\vec{\theta^2})\rangle = \lambda_2$. De forma análoga, para obtener el $k$-ésimo valor propio $\lambda_{k-1}$, minimizarías la función de coste:

$$C_{k-1}(\vec{\theta}) := C_{k-2}(\vec{\theta}) + \beta_{k-2} |\langle \psi(\vec{\theta})| \psi(\vec{\theta^{k-2}})\rangle |^2,$$

Recuerda que definimos $\vec{\theta^j}$ de forma que $\langle \psi(\vec{\theta^j}) | \hat{H} | \psi(\vec{\theta^j})\rangle = \lambda_j, \forall j<k$. Este problema es equivalente a minimizar $C(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ pero con la restricción de que el estado debe ser ortogonal a $|\psi(\vec{\theta^j})\rangle ; \forall j \in {0, \cdots, k-1}$, restringiendo así el espacio de búsqueda al subespacio correspondiente a los valores propios $\{\lambda_{k-1},\cdots,\lambda_{N-1}\}$.

Este problema es equivalente a un VQE con el observable:

$$\hat{H}_{k-1} := \hat{H}_{k-2} + \beta_{k-2} |\psi(\vec{\theta^{k-2}})\rangle \langle \psi(\vec{\theta^{k-2}})| \quad \Rightarrow \quad C_{k-1}(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H}_{k-1} | \psi(\vec{\theta})\rangle,$$

Como puedes ver en el proceso, para obtener el $k$-ésimo valor propio necesitas los estados propios (aproximados) de los $k-1$ valores propios anteriores, por lo que tendrías que ejecutar VQE un total de $k$ veces. Por tanto, la función de coste de VQD es la siguiente:

$$C_k(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle + \sum_{j=0}^{k-1}\beta_j |\langle \psi(\vec{\theta})| \psi(\vec{\theta^j})\rangle |^2$$

Esquema teórico

El esquema de VQD puede resumirse así:

• Preparar un operador de referencia $U_R$
• Aplicar la forma variacional $U_V(\vec\theta_{i,j})$ al estado de referencia, creando los siguientes ansätze $U_A(\vec\theta_{i,j})$
• Arrancar en $i=0$ si se dispone de un problema similar (típicamente hallado mediante simulación clásica o muestreo).
• Evaluar la función de coste $C_k(\vec{\theta})$, que implica calcular $k$ estados excitados y un array de $\beta$'s que define la penalización por solapamiento para cada término.
  • Calcular el valor esperado de un observable $\langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ para cada $k$
  • Calcular la penalización $\sum_{j=0}^{k-1}\beta_j |\langle \psi(\vec{\theta})| \psi(\vec{\theta^j})\rangle |^2$.
  • La función de coste debe devolver la suma de estos dos términos
• Usar un optimizador clásico para elegir el siguiente conjunto de parámetros $\Theta_{i+1}$.
• Repetir este proceso hasta alcanzar la convergencia.

Implementación

Para esta implementación, crearemos una función para una penalización por solapamiento. Esta penalización se usará en la función de coste en cada iteración. Este proceso se repetirá para cada estado excitado.

from qiskit.circuit.library import TwoLocal

ansatz = TwoLocal(2, rotation_blocks=["ry", "rz"], entanglement_blocks="cz", reps=1)

ansatz.decompose().draw("mpl")

Primero, definiremos una función que calcula la fidelidad del estado — un porcentaje de solapamiento entre dos estados que usaremos como penalización para VQD:

import numpy as np

def calculate_overlaps(ansatz, prev_circuits, parameters, sampler):
    def create_fidelity_circuit(circuit_1, circuit_2):
        """
        Constructs the list of fidelity circuits to be evaluated.
        These circuits represent the state overlap between pairs of input circuits,
        and their construction depends on the fidelity method implementations.
        """

        if len(circuit_1.clbits) > 0:
            circuit_1.remove_final_measurements()
        if len(circuit_2.clbits) > 0:
            circuit_2.remove_final_measurements()

        circuit = circuit_1.compose(circuit_2.inverse())
        circuit.measure_all()
        return circuit

    overlaps = []

    for prev_circuit in prev_circuits:
        fidelity_circuit = create_fidelity_circuit(ansatz, prev_circuit)
        sampler_job = sampler.run([(fidelity_circuit, parameters)])
        meas_data = sampler_job.result()[0].data.meas

        counts_0 = meas_data.get_int_counts().get(0, 0)
        shots = meas_data.num_shots
        overlap = counts_0 / shots
        overlaps.append(overlap)

    return np.array(overlaps)

Es el momento de escribir la función de coste de VQD. Como antes, cuando calculamos solo el estado fundamental, determinaremos el estado de menor energía usando la primitiva Estimator. Sin embargo, tal como se describió más arriba, ahora añadiremos un término de penalización para garantizar la ortogonalidad de los estados de mayor energía. Es decir, para cada nuevo estado excitado, se añade una penalización por cualquier solapamiento entre el estado variacional actual y los estados propios de menor energía ya encontrados.

def cost_func_vqd(
    parameters, ansatz, prev_states, step, betas, estimator, sampler, hamiltonian
):
    estimator_job = estimator.run([(ansatz, hamiltonian, [parameters])])

    total_cost = 0

    if step > 1:
        overlaps = calculate_overlaps(ansatz, prev_states, parameters, sampler)
        total_cost = np.sum(
            [np.real(betas[state] * overlap) for state, overlap in enumerate(overlaps)]
        )

    estimator_result = estimator_job.result()[0]

    value = estimator_result.data.evs[0] + total_cost

    return value

Nota especialmente que la función de coste anterior hace referencia a la función calculate_overlaps, que en realidad crea un nuevo circuito cuántico. Si queremos ejecutar esto en hardware real, ese nuevo circuito también debe ser transpilado, idealmente de forma óptima, para ejecutarse en el backend que seleccionemos. Observa que la transpilación no está incorporada en las funciones calculate_overlaps ni cost_func_vqd. Puedes intentar modificar el código tú mismo para añadir esta transpilación adicional (y condicional) — pero esto también se hará por ti en la siguiente lección.

En esta lección, ejecutaremos el algoritmo VQD usando el Statevector Sampler y el Statevector Estimator:

from qiskit.primitives import StatevectorEstimator as Estimator

sampler = Sampler()
estimator = Estimator()

Introduciremos un observable a estimar. En la siguiente lección le daremos un contexto físico, como el estado excitado de una molécula. Puede ser útil pensar en este observable como el hamiltoniano de un sistema que puede tener estados excitados, aunque este observable no se ha elegido para representar ninguna molécula o átomo en particular.

from qiskit.quantum_info import SparsePauliOp

observable = SparsePauliOp.from_list([("II", 2), ("XX", -2), ("YY", 3), ("ZZ", -3)])

Aquí definimos el número total de estados que deseamos calcular (estado fundamental y estados excitados, k), y las penalizaciones (betas) por solapamiento entre vectores de estado que deberían ser ortogonales. Las consecuencias de elegir betas demasiado altas o demasiado bajas se explorarán un poco en la siguiente lección. Por ahora, usaremos simplemente las que se proporcionan a continuación. Comenzaremos usando ceros como parámetros iniciales. En tus propios cálculos, quizás quieras usar parámetros iniciales más inteligentes basados en tu conocimiento del sistema o en cálculos anteriores.

k = 3
betas = [33, 33, 33]
x0 = np.zeros(8)

Ya podemos ejecutar el cálculo:

from scipy.optimize import minimize

prev_states = []
prev_opt_parameters = []
eigenvalues = []

for step in range(1, k + 1):
    if step > 1:
        prev_states.append(ansatz.assign_parameters(prev_opt_parameters))

    result = minimize(
        cost_func_vqd,
        x0,
        args=(ansatz, prev_states, step, betas, estimator, sampler, observable),
        method="COBYLA",
        options={
            "maxiter": 200,
        },
    )
    print(result)

    prev_opt_parameters = result.x
    eigenvalues.append(result.fun)

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -5.999999979545955
       x: [-5.150e-01 -5.452e-02 -1.571e+00 -2.853e-05  2.671e-01
           -2.672e-01 -8.509e-01 -8.510e-01]
    nfev: 131
   maxcv: 0.0
 message: Optimization terminated successfully.
 success: True
  status: 1
     fun: 4.024550284767612
       x: [-3.745e-01  1.041e+00  8.637e-01  1.202e+00 -8.847e-02
            1.181e-02  7.611e-01 -3.006e-01]
    nfev: 110
   maxcv: 0.0
 message: Optimization terminated successfully.
 success: True
  status: 1
     fun: 5.608925562838559
       x: [-2.670e-01  1.280e+00  1.070e+00 -8.031e-01 -1.524e-01
           -6.956e-02  7.018e-01  1.514e+00]
    nfev: 90
   maxcv: 0.0

Los valores que obtuvimos de la función de coste son aproximadamente -6,00, 4,02 y 5,61. Lo importante de estos resultados es que los valores de la función van aumentando. Si hubiéramos obtenido un primer estado excitado con una energía menor que el cálculo inicial sin restricciones del estado fundamental, eso habría indicado un error en algún lugar de nuestro código.

Los valores de x son los parámetros que produjeron un vector de estado correspondiente a cada uno de estos costes (energías).

Por último, observamos que las tres minimizaciones convergieron dentro de la tolerancia predeterminada del optimizador clásico (aquí COBYLA). Requirieron 131, 110 y 90 evaluaciones de función, respectivamente.

Regresión por muestreo cuántico (QSR)

Uno de los principales problemas de VQE es la gran cantidad de llamadas a un ordenador cuántico necesarias para obtener los parámetros en cada paso, por ejemplo, $k$, $k-1$, etc. Esto es especialmente problemático cuando el acceso a los dispositivos cuánticos está en cola. Aunque se puede usar una Session para agrupar múltiples llamadas iterativas, un enfoque alternativo es el muestreo. Al aprovechar más recursos clásicos, podemos completar el proceso completo de optimización en una sola llamada. Aquí es donde entra en juego la Regresión por muestreo cuántico. Dado que el acceso a los ordenadores cuánticos sigue siendo un recurso escaso con alta demanda, este intercambio nos parece tanto factible como conveniente para muchos estudios actuales. Este enfoque aprovecha toda la capacidad clásica disponible y, al mismo tiempo, captura muchos de los mecanismos internos y propiedades intrínsecas de los cálculos cuánticos que no aparecen en la simulación.

La idea detrás de QSR es que la función de coste $C(\theta) := \langle \psi(\theta) | \hat{H} | \psi(\theta)\rangle$ puede expresarse como una serie de Fourier de la siguiente manera:

$$\begin{aligned} C(\vec{\theta}) & := \langle \psi(\theta) | \hat{H} | \psi(\theta)\rangle \\[1mm] & := a_0 + \sum_{k=1}^S[a_k\cos(k\theta)+ b_k\sin(k\theta)] \\[1mm] \end{aligned}$$

Dependiendo de la periodicidad y el ancho de banda de la función original, el conjunto $S$ puede ser finito o infinito. A efectos de esta discusión, asumiremos que es infinito. El siguiente paso es muestrear la función de coste $C(\theta)$ múltiples veces para obtener los coeficientes de Fourier $\{a_0, a_k, b_k\}_{k=1}^S$. En concreto, como tenemos $2S+1$ incógnitas, necesitaremos muestrear la función de coste $2S+1$ veces.

Si muestreamos la función de coste para $2S+1$ valores de parámetro $\{\theta_1,...,\theta_{2S+1}\}$, obtenemos el siguiente sistema:

$$\begin{pmatrix} 1 & \cos(\theta_1) & \sin(\theta_1) & \cos(2\theta_1) & ... & \sin(S\theta_1) \\ 1 & \cos(\theta_2) & \sin(\theta_2) & \cos(2\theta_2) & \cdots & \sin(S\theta_2)\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & \cos(\theta_{2S+1}) & \sin(\theta_{2S+1}) & \cos(2\theta_{2S+1}) & \cdots & \sin(S\theta_{2S+1}) \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ b_1 \\ a_2 \\ \vdots \\ b_S \end{pmatrix} = \begin{pmatrix} C(\theta_1) \\ C(\theta_2) \\ \vdots \\ C(\theta_{2S+1}) \end{pmatrix},$$

que reescribiremos como

$$Fa=c.$$

En la práctica, este sistema generalmente no es consistente porque los valores de la función de coste $c$ no son exactos. Por ello, suele ser buena idea normalizarlos multiplicando por $F^\dagger$ por la izquierda, lo que da como resultado:

$$F^\dagger Fa = F^\dagger c.$$

Este nuevo sistema siempre es consistente, y su solución es la solución de mínimos cuadrados al problema original. Si tenemos $k$ parámetros en lugar de solo uno, y cada parámetro $\theta^i$ tiene su propio $S_i$ para $i \in {1,...,k}$, entonces el número total de muestras requeridas es:

$$T=\prod_{i=1}^k(2S_i+1)\leq \prod_{i=1}^k(2S_{max}+1) = (2S_{max}+1)^n,$$

donde $S_{\max} = \max_i(S_i)$. Además, ajustar $S_{\max}$ como parámetro ajustable (en lugar de inferirlo) abre nuevas posibilidades, como:

• Sobremuestreo para mejorar la precisión.
• Submuestreo para aumentar el rendimiento reduciendo la sobrecarga de tiempo de ejecución o eliminando mínimos locales.

Esquema teórico

El esquema de QSR puede resumirse así:

• Preparar operadores de referencia $U_R$.
  • Iremos del estado $|0\rangle$ al estado de referencia $|\rho\rangle$
• Aplicar la forma variacional $U_V(\vec\theta_{i,j})$ para crear un ansatz $U_A(\vec\theta_{i,j})$.
  • Determinar el ancho de banda asociado a cada parámetro del ansatz. Una cota superior es suficiente.
• Arrancar en $i=0$ si se dispone de un problema similar (típicamente hallado mediante simulación clásica o muestreo).
• Muestrear la función de coste $C(\vec\theta) := a_0 + \sum_{k=1}^S[a_k\cos(k\theta)+ b_k\sin(k\theta)]$ al menos $T$ veces.
  • $T=\prod_{i=1}^k(2S_i+1)\leq \prod_{i=1}^k(2S_{max}+1) = (2S_{max}+1)^n$
  • Decidir si realizar sobremuestreo/submuestreo para equilibrar velocidad frente a precisión ajustando $T$.
• Calcular los coeficientes de Fourier a partir de las muestras (es decir, resolver el sistema de ecuaciones lineales normalizado).
• Resolver el mínimo global de la función de regresión resultante en una máquina clásica.

Resumen

Con esta lección, aprendiste sobre múltiples instancias variacionales disponibles:

• Esquema general
• Introducción de pesos y penalizaciones para ajustar una función de coste
• Exploración del submuestreo frente al sobremuestreo para intercambiar velocidad por precisión