อินสแตนซ์และส่วนขยาย

บทนี้จะครอบคลุมอัลกอริทึม variational เชิงควอนตัมหลายแบบ ได้แก่

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

ด้วยการใช้อัลกอริทึมเหล่านี้ เราจะเรียนรู้เกี่ยวกับแนวคิดการออกแบบหลายแบบที่สามารถนำไปใช้ในอัลกอริทึม variational แบบกำหนดเอง เช่น weights, penalties, over-sampling และ under-sampling เราสนับสนุนให้คุณทดลองกับแนวคิดเหล่านี้และแบ่งปันสิ่งที่ค้นพบกับชุมชน

Qiskit patterns framework ใช้กับอัลกอริทึมทั้งหมดเหล่านี้ แต่เราจะระบุขั้นตอนอย่างชัดเจนในตัวอย่างแรกเท่านั้น

Variational Quantum Eigensolver (VQE)

VQE เป็นหนึ่งในอัลกอริทึม variational เชิงควอนตัมที่ใช้กันมากที่สุด โดยเป็นต้นแบบสำหรับอัลกอริทึมอื่น ๆ

ขั้นตอนที่ 1: แมป input แบบคลาสสิกไปยังปัญหาเชิงควอนตัม

โครงร่างทางทฤษฎี

โครงร่างของ VQE นั้นเรียบง่าย:

• เตรียม reference operator $U_R$
  • เราเริ่มจากสถานะ $|0\rangle$ และไปยัง reference state $|\rho\rangle$
• ใช้ variational form $U_V(\vec\theta_{i,j})$ เพื่อสร้าง ansatz $U_A(\vec\theta_{i,j})$
  • เราไปจากสถานะ $|\rho\rangle$ ไปยัง $U_V(\vec\theta_{i,j})|\rho\rangle = |\psi(\vec\theta_{i,j})\rangle$
• Bootstrap ที่ $i=0$ ถ้ามีปัญหาที่คล้ายกัน (ปกติหาได้จากการจำลองแบบคลาสสิกหรือการสุ่มตัวอย่าง)
  • optimizer แต่ละตัวจะ bootstrap ต่างกัน ทำให้ได้ชุดเวกเตอร์พารามิเตอร์เริ่มต้น $\Theta_0 := \\{ {\vec\theta_{0,j} | j \in \mathcal{J}_\text{opt}^0} \\}$ (เช่น จากจุดเริ่มต้น $\vec\theta_0$)
• ประเมิน cost function $C(\vec\theta_{i,j}) := \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ สำหรับสถานะที่เตรียมไว้ทั้งหมดบนคอมพิวเตอร์ควอนตัม
• ใช้ classical optimizer เพื่อเลือกชุดพารามิเตอร์ถัดไป $\Theta_{i+1}$
• ทำซ้ำกระบวนการจนกว่าจะ converge

นี่คือ optimization loop แบบคลาสสิกง่าย ๆ ที่เราประเมิน cost function optimizer บางตัวอาจต้องการการประเมินหลายครั้งเพื่อคำนวณ gradient กำหนดการ iterate ถัดไป หรือประเมินการ converge

ต่อไปนี้เป็นตัวอย่างสำหรับ observable ดังต่อไปนี้:

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

การนำไปใช้

# 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()

เราสามารถใช้ cost function นี้เพื่อคำนวณพารามิเตอร์ที่เหมาะสมที่สุด

# 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

ขั้นตอนที่ 2: ปรับปรุงปัญหาสำหรับการรันบนควอนตัม

เราจะเลือก Backend ที่ว่างที่สุด และ import component ที่จำเป็นจาก 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')>

เราจะ transpile Circuit โดยใช้ preset pass manager ที่ optimization level 3 และใช้ layout ที่สอดคล้องกับ 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)

ขั้นตอนที่ 3: รันโดยใช้ Qiskit Runtime primitives

ตอนนี้เราพร้อมรันการคำนวณบนฮาร์ดแวร์ IBM Quantum® แล้ว เนื่องจากการ minimize cost function เป็นแบบ iterative สูง เราจะเริ่ม Runtime session เพื่อให้รอคิวเพียงครั้งเดียว เมื่อ job เริ่มรัน การ iterate แต่ละครั้งพร้อมการอัปเดตพารามิเตอร์จะรันทันที

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)

ขั้นตอนที่ 4: ประมวลผล ส่งคืนผลลัพธ์ในรูปแบบคลาสสิก

เราเห็นว่า minimization routine สำเร็จ หมายความว่าเราถึง tolerance เริ่มต้นของ COBYLA classical optimizer แล้ว ถ้าต้องการผลลัพธ์ที่แม่นยำกว่า เราสามารถระบุ tolerance ที่เล็กกว่าได้ ซึ่งอาจจำเป็น เนื่องจากผลลัพธ์ต่างจากที่ได้จาก simulator ข้างต้นหลายเปอร์เซ็นต์

ค่า x ที่ได้คือค่าประมาณที่ดีที่สุดในปัจจุบันของพารามิเตอร์ที่ minimize cost function ถ้าต้อง iterate เพื่อให้ได้ความแม่นยำสูงขึ้น ควรใช้ค่าเหล่านั้นแทน x0 ที่ใช้เริ่มต้น (เวกเตอร์ของค่า 1)

สุดท้าย เราสังเกตว่า function ถูกประเมิน 96 ครั้งในระหว่างการ optimize ซึ่งอาจต่างจากจำนวน optimization step เนื่องจาก optimizer บางตัวต้องการการประเมิน function หลายครั้งในขั้นตอนเดียว เช่น เมื่อประมาณ gradient

Subspace Search VQE (SSVQE)

SSVQE คือตัวแปรของ VQE ที่ช่วยให้ได้ $k$ eigenvalue แรกของ observable $\hat{H}$ ที่มี eigenvalue $\{\lambda_0, \lambda_1,...,\lambda_{N-1}\}$ โดยที่ $N\geq k$ โดยไม่เสียความทั่วไป เราสมมติว่า $\lambda_0<\lambda_1<...<\lambda_{N-1}$ SSVQE แนะนำแนวคิดใหม่โดยการเพิ่ม weights เพื่อช่วยจัดลำดับความสำคัญในการ optimize สำหรับ term ที่มี weight มากที่สุด

เพื่อนำอัลกอริทึมนี้ไปใช้ เราต้องการ reference state $\{ |\rho_j\rangle \}_{j=0}^{k-1}$ จำนวน $k$ ตัวที่ orthogonal กัน หมายความว่า $\langle \rho_j | \rho_l \rangle = \delta_{jl}$ สำหรับ $j,l<k$ สถานะเหล่านี้สามารถสร้างได้โดยใช้ Pauli operator จากนั้น cost function ของอัลกอริทึมนี้คือ:

$$\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}$$

โดยที่ $w_j$ คือจำนวนบวกใด ๆ ที่ถ้า $j<l<k$ แล้ว $w_j>w_l$ และ $U_V(\vec{\theta})$ คือ variational form ที่ผู้ใช้กำหนด

อัลกอริทึม SSVQE อาศัยข้อเท็จจริงที่ว่า eigenstate ที่สอดคล้องกับ eigenvalue ต่างกันนั้น orthogonal กัน โดยเฉพาะ inner product ของ $U_V(\vec{\theta})|\rho_j\rangle$ และ $U_V(\vec{\theta})|\rho_l\rangle$ สามารถเขียนได้เป็น:

$$\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}$$

ความเท่ากันแรกเป็นจริงเพราะ $U_{V}(\vec{\theta})$ เป็น quantum operator จึงเป็น unitary ความเท่ากันสุดท้ายเป็นจริงเพราะ orthogonality ของ reference state $|\rho_j\rangle$ ข้อเท็จจริงที่ว่า orthogonality ได้รับการรักษาผ่าน unitary transformation นั้นเกี่ยวข้องอย่างลึกซึ้งกับหลักการของการอนุรักษ์ข้อมูล ตามที่แสดงในศาสตร์ข้อมูลควอนตัม ในมุมมองนี้ การแปลงแบบ non-unitary แทนกระบวนการที่ข้อมูลถูกสูญหายหรือถูกฉีดเข้าไป

Weights $w_j$ ช่วยให้แน่ใจว่าสถานะทั้งหมดเป็น eigenstate ถ้า weights แตกต่างกันมากพอ term ที่มี weight มากที่สุด (นั่นคือ $w_0$) จะได้รับความสำคัญในระหว่างการ optimize มากกว่าตัวอื่น ส่งผลให้สถานะ $U_{V}(\vec{\theta})|\rho_0 \rangle$ จะกลายเป็น eigenstate ที่สอดคล้องกับ $\lambda_0$ เนื่องจาก $\{ U_{V}(\vec{\theta})|\rho_j\rangle \}_{j=0}^{k-1}$ orthogonal กัน สถานะที่เหลือจะ orthogonal กับมัน และด้วยเหตุนี้จึงอยู่ใน subspace ที่สอดคล้องกับ eigenvalue $\{\lambda_1,...,\lambda_{N-1}\}$

ด้วยการใช้ argument เดียวกันกับ term ที่เหลือ ลำดับความสำคัญถัดไปจะเป็น term ที่มี weight $w_1$ ดังนั้น $U_{V}(\vec{\theta})|\rho_1 \rangle$ จะเป็น eigenstate ที่สอดคล้องกับ $\lambda_1$ และ term อื่น ๆ จะอยู่ใน eigenspace ของ $\{\lambda_2,...,\lambda_{N-1}\}$

ด้วยการให้เหตุผลแบบ inductive เราสรุปได้ว่า $U_{V}(\vec{\theta})|\rho_j \rangle$ จะเป็น eigenstate โดยประมาณของ $\lambda_j$ สำหรับ $0\leq j < k$

โครงร่างทางทฤษฎี

SSVQE สรุปได้ดังนี้:

• เตรียม reference state หลายตัวโดยใช้ unitary U_R กับ k computational basis state ต่างกัน
  • อัลกอริทึมนี้ต้องใช้ reference state $\{ |\rho_j\rangle \}_{j=0}^{k-1}$ จำนวน $k$ ตัวที่ orthogonal กัน ซึ่ง $\langle \rho_j | \rho_l \rangle = \delta_{jl}$ สำหรับ $j,l<k$
• ใช้ variational form $U_V(\vec\theta_{i,j})$ กับแต่ละ reference state ได้ ansatz $U_A(\vec\theta_{i,j})$
• Bootstrap ที่ $i=0$ ถ้ามีปัญหาที่คล้ายกัน (ปกติหาได้จากการจำลองแบบคลาสสิกหรือการสุ่มตัวอย่าง)
• ประเมิน cost function $C(\vec\theta_{i,j}) := \sum_{j=0}^{k-1} w_j \langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ สำหรับสถานะที่เตรียมไว้ทั้งหมดบนคอมพิวเตอร์ควอนตัม
  • สามารถแยกออกเป็นการคำนวณ expectation value สำหรับ observable $\langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ และคูณผลลัพธ์ด้วย $w_j$
  • หลังจากนั้น cost function จะส่งคืนผลรวมของ weighted expectation value ทั้งหมด
• ใช้ classical optimizer เพื่อกำหนดชุดพารามิเตอร์ถัดไป $\Theta_{i+1}$
• ทำซ้ำขั้นตอนข้างต้นจนกว่าจะ converge

คุณจะสร้าง SSVQE cost function ในการประเมิน แต่เรามี snippet นี้เพื่อแรงบันดาลใจในการแก้ปัญหาของคุณ:

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

Variational Quantum Deflation (VQD)

VQD คือวิธีการ iterative ที่ขยาย VQE ให้ได้ $k$ eigenvalue แรกของ observable $\hat{H}$ ที่มี eigenvalue $\{\lambda_0, \lambda_1,...,\lambda_{N-1}\}$ โดยที่ $N\geq k$ แทนที่จะได้แค่ค่าแรก สำหรับส่วนที่เหลือ เราจะสมมติโดยไม่เสียความทั่วไปว่า $\lambda_0\leq\lambda_1\leq...\leq\lambda_{N-1}$ VQD แนะนำแนวคิดของ penalty cost เพื่อนำทางกระบวนการ optimization

VQD แนะนำ penalty term ที่แทนด้วย $\beta$ เพื่อสมดุลการมีส่วนร่วมของแต่ละ overlap term ต่อ cost Penalty term นี้ทำหน้าที่ลงโทษกระบวนการ optimization ถ้า orthogonality ไม่เกิดขึ้น เราบังคับใช้ข้อจำกัดนี้เพราะ eigenstate ของ observable หรือ Hermitian operator ที่สอดคล้องกับ eigenvalue ต่างกันนั้นเสมอ orthogonal กัน หรือสามารถทำให้เป็นเช่นนั้นได้ในกรณีของ degeneracy หรือ eigenvalue ซ้ำ ดังนั้น ด้วยการบังคับ orthogonality กับ eigenstate ที่สอดคล้องกับ $\lambda_0$ เราจึง optimize ใน subspace ที่สอดคล้องกับ eigenvalue ที่เหลือ $\{\lambda_1, \lambda_2,..., \lambda_{N-1}\}$ ซึ่ง $\lambda_1$ เป็น eigenvalue ต่ำสุดจาก eigenvalue ที่เหลือ และด้วยเหตุนี้ optimal solution ของปัญหาใหม่จึงหาได้โดยใช้ variational theorem

แนวคิดหลักของ VQD คือการใช้ VQE ตามปกติเพื่อหา lowest eigenvalue $\lambda_0 := C_0(\vec\theta^0) \equiv C_\text{VQE}(\vec\theta^0)$ พร้อมกับ (approximate) eigenstate $|\psi(\vec{\theta^0})\rangle$ ที่สอดคล้องกัน สำหรับ optimal parameter vector $\vec{\theta^0}$ บางค่า จากนั้น เพื่อหา eigenvalue ถัดไป $\lambda_1 > \lambda_0$ แทนที่จะ minimize cost function $C_0(\vec{\theta}) := \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ เราจะ optimize:

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

ค่าบวก $\beta_0$ ควรมีค่ามากกว่า $\lambda_1-\lambda_0$

สิ่งนี้แนะนำ cost function ใหม่ที่มองได้เป็นปัญหา constrained ซึ่งเรา minimize $C_\text{VQE}(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ โดยมีข้อจำกัดว่าสถานะต้อง orthogonal กับ $|\psi(\vec{\theta^0})\rangle$ ที่ได้มาก่อนหน้า โดย $\beta_0$ ทำหน้าที่เป็น penalty term ถ้าข้อจำกัดไม่ได้รับการปฏิบัติตาม

อีกทางหนึ่ง ปัญหาใหม่นี้สามารถตีความได้เป็นการรัน VQE บน 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,$$

สมมติว่า solution ของปัญหาใหม่คือ $|\psi(\vec{\theta^1})\rangle$ expectation value ของ $\hat{H}$ (ไม่ใช่ $\hat{H_1}$) ควรเป็น $\langle \psi(\vec{\theta^1}) | \hat{H} | \psi(\vec{\theta^1})\rangle = \lambda_1$ เพื่อหา eigenvalue ที่สาม $\lambda_2$ cost function ที่ต้อง optimize คือ:

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

โดยที่ $\beta_1$ คือค่าคงที่บวกที่มากพอเพื่อบังคับ orthogonality ของสถานะ solution กับทั้ง $|\psi(\vec{\theta^0})\rangle$ และ $|\psi(\vec{\theta^1})\rangle$ สิ่งนี้ลงโทษสถานะใน search space ที่ไม่ตรงตามเงื่อนไข ซึ่งจำกัด search space ได้อย่างมีประสิทธิภาพ ดังนั้น optimal solution ของปัญหาใหม่ควรเป็น eigenstate ที่สอดคล้องกับ $\lambda_2$

เช่นเดียวกับกรณีก่อนหน้า ปัญหาใหม่นี้ยังสามารถตีความเป็น VQE กับ 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.$$

ถ้า solution ของปัญหาใหม่นี้คือ $|\psi(\vec{\theta^2})\rangle$ expectation value ของ $\hat{H}$ (ไม่ใช่ $\hat{H_2}$) ควรเป็น $\langle \psi(\vec{\theta^2}) | \hat{H} | \psi(\vec{\theta^2})\rangle = \lambda_2$ ในทำนองเดียวกัน เพื่อหา $k$-th eigenvalue $\lambda_{k-1}$ คุณจะ minimize cost function:

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

โปรดจำไว้ว่าเราสร้าง $\vec{\theta^j}$ ให้ $\langle \psi(\vec{\theta^j}) | \hat{H} | \psi(\vec{\theta^j})\rangle = \lambda_j, \forall j<k$ ปัญหานี้เทียบเท่ากับการ minimize $C(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ แต่มีข้อจำกัดว่าสถานะต้อง orthogonal กับ $|\psi(\vec{\theta^j})\rangle ; \forall j \in {0, \cdots, k-1}$ ซึ่งจำกัด search space ให้อยู่ใน subspace ที่สอดคล้องกับ eigenvalue $\{\lambda_{k-1},\cdots,\lambda_{N-1}\}$

ปัญหานี้เทียบเท่ากับ VQE ที่มี 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,$$

ดังที่เห็นจากกระบวนการ เพื่อหา $k$-th eigenvalue คุณต้องการ (approximate) eigenstate ของ $k-1$ eigenvalue ก่อนหน้า ดังนั้นคุณจะต้องรัน VQE ทั้งหมด $k$ ครั้ง ดังนั้น cost function ของ VQD คือ:

$$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$$

โครงร่างทางทฤษฎี

โครงร่างของ VQD สรุปได้ดังนี้:

• เตรียม reference operator $U_R$
• ใช้ variational form $U_V(\vec\theta_{i,j})$ กับ reference state สร้าง ansatz $U_A(\vec\theta_{i,j})$
• Bootstrap ที่ $i=0$ ถ้ามีปัญหาที่คล้ายกัน (ปกติหาได้จากการจำลองแบบคลาสสิกหรือการสุ่มตัวอย่าง)
• ประเมิน cost function $C_k(\vec{\theta})$ ซึ่งเกี่ยวข้องกับการคำนวณ $k$ excited state และอาร์เรย์ของ $\beta$ ที่กำหนด overlap penalty สำหรับแต่ละ overlap term
  • คำนวณ expectation value สำหรับ observable $\langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ สำหรับแต่ละ $k$
  • คำนวณ penalty $\sum_{j=0}^{k-1}\beta_j |\langle \psi(\vec{\theta})| \psi(\vec{\theta^j})\rangle |^2$
  • cost function ควรส่งคืนผลรวมของสองส่วนนี้
• ใช้ classical optimizer เพื่อเลือกชุดพารามิเตอร์ถัดไป $\Theta_{i+1}$
• ทำซ้ำกระบวนการนี้จนกว่าจะ converge

การนำไปใช้

สำหรับการนำไปใช้นี้ เราจะสร้างฟังก์ชันสำหรับ overlap penalty ซึ่งจะใช้ใน cost function ในแต่ละ iteration กระบวนการนี้จะถูกทำซ้ำสำหรับแต่ละ excited state

from qiskit.circuit.library import TwoLocal

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

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

ก่อนอื่น เราจะตั้งค่าฟังก์ชันที่คำนวณ state fidelity ซึ่งเป็นเปอร์เซ็นต์ของ overlap ระหว่างสองสถานะที่จะใช้เป็น penalty สำหรับ 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)

ถึงเวลาเขียน VQD cost function แล้ว เช่นเดียวกับก่อนหน้าที่เราคำนวณแค่ ground state เราจะหา lowest energy state โดยใช้ Estimator primitive อย่างไรก็ตาม ดังที่อธิบายข้างต้น เราจะเพิ่ม penalty term เพื่อให้แน่ใจว่า higher-energy state นั้น orthogonal กัน นั่นคือ สำหรับแต่ละ excited state ใหม่ penalty จะถูกเพิ่มสำหรับ overlap ระหว่างสถานะ variational ปัจจุบันกับ lower-energy eigenstate ที่พบแล้ว

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

โปรดสังเกตโดยเฉพาะว่า cost function ข้างต้นอ้างอิงฟังก์ชัน calculate_overlaps ซึ่งสร้าง quantum circuit ใหม่จริง ๆ ถ้าต้องการรันบนฮาร์ดแวร์จริง circuit ใหม่นั้นต้องถูก transpile ด้วย ควรเป็นวิธีที่ optimal เพื่อรันบน Backend ที่เลือก โปรดทราบว่า transpilation ได้ถูก built in ไว้ในฟังก์ชัน calculate_overlaps หรือ cost_func_vqd แล้ว ลองปรับแต่งโค้ดเองเพื่อเพิ่ม transpilation เพิ่มเติม (และแบบมีเงื่อนไข) แต่สิ่งนี้จะถูกดำเนินการให้ในบทเรียนถัดไปด้วย

ในบทเรียนนี้ เราจะรันอัลกอริทึม VQD โดยใช้ Statevector Sampler และ Statevector Estimator:

from qiskit.primitives import StatevectorEstimator as Estimator

sampler = Sampler()
estimator = Estimator()

เราจะแนะนำ observable เพื่อประมาณ ในบทเรียนถัดไปเราจะเพิ่มบริบททางฟิสิกส์ให้กับมัน เช่น excited state ของโมเลกุล อาจมีประโยชน์ที่จะคิดว่า observable นี้เป็น Hamiltonian ของระบบที่มี excited state แม้ว่า observable นี้ไม่ได้ถูกเลือกให้ตรงกับโมเลกุลหรืออะตอมใดเฉพาะ

from qiskit.quantum_info import SparsePauliOp

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

ที่นี่ เราตั้งจำนวนสถานะทั้งหมดที่ต้องการคำนวณ (ground state และ excited state, k) และ penalties (betas) สำหรับ overlap ระหว่าง statevector ที่ควร orthogonal กัน ผลของการเลือก betas สูงหรือต่ำเกินไปจะถูกสำรวจในบทเรียนถัดไป สำหรับตอนนี้ เราจะใช้ค่าที่ให้ไว้ด้านล่าง เราจะเริ่มด้วยค่า zeros ทั้งหมดเป็นพารามิเตอร์ ในการคำนวณของคุณเอง คุณอาจต้องการใช้พารามิเตอร์เริ่มต้นที่ฉลาดกว่านี้โดยอาศัยความรู้เกี่ยวกับระบบหรือการคำนวณก่อนหน้า

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

ตอนนี้เราสามารถรันการคำนวณได้:

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

ค่าที่ได้จาก cost function คือประมาณ -6.00, 4.02 และ 5.61 สิ่งสำคัญเกี่ยวกับผลลัพธ์เหล่านี้คือ function value มีค่าเพิ่มขึ้น ถ้าเราได้ first excited state ที่มีพลังงานต่ำกว่าการคำนวณ ground state เริ่มต้นที่ไม่มีข้อจำกัด นั่นจะบ่งชี้ว่ามีข้อผิดพลาดในโค้ดของเรา

ค่า x คือพารามิเตอร์ที่ให้ statevector ที่สอดคล้องกับแต่ละ cost (พลังงาน)

สุดท้าย เราสังเกตว่า minimization ทั้งสามครั้ง converge ภายใน tolerance เริ่มต้นของ classical optimizer (COBYLA ในที่นี้) ต้องการการประเมิน function 131, 110 และ 90 ครั้งตามลำดับ

Quantum Sampling Regression (QSR)

ปัญหาหลักของ VQE คือการเรียก quantum computer หลายครั้งที่จำเป็นเพื่อหาพารามิเตอร์สำหรับแต่ละขั้นตอน เช่น $k$, $k-1$ และต่อไปเรื่อย ๆ สิ่งนี้เป็นปัญหาอย่างยิ่งเมื่อการเข้าถึงอุปกรณ์ควอนตัมต้องรอคิว แม้ว่า Session สามารถใช้จัดกลุ่มการเรียกแบบ iterative หลาย ๆ ครั้ง แต่ทางเลือกอื่นคือการใช้การสุ่มตัวอย่าง ด้วยการใช้ทรัพยากรคลาสสิกมากขึ้น เราสามารถทำกระบวนการ optimization ทั้งหมดให้เสร็จในการเรียกครั้งเดียว นี่คือจุดที่ Quantum Sampling Regression เข้ามา เนื่องจากการเข้าถึงคอมพิวเตอร์ควอนตัมยังคงเป็น commodity ที่มีอุปทานน้อย/อุปสงค์สูง เราจึงเห็นว่า trade-off นี้เป็นไปได้และสะดวกสำหรับการศึกษาปัจจุบันหลายแบบ วิธีนี้ใช้ประโยชน์จากความสามารถคลาสสิกที่มีอยู่ทั้งหมด ขณะที่ยังคงจับเอาการทำงานภายในและคุณสมบัติเนื้อในของการคำนวณเชิงควอนตัมที่ไม่ปรากฏในการจำลอง

แนวคิดของ QSR คือ cost function $C(\theta) := \langle \psi(\theta) | \hat{H} | \psi(\theta)\rangle$ สามารถแสดงเป็น Fourier series ในลักษณะดังนี้:

$$\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}$$

ขึ้นอยู่กับ periodicity และ bandwidth ของฟังก์ชันดั้งเดิม ชุด $S$ อาจมีจำกัดหรือไม่มีจำกัด สำหรับการอภิปรายนี้ เราจะสมมติว่ามันไม่มีจำกัด ขั้นตอนถัดไปคือการสุ่มตัวอย่าง cost function $C(\theta)$ หลายครั้งเพื่อหา Fourier coefficient $\{a_0, a_k, b_k\}_{k=1}^S$ โดยเฉพาะ เนื่องจากเรามี $2S+1$ unknown เราต้องสุ่มตัวอย่าง cost function $2S+1$ ครั้ง

ถ้าเราสุ่มตัวอย่าง cost function สำหรับ $2S+1$ ค่าพารามิเตอร์ $\{\theta_1,...,\theta_{2S+1}\}$ เราสามารถหาระบบดังต่อไปนี้:

$$\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},$$

ซึ่งเราจะเขียนใหม่เป็น

$$Fa=c.$$

ในทางปฏิบัติ ระบบนี้มักไม่ consistent เพราะค่า cost function $c$ ไม่แน่นอน ดังนั้น มักเป็นความคิดที่ดีที่จะ normalize โดยการคูณด้วย $F^\dagger$ ทางซ้าย ซึ่งให้ผล:

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

ระบบใหม่นี้ consistent เสมอ และ solution คือ least-squares solution ของปัญหาเดิม ถ้าเรามี $k$ พารามิเตอร์แทนที่จะเป็นหนึ่งตัว และแต่ละพารามิเตอร์ $\theta^i$ มี $S_i$ ของตัวเองสำหรับ $i \in {1,...,k}$ จำนวนตัวอย่างทั้งหมดที่จำเป็นคือ:

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

โดยที่ $S_{\max} = \max_i(S_i)$ นอกจากนี้ การปรับ $S_{\max}$ เป็นพารามิเตอร์ที่ปรับได้ (แทนที่จะอนุมาน) เปิดความเป็นไปได้ใหม่ เช่น:

• Over-sampling เพื่อเพิ่มความแม่นยำ
• Under-sampling เพื่อเพิ่มประสิทธิภาพโดยลด runtime overhead หรือกำจัด local minima

โครงร่างทางทฤษฎี

โครงร่างของ QSR สรุปได้ดังนี้:

• เตรียม reference operator $U_R$
  • เราไปจากสถานะ $|0\rangle$ ไปยัง reference state $|\rho\rangle$
• ใช้ variational form $U_V(\vec\theta_{i,j})$ เพื่อสร้าง ansatz $U_A(\vec\theta_{i,j})$
  • กำหนด bandwidth ที่เชื่อมโยงกับแต่ละพารามิเตอร์ใน ansatz ค่าขอบเขตบนก็เพียงพอ
• Bootstrap ที่ $i=0$ ถ้ามีปัญหาที่คล้ายกัน (ปกติหาได้จากการจำลองแบบคลาสสิกหรือการสุ่มตัวอย่าง)
• สุ่มตัวอย่าง cost function $C(\vec\theta) := a_0 + \sum_{k=1}^S[a_k\cos(k\theta)+ b_k\sin(k\theta)]$ อย่างน้อย $T$ ครั้ง
  • $T=\prod_{i=1}^k(2S_i+1)\leq \prod_{i=1}^k(2S_{max}+1) = (2S_{max}+1)^n$
  • ตัดสินใจว่าจะ over-sample/under-sample เพื่อสมดุลระหว่างความเร็วและความแม่นยำโดยปรับ $T$
• คำนวณ Fourier coefficient จากตัวอย่าง (นั่นคือ แก้ระบบสมการเชิงเส้นที่ normalize แล้ว)
• แก้หา global minimum ของ regression function ที่ได้บนเครื่องคลาสสิก

สรุป

จากบทเรียนนี้ คุณได้เรียนรู้เกี่ยวกับ variational instance ต่าง ๆ ที่มีอยู่:

• โครงร่างทั่วไป
• การแนะนำ weights และ penalties เพื่อปรับ cost function
• การสำรวจ under-sampling และ over-sampling เพื่อ trade-off ระหว่างความเร็วและความแม่นยำ