ฟังก์ชันค่าใช้จ่าย

ในบทเรียนนี้ เราจะเรียนรู้วิธีประเมิน ฟังก์ชันค่าใช้จ่าย:

• ขั้นแรก เราจะทำความรู้จักกับ Qiskit Runtime primitives
• นิยาม ฟังก์ชันค่าใช้จ่าย $C(\vec\theta)$ ซึ่งเป็นฟังก์ชันเฉพาะสำหรับปัญหาที่กำหนดเป้าหมายให้ optimizer นำไปลดค่า (หรือเพิ่มค่า)
• กำหนดกลยุทธ์การวัดด้วย Qiskit Runtime primitives เพื่อสมดุลระหว่างความเร็วและความแม่นยำ

 

Primitives

ระบบฟิสิกส์ทุกอย่าง ไม่ว่าจะเป็นคลาสสิกหรือควอนตัม ล้วนสามารถอยู่ในสถานะต่างๆ ได้ ตัวอย่างเช่น รถบนถนนสามารถมีมวล ตำแหน่ง ความเร็ว หรือความเร่งที่บ่งบอกสถานะได้ เช่นเดียวกัน ระบบควอนตัมก็สามารถมีการจัดเรียงหรือสถานะที่แตกต่างกันได้ แต่ต่างจากระบบคลาสสิกในวิธีที่เราจัดการกับการวัดและการเปลี่ยนแปลงสถานะ สิ่งนี้นำไปสู่คุณสมบัติเฉพาะตัว เช่น superposition และ entanglement ที่มีเฉพาะในกลศาสตร์ควอนตัม เช่นเดียวกับที่เราอธิบายสถานะของรถด้วยคุณสมบัติทางฟิสิกส์อย่างความเร็วหรือความเร่ง เราก็สามารถอธิบายสถานะของระบบควอนตัมด้วย ออบเซอร์วาเบิล ซึ่งเป็นวัตถุทางคณิตศาสตร์

ในกลศาสตร์ควอนตัม สถานะถูกแทนด้วยเวกเตอร์คอลัมน์จำนวนเชิงซ้อนที่ normalize แล้ว หรือ kets ($|\psi\rangle$) และออบเซอร์วาเบิลคือตัวดำเนินการเชิงเส้นแบบ Hermitian ($\hat{H}=\hat{H}^{\dagger}$) ที่กระทำต่อ kets ไอเกนเวกเตอร์ ($|\lambda\rangle$) ของออบเซอร์วาเบิลเรียกว่า eigenstate การวัดออบเซอร์วาเบิลสำหรับ eigenstate หนึ่งๆ ($|\lambda\rangle$) จะให้ค่าไอเกนวาลู ($\lambda$) ที่สอดคล้องกันเป็นผลลัพธ์

หากสงสัยว่าจะวัดระบบควอนตัมอย่างไรและวัดอะไรได้บ้าง Qiskit มี primitives สองตัวที่ช่วยได้:

• Sampler: กำหนดสถานะควอนตัม $|\psi\rangle$ primitive นี้จะหาความน่าจะเป็นของแต่ละสถานะฐานการคำนวณที่เป็นไปได้
• Estimator: กำหนดออบเซอร์วาเบิลควอนตัม $\hat{H}$ และสถานะ $|\psi\rangle$ primitive นี้จะคำนวณค่าคาดหวังของ $\hat{H}$

The Sampler primitive

Sampler primitive คำนวณความน่าจะเป็นของการได้สถานะ $|k\rangle$ แต่ละสถานะจากฐานการคำนวณ โดยกำหนด quantum circuit ที่เตรียมสถานะ $|\psi\rangle$ โดยคำนวณ

$$p_k = |\langle k | \psi \rangle|^2 \quad \forall k \in \mathbb{Z}_2^n \equiv \{0,1,\cdots,2^n-1\},$$

โดยที่ $n$ คือจำนวน Qubit และ $k$ คือการแทนค่าจำนวนเต็มของ binary string ผลลัพธ์ที่เป็นไปได้ใดๆ $\{0,1\}^n$ (นั่นคือ จำนวนเต็มฐาน $2$)

Qiskit Runtime Sampler รัน Circuit หลายครั้งบนอุปกรณ์ควอนตัม ทำการวัดในแต่ละรอบ และสร้างการกระจายความน่าจะเป็นขึ้นใหม่จาก bitstrings ที่ได้ ยิ่งรันมากครั้ง (หรือ shots) ผลลัพธ์ก็ยิ่งแม่นยำ แต่ต้องใช้เวลาและทรัพยากรควอนตัมมากขึ้น

อย่างไรก็ตาม เนื่องจากจำนวนผลลัพธ์ที่เป็นไปได้เติบโตแบบ exponential ตามจำนวน Qubit $n$ (คือ $2^n$) จำนวน shots ก็จะต้องเติบโตแบบ exponential เช่นกันเพื่อให้ได้การกระจายความน่าจะเป็นแบบ dense ดังนั้น Sampler จึงมีประสิทธิภาพเฉพาะสำหรับการกระจายความน่าจะเป็นแบบ sparse โดยสถานะเป้าหมาย $|\psi\rangle$ ต้องสามารถแสดงเป็นการรวมเชิงเส้นของสถานะฐานการคำนวณ โดยมีจำนวนพจน์เติบโตได้มากที่สุดแบบ polynomial ตามจำนวน Qubit:

$$|\psi\rangle = \sum^{\text{Poly}(n)}_k w_k |k\rangle.$$

Sampler สามารถกำหนดค่าให้ดึงความน่าจะเป็นจากส่วนย่อยของ Circuit ซึ่งแทนชุดย่อยของสถานะทั้งหมดที่เป็นไปได้ได้ด้วย

The Estimator primitive

Estimator primitive คำนวณค่าคาดหวังของออบเซอร์วาเบิล $\hat{H}$ สำหรับสถานะควอนตัม $|\psi\rangle$ โดยที่ความน่าจะเป็นของออบเซอร์วาเบิลสามารถแสดงเป็น $p_\lambda = |\langle\lambda|\psi\rangle|^2$ โดย $|\lambda\rangle$ คือ eigenstate ของออบเซอร์วาเบิล $\hat{H}$ ค่าคาดหวังถูกนิยามเป็นค่าเฉลี่ยของผลลัพธ์ที่เป็นไปได้ทั้งหมด $\lambda$ (นั่นคือ ไอเกนวาลูของออบเซอร์วาเบิล) จากการวัดสถานะ $|\psi\rangle$ ถ่วงน้ำหนักด้วยความน่าจะเป็นที่สอดคล้องกัน:

$$\langle\hat{H}\rangle_\psi := \sum_\lambda p_\lambda \lambda = \langle \psi | \hat{H} | \psi \rangle$$

อย่างไรก็ตาม การคำนวณค่าคาดหวังของออบเซอร์วาเบิลไม่สามารถทำได้เสมอไป เพราะเรา종종ไม่รู้ eigenbasis ของมัน Qiskit Runtime Estimator ใช้กระบวนการพีชคณิตที่ซับซ้อนในการประมาณค่าคาดหวังบนอุปกรณ์ควอนตัมจริง โดยแยกออบเซอร์วาเบิลออกเป็นออบเซอร์วาเบิลอื่นๆ ที่รู้ eigenbasis

พูดง่ายๆ คือ Estimator แยกออบเซอร์วาเบิลใดๆ ที่ไม่รู้วิธีวัดออกเป็นออบเซอร์วาเบิลที่วัดได้ง่ายกว่าเรียกว่า Pauli operators ตัวดำเนินการใดๆ สามารถแสดงเป็นการรวมกันของ Pauli operators จำนวน $4^n$ ตัว

$$\hat{P}_k := \sigma_{k_{n-1}}\otimes \cdots \otimes \sigma_{k_0} \quad \forall k \in \mathbb{Z}_4^n \equiv \{0,1,\cdots,4^n-1\}, \\$$

ซึ่งทำให้

$$\hat{H} = \sum^{4^n-1}_{k=0} w_k \hat{P}_k$$

โดยที่ $n$ คือจำนวน Qubit, $k \equiv k_{n-1} \cdots k_0$ สำหรับ $k_l \in \mathbb{Z}_4 \equiv \{0, 1, 2, 3\}$ (นั่นคือ จำนวนเต็มฐาน $4$) และ $(\sigma_0, \sigma_1, \sigma_2, \sigma_3) := (I, X, Y, Z)$

หลังจากทำการแยกนี้แล้ว Estimator จะสร้าง Circuit ใหม่ $V_k|\psi\rangle$ สำหรับแต่ละออบเซอร์วาเบิล $\hat{P}_k$ (จาก Circuit เดิม) เพื่อ diagonalize Pauli observable ในฐานการคำนวณและวัดมัน เราสามารถวัด Pauli observables ได้ง่ายเพราะรู้ $V_k$ ล่วงหน้า ซึ่งโดยทั่วไปไม่เป็นเช่นนั้นสำหรับออบเซอร์วาเบิลอื่นๆ

สำหรับแต่ละ $\hat{P}_{k}$ นั้น Estimator จะรัน Circuit ที่สอดคล้องกันบนอุปกรณ์ควอนตัมหลายครั้ง วัดสถานะผลลัพธ์ในฐานการคำนวณ และคำนวณความน่าจะเป็น $p_{kj}$ ของการได้ผลลัพธ์ที่เป็นไปได้แต่ละตัว $j$ จากนั้นค้นหาไอเกนวาลู $\lambda_{kj}$ ของ $P_k$ ที่สอดคล้องกับผลลัพธ์แต่ละตัว $j$ คูณด้วย $w_k$ และรวมผลลัพธ์ทั้งหมดเข้าด้วยกันเพื่อให้ได้ค่าคาดหวังของออบเซอร์วาเบิล $\hat{H}$ สำหรับสถานะ $|\psi\rangle$ ที่กำหนด

$$\langle\hat{H}\rangle_\psi = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}p_{kj} \lambda_{kj},$$

เนื่องจากการคำนวณค่าคาดหวังของ Paulis จำนวน $4^n$ ตัวนั้นไม่ปฏิบัติได้จริง (คือ เติบโตแบบ exponential) Estimator จึงมีประสิทธิภาพก็ต่อเมื่อ $w_k$ จำนวนมากเป็นศูนย์ (คือ การแยก Pauli แบบ sparse แทนที่จะเป็น dense) โดยในทางการเราบอกว่า สำหรับการคำนวณนี้จะ แก้ได้อย่างมีประสิทธิภาพ จำนวนพจน์ที่ไม่เป็นศูนย์ต้องเติบโตได้มากที่สุดแบบ polynomial ตามจำนวน Qubit $n$: $\hat{H} = \sum^{\text{Poly}(n)}_k w_k \hat{P}_k.$

ผู้อ่านอาจสังเกตสมมติฐานโดยนัยว่า การสุ่มตัวอย่างความน่าจะเป็นก็ต้องมีประสิทธิภาพด้วย ตามที่อธิบายสำหรับ Sampler ซึ่งหมายความว่า

$$\langle\hat{H}\rangle_\psi = \sum_{k}^{\text{Poly}(n)} w_k \sum_{j}^{\text{Poly}(n)}p_{kj} \lambda_{kj}.$$

ตัวอย่างนำทางสำหรับการคำนวณค่าคาดหวัง

สมมติสถานะ single-qubit $|+\rangle := H|0\rangle = \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)$ และออบเซอร์วาเบิล

$$\begin{aligned} \hat{H} & = \begin{pmatrix} -1 & 2 \\ 2 & 1 \\ \end{pmatrix}\\[1mm] & = 2X - Z \end{aligned}$$

โดยมีค่าคาดหวังตามทฤษฎี $\langle\hat{H}\rangle_+ = \langle+|\hat{H}|+\rangle = 2.$

เนื่องจากเราไม่รู้วิธีวัดออบเซอร์วาเบิลนี้ เราจึงไม่สามารถคำนวณค่าคาดหวังได้โดยตรง และต้องแสดงใหม่เป็น $\langle\hat{H}\rangle_+ = 2\langle X \rangle_+ - \langle Z \rangle_+$ ซึ่งสามารถแสดงให้ได้ผลลัพธ์เดิมโดยสังเกตว่า $\langle+|X|+\rangle = 1$ และ $\langle+|Z|+\rangle = 0$

มาดูวิธีคำนวณ $\langle X \rangle_+$ และ $\langle Z \rangle_+$ โดยตรงกัน เนื่องจาก $X$ และ $Z$ ไม่ commute (นั่นคือ ไม่มี eigenbasis ร่วมกัน) จึงวัดพร้อมกันไม่ได้ เราจึงต้องใช้ auxiliary circuits:

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

from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

# The following code will work for any other initial single-qubit state and observable
original_circuit = QuantumCircuit(1)
original_circuit.h(0)

H = SparsePauliOp(["X", "Z"], [2, -1])

aux_circuits = []
for pauli in H.paulis:
    aux_circ = original_circuit.copy()
    aux_circ.barrier()
    if str(pauli) == "X":
        aux_circ.h(0)
    elif str(pauli) == "Y":
        aux_circ.sdg(0)
        aux_circ.h(0)
    else:
        aux_circ.id(0)
    aux_circ.measure_all()
    aux_circuits.append(aux_circ)

original_circuit.draw("mpl")

# Auxiliary circuit for X
aux_circuits[0].draw("mpl")

# Auxiliary circuit for Z
aux_circuits[1].draw("mpl")

ตอนนี้เราสามารถทำการคำนวณด้วยตนเองโดยใช้ Sampler และตรวจสอบผลลัพธ์บน Estimator ได้:

from qiskit.primitives import StatevectorSampler, StatevectorEstimator
from qiskit.result import QuasiDistribution
import numpy as np

## SAMPLER
shots = 10000
sampler = StatevectorSampler()
job = sampler.run(aux_circuits, shots=shots)

# Run the sampler job and step through results
expvals = []
for index, pauli in enumerate(H.paulis):
    data_pub = job.result()[index].data
    bitstrings = data_pub.meas.get_bitstrings()
    counts = data_pub.meas.get_counts()
    quasi_dist = QuasiDistribution(
        {outcome: freq / shots for outcome, freq in counts.items()}
    )

    # Use the probabilities and known eigenvalues of Pauli operators to estimate the expectation value.
    val = 0

    if str(pauli) == "X":
        val += -1 * quasi_dist.get(1, 0)
        val += 1 * quasi_dist.get(0, 0)

    if str(pauli) == "Y":
        val += -1 * quasi_dist.get(1, 0)
        val += 1 * quasi_dist.get(0, 0)

    if str(pauli) == "Z":
        val += 1 * quasi_dist.get(0, 0)
        val += -1 * quasi_dist.get(1, 0)

    expvals.append(val)

# Print expectation values

print("Sampler results:")
for pauli, expval in zip(H.paulis, expvals):
    print(f"  >> Expected value of {str(pauli)}: {expval:.5f}")

total_expval = np.sum(H.coeffs * expvals).real
print(f"  >> Total expected value: {total_expval:.5f}")

# Use estimator for comparison
observables = [
    *H.paulis,
    H,
]  # Note: run for individual Paulis as well as full observable H

estimator = StatevectorEstimator()
job = estimator.run([(original_circuit, observables)])
estimator_expvals = job.result()[0].data.evs

# Print results
print("Estimator results:")
for obs, expval in zip(observables, estimator_expvals):
    if obs is not H:
        print(f"  >> Expected value of {str(obs)}: {expval:.5f}")
    else:
        print(f"  >> Total expected value: {expval:.5f}")

Sampler results:
  >> Expected value of X: 1.00000
  >> Expected value of Z: 0.00420
  >> Total expected value: 1.99580
Estimator results:
  >> Expected value of X: 1.00000
  >> Expected value of Z: 0.00000
  >> Total expected value: 2.00000

ความเข้มงวดทางคณิตศาสตร์ (ทางเลือก)

การแสดง $|\psi\rangle$ ในรูปของ eigenstates ของ $\hat{H}$ คือ $|\psi\rangle = \sum_\lambda a_\lambda |\lambda\rangle$ จะได้:

$$\begin{aligned} \langle \psi | \hat{H} | \psi \rangle & = \bigg(\sum_{\lambda'}a^*_{\lambda'} \langle \lambda'|\bigg) \hat{H} \bigg(\sum_{\lambda} a_\lambda | \lambda\rangle\bigg)\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \langle \lambda'|\hat{H}| \lambda\rangle\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \lambda \langle \lambda'| \lambda\rangle\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \lambda \cdot \delta_{\lambda, \lambda'}\\[1mm] & = \sum_\lambda |a_\lambda|^2 \lambda\\[1mm] & = \sum_\lambda p_\lambda \lambda\\[1mm] \end{aligned}$$

เนื่องจากเราไม่รู้ไอเกนวาลูหรือ eigenstates ของออบเซอร์วาเบิลเป้าหมาย $\hat{H}$ เราจึงต้องพิจารณา diagonalization ของมันก่อน เนื่องจาก $\hat{H}$ เป็น Hermitian จึงมีการแปลง unitary $V$ ซึ่ง $\hat{H}=V^\dagger \Lambda V$ โดยที่ $\Lambda$ คือเมทริกซ์ diagonal ของไอเกนวาลู ดังนั้น $\langle j | \Lambda | k \rangle = 0$ ถ้า $j\neq k$ และ $\langle j | \Lambda | j \rangle = \lambda_j$

ซึ่งหมายความว่าค่าคาดหวังสามารถเขียนใหม่ได้เป็น:

$$\begin{aligned} \langle\psi|\hat{H}|\psi\rangle & = \langle\psi|V^\dagger \Lambda V|\psi\rangle\\[1mm] & = \langle\psi|V^\dagger \bigg(\sum_{j=0}^{2^n-1} |j\rangle \langle j|\bigg) \Lambda \bigg(\sum_{k=0}^{2^n-1} |k\rangle \langle k|\bigg) V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1} \sum_{k=0}^{2^n-1}\langle\psi|V^\dagger |j\rangle \langle j| \Lambda |k\rangle \langle k| V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1}\langle\psi|V^\dagger |j\rangle \langle j| \Lambda |j\rangle \langle j| V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1}|\langle j| V|\psi\rangle|^2 \lambda_j\\[1mm] \end{aligned}$$

เนื่องจากถ้าระบบอยู่ในสถานะ $|\phi\rangle = V |\psi\rangle$ ความน่าจะเป็นของการวัดได้ $| j\rangle$ คือ $p_j = |\langle j|\phi \rangle|^2$ ค่าคาดหวังข้างต้นสามารถแสดงได้เป็น:

$$\langle\psi|\hat{H}|\psi\rangle = \sum_{j=0}^{2^n-1} p_j \lambda_j.$$

สิ่งสำคัญมากที่ต้องสังเกตคือ ความน่าจะเป็นถูกนำมาจากสถานะ $V |\psi\rangle$ แทนที่จะเป็น $|\psi\rangle$ นี่คือเหตุผลที่เมทริกซ์ $V$ มีความจำเป็นอย่างยิ่ง คุณอาจสงสัยว่าจะหาเมทริกซ์ $V$ และไอเกนวาลู $\Lambda$ ได้อย่างไร ถ้าคุณมีไอเกนวาลูอยู่แล้ว ก็ไม่จำเป็นต้องใช้คอมพิวเตอร์ควอนตัมเลย เพราะเป้าหมายของอัลกอริทึมแบบ variational คือการหาไอเกนวาลูเหล่านี้ของ $\hat{H}$

โชคดีที่มีวิธีแก้: เมทริกซ์ขนาด $2^n \times 2^n$ ใดๆ สามารถเขียนเป็นการรวมเชิงเส้นของ tensor products ของ Pauli matrices และ identities จำนวน $4^n$ ตัว ทั้งหมดล้วน Hermitian และ unitary โดยมี $V$ และ $\Lambda$ ที่รู้จัก นี่คือสิ่งที่ Runtime's Estimator ทำภายในโดยแยก object Operator ใดๆ ออกเป็น SparsePauliOp

นี่คือ Operators ที่สามารถใช้ได้:

$$\begin{array}{c|c|c|c} \text{Operator} & \sigma & V & \Lambda \\[1mm] \hline I & \sigma_0 = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} & V_0 = I & \Lambda_0 = I = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \\[4mm] X & \sigma_1 = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} & V_1 = H =\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} & \Lambda_1 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \\[4mm] Y & \sigma_2 = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} & V_2 = HS^\dagger =\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}\cdot \begin{pmatrix} 1 & 0 \\ 0 & -i \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & -i \\ 1 & i \end{pmatrix}\quad & \Lambda_2 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \\[4mm] Z & \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} & V_3 = I & \Lambda_3 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \end{array}$$

ดังนั้นมาเขียน $\hat{H}$ ใหม่ในรูปของ Paulis และ identities:

$$\hat{H} = \sum_{k_{n-1}=0}^3... \sum_{k_0=0}^3 w_{k_{n-1}...k_0} \sigma_{k_{n-1}}\otimes ... \otimes \sigma_{k_0} = \sum_{k=0}^{4^n-1} w_k \hat{P}_k,$$

โดยที่ $k = \sum_{l=0}^{n-1} 4^l k_l \equiv k_{n-1}...k_0$ สำหรับ $k_{n-1},...,k_0\in \{0,1,2,3\}$ (คือ ฐาน $4$) และ $\hat{P}_{k} := \sigma_{k_{n-1}}\otimes ... \otimes \sigma_{k_0}$:

$$\begin{aligned} \langle\psi|\hat{H}|\psi\rangle & = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}|\langle j| V_k|\psi\rangle|^2 \langle j| \Lambda_k |j\rangle \\[1mm] & = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}p_{kj} \lambda_{kj}, \\[1mm] \end{aligned}$$

โดยที่ $V_k := V_{k_{n-1}}\otimes ... \otimes V_{k_0}$ และ $\Lambda_k := \Lambda_{k_{n-1}}\otimes ... \otimes \Lambda_{k_0}$ ซึ่งทำให้: $\hat{P_k}=V_k^\dagger \Lambda_k V_k.$

Cost functions

โดยทั่วไปแล้ว cost function ถูกใช้เพื่ออธิบายเป้าหมายของปัญหาและวัดว่า trial state ทำงานได้ดีแค่ไหนเมื่อเทียบกับเป้าหมายนั้น นิยามนี้สามารถนำไปใช้กับตัวอย่างต่าง ๆ ในสาขาเคมี, machine learning, การเงิน, การปรับแต่ง (optimization) และอื่น ๆ

มาดูตัวอย่างง่าย ๆ ของการหา ground state ของระบบกัน เป้าหมายของเราคือลด expected value ของ observable ที่แทนพลังงาน (Hamiltonian $\hat{\mathcal{H}}$) ให้ต่ำที่สุด:

$$\min_{\vec\theta} \langle\psi(\vec\theta)|\hat{\mathcal{H}}|\psi(\vec\theta)\rangle$$

เราสามารถใช้ Estimator เพื่อคำนวณ expected value แล้วส่งค่านี้ไปให้ optimizer เพื่อทำการ minimize ถ้า optimization สำเร็จ มันจะคืนค่าพารามิเตอร์ที่เหมาะสมที่สุด $\vec\theta^*$ ซึ่งเราจะสามารถสร้าง solution state $|\psi(\vec\theta^*)\rangle$ ที่เสนอมาได้ และคำนวณ observed expected value เป็น $C(\vec\theta^*)$

สังเกตว่าเราจะสามารถ minimize cost function ได้เฉพาะสำหรับชุดของ state ที่จำกัดที่เรากำลังพิจารณาเท่านั้น ซึ่งนำไปสู่ความเป็นไปได้สองแบบ:

• ansatz ของเราไม่ได้กำหนด solution state ครอบคลุม search space: ถ้าเป็นเช่นนี้ optimizer จะไม่มีทางหาคำตอบได้ และเราต้องทดลองใช้ ansatz อื่น ๆ ที่อาจจะแทนค่า search space ของเราได้แม่นยำกว่า
• optimizer ของเราไม่สามารถหาคำตอบที่ถูกต้องนี้ได้: Optimization สามารถนิยามแบบ globally และ locally ได้ เราจะสำรวจว่าหมายความว่าอย่างไรในหัวข้อถัดไป

โดยรวมแล้ว เราจะทำ classical optimization loop แต่พึ่งพาการประเมิน cost function บนคอมพิวเตอร์ควอนตัม จากมุมมองนี้ อาจมองว่า optimization เป็นงาน classical ล้วน ๆ ที่เราเรียก black-box quantum oracle ทุกครั้งที่ optimizer ต้องการประเมิน cost function

def cost_func_vqe(params, circuit, 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
    """
    pub = (circuit, hamiltonian, params)
    cost = estimator.run([pub]).result()[0].data.evs
    return cost

from qiskit.circuit.library import TwoLocal

observable = SparsePauliOp.from_list([("XX", 1), ("YY", -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)

theta_list = (2 * np.pi * np.random.rand(1, 8)).tolist()
ansatz.decompose().draw("mpl")

ขั้นแรกเราจะทำสิ่งนี้โดยใช้ simulator: StatevectorEstimator ซึ่งมักแนะนำให้ใช้สำหรับ debugging แต่เราจะทำการ debug แล้วตามด้วยการคำนวณบนฮาร์ดแวร์ควอนตัมจริงทันที ในปัจจุบัน ปัญหาที่น่าสนใจมากขึ้นเรื่อย ๆ ไม่สามารถ simulate แบบ classical ได้หากไม่มีสิ่งอำนวยความสะดวก supercomputing ระดับสูง

estimator = StatevectorEstimator()
cost = cost_func_vqe(theta_list, ansatz, observable, estimator)
print(cost)

[-0.58744589]

ตอนนี้เราจะดำเนินการรันบนคอมพิวเตอร์ควอนตัมจริง สังเกต syntax ที่เปลี่ยนไป ขั้นตอนที่เกี่ยวข้องกับ pass_manager จะถูกพูดถึงเพิ่มเติมในตัวอย่างถัดไป หนึ่งขั้นตอนที่สำคัญอย่างยิ่งใน variational algorithm คือการใช้ Qiskit Runtime Session การเปิด Session ช่วยให้คุณรัน iteration หลาย ๆ ครั้งของ variational algorithm โดยไม่ต้องรอคิวใหม่ทุกครั้งที่พารามิเตอร์ถูกอัปเดต ซึ่งสำคัญมากถ้าเวลารอคิวยาวนานและ/หรือต้องใช้หลาย iteration เฉพาะพาร์ตเนอร์ใน IBM Quantum® Network เท่านั้นที่สามารถใช้ Runtime Session ได้ ถ้าคุณไม่มีสิทธิ์เข้าถึง Session คุณสามารถลดจำนวน iteration ที่ส่งในแต่ละครั้ง และบันทึกพารามิเตอร์ล่าสุดเพื่อใช้ในการรันครั้งต่อไป ถ้าส่ง iteration มากเกินไปหรือเจอเวลาคิวที่นานเกินไป คุณอาจพบ error code 1217 ซึ่งหมายถึงความล่าช้าระหว่างการส่ง job นาน

# Estimated usage: < 1 min. Benchmarked at 7 seconds on an Eagle processor
# Load necessary packages:

from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    Session,
    EstimatorOptions,
    EstimatorV2 as Estimator,
)
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

# Select the least busy backend:

service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
)
# Or get a specific backend:
# backend = service.backend("ibm_brisbane")

# Use a pass manager to transpile the circuit and observable for the specific backend being used:

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

# Set estimator options
estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

# Open a Runtime session:

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)
    cost = cost_func_vqe(theta_list, isa_ansatz, isa_observable, estimator)

session.close()
print(cost)

สังเกตว่าค่าที่ได้จากการคำนวณทั้งสองข้างต้นมีความคล้ายกันมาก เทคนิคสำหรับการปรับปรุงผลลัพธ์จะถูกพูดถึงเพิ่มเติมด้านล่าง

Example mapping to non-physical systems

ปัญหา maximum cut (Max-Cut) เป็นปัญหา combinatorial optimization ที่เกี่ยวข้องกับการแบ่ง vertex ของกราฟออกเป็นสองชุดที่ไม่ซ้อนทับกัน เพื่อให้จำนวน edge ระหว่างสองชุดมีค่าสูงสุด พูดอย่างเป็นทางการกว่านั้น กำหนดกราฟไม่มีทิศทาง $G=(V,E)$ โดยที่ $V$ คือชุดของ vertex และ $E$ คือชุดของ edge ปัญหา Max-Cut ขอให้แบ่ง vertex ออกเป็นสองชุดย่อยที่ไม่ซ้อนทับกัน $S$ และ $T$ เพื่อให้จำนวน edge ที่มีปลายด้านหนึ่งใน $S$ และอีกด้านใน $T$ มีค่าสูงสุด

เราสามารถนำ Max-Cut ไปใช้แก้ปัญหาต่าง ๆ ได้ เช่น: clustering, network design, phase transitions และอื่น ๆ เริ่มต้นด้วยการสร้างกราฟของปัญหา:

import rustworkx as rx
from rustworkx.visualization import mpl_draw

n = 4
G = rx.PyGraph()
G.add_nodes_from(range(n))
# The edge syntax is (start, end, weight)
edges = [(0, 1, 1.0), (0, 2, 1.0), (0, 3, 1.0), (1, 2, 1.0), (2, 3, 1.0)]
G.add_edges_from(edges)

mpl_draw(
    G, pos=rx.shell_layout(G), with_labels=True, edge_labels=str, node_color="#1192E8"
)

ปัญหานี้สามารถแสดงในรูปแบบ binary optimization ได้ สำหรับแต่ละ node $0 \leq i < n$ โดยที่ $n$ คือจำนวน node ของกราฟ (ในกรณีนี้ $n=4$) เราจะพิจารณาตัวแปร binary $x_i$ ตัวแปรนี้จะมีค่า $1$ ถ้า node $i$ อยู่ในกลุ่มที่เราจะเรียกว่า $1$ และ $0$ ถ้าอยู่ในกลุ่มอื่นที่เราจะเรียกว่า $0$ นอกจากนี้เราจะแทน $w_{ij}$ (ตำแหน่ง $(i,j)$ ของ adjacency matrix $w$) เป็นน้ำหนักของ edge ที่เชื่อมจาก node $i$ ถึง node $j$ เนื่องจากกราฟไม่มีทิศทาง $w_{ij}=w_{ji}$ จากนั้นเราสามารถกำหนดปัญหาเป็นการ maximize cost function ต่อไปนี้:

$$\begin{aligned} C(\vec{x}) & =\sum_{i,j=0}^n w_{ij} x_i(1-x_j)\\[1mm] & = \sum_{i,j=0}^n w_{ij} x_i - \sum_{i,j=0}^n w_{ij} x_ix_j\\[1mm] & = \sum_{i,j=0}^n w_{ij} x_i - \sum_{i=0}^n \sum_{j=0}^i 2w_{ij} x_ix_j \end{aligned}$$

เพื่อแก้ปัญหานี้ด้วยคอมพิวเตอร์ควอนตัม เราจะแสดง cost function ในรูปของ expected value ของ observable อย่างไรก็ตาม observable ที่ Qiskit รองรับโดยตรงประกอบด้วย Pauli operator ที่มี eigenvalue $1$ และ $-1$ แทนที่จะเป็น $0$ และ $1$ นั่นคือเหตุผลที่เราจะทำการเปลี่ยนตัวแปรดังต่อไปนี้:

โดยที่ $\vec{x}=(x_0,x_1,\cdots ,x_{n-1})$ เราสามารถใช้ adjacency matrix $w$ เพื่อเข้าถึงน้ำหนักของ edge ทั้งหมดได้อย่างสะดวก ซึ่งจะใช้ในการหา cost function ของเรา:

$$z_i = 1-2x_i \rightarrow x_i = \frac{1-z_i}{2}$$

ซึ่งหมายความว่า:

$$\begin{array}{lcl} x_i=0 & \rightarrow & z_i=1 \\ x_i=1 & \rightarrow & z_i=-1.\end{array}$$

ดังนั้น cost function ใหม่ที่เราต้องการ maximize คือ:

$$\begin{aligned} C(\vec{z}) & = \sum_{i,j=0}^n w_{ij} \bigg(\frac{1-z_i}{2}\bigg)\bigg(1-\frac{1-z_j}{2}\bigg)\\[1mm] & = \sum_{i,j=0}^n \frac{w_{ij}}{4} - \sum_{i,j=0}^n \frac{w_{ij}}{4} z_iz_j\\[1mm] & = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} z_iz_j \end{aligned}$$

ยิ่งกว่านั้น แนวโน้มตามธรรมชาติของคอมพิวเตอร์ควอนตัมคือการหาค่าต่ำสุด (โดยปกติคือพลังงานต่ำสุด) แทนที่จะหาค่าสูงสุด ดังนั้นแทนที่จะ maximize $C(\vec{z})$ เราจะ minimize:

$$-C(\vec{z}) = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} z_iz_j - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2}$$

เมื่อเรามี cost function ที่จะ minimize แล้ว โดยที่ตัวแปรมีค่าได้ $-1$ และ $1$ เราสามารถทำการเปรียบเทียบกับ Pauli $Z$ ดังนี้:

$$z_i \equiv Z_i = \overbrace{I}^{n-1}\otimes ... \otimes \overbrace{Z}^{i} \otimes ... \otimes \overbrace{I}^{0}$$

กล่าวอีกนัยหนึ่ง ตัวแปร $z_i$ จะเทียบเท่ากับ $Z$ Gate ที่กระทำบน Qubit $i$ ยิ่งไปกว่านั้น:

$$Z_i|x_{n-1}\cdots x_0\rangle = z_i|x_{n-1}\cdots x_0\rangle \rightarrow \langle x_{n-1}\cdots x_0 |Z_i|x_{n-1}\cdots x_0\rangle = z_i$$

จากนั้น observable ที่เราจะพิจารณาคือ:

$$\hat{H} = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} Z_iZ_j$$

ซึ่งเราจะต้องบวก independent term ทีหลัง:

$$\texttt{offset} = - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2}$$

operator คือ linear combination ของ term ที่มี Z operator บน node ที่เชื่อมต่อกันด้วย edge (จำว่า Qubit ที่ 0 อยู่ขวาสุด): $IIZZ + IZIZ + IZZI + ZIIZ + ZZII$ เมื่อสร้าง operator แล้ว ansatz สำหรับ algorithm QAOA สามารถสร้างได้ง่าย ๆ โดยใช้ QAOAAnsatz Circuit จาก Qiskit circuit library

from qiskit.circuit.library import QAOAAnsatz
from qiskit.quantum_info import SparsePauliOp

hamiltonian = SparsePauliOp.from_list(
    [("IIZZ", 1), ("IZIZ", 1), ("IZZI", 1), ("ZIIZ", 1), ("ZZII", 1)]
)

ansatz = QAOAAnsatz(hamiltonian, reps=2)
# Draw
ansatz.decompose(reps=3).draw("mpl")

# Sum the weights, and divide by 2

offset = -sum(edge[2] for edge in edges) / 2
print(f"""Offset: {offset}""")

Offset: -2.5

เนื่องจาก Runtime Estimator รับ Hamiltonian และ parameterized ansatz โดยตรงและคืนค่าพลังงานที่จำเป็น cost function สำหรับ QAOA instance จึงค่อนข้างง่าย:

def cost_func(params, 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
    """
    pub = (ansatz, hamiltonian, params)
    cost = estimator.run([pub]).result()[0].data.evs
    #    cost = estimator.run(ansatz, hamiltonian, parameter_values=params).result().values[0]
    return cost

import numpy as np

x0 = 2 * np.pi * np.random.rand(ansatz.num_parameters)

estimator = StatevectorEstimator()
cost = cost_func_vqe(x0, ansatz, hamiltonian, estimator)
print(cost)

1.473098768180865

# Estimated usage: < 1 min, benchmarked at 6 seconds on ibm_osaka, 5-23-24
# Load some necessary packages:

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import Session, EstimatorV2 as Estimator

# Select the least busy backend:

backend = service.least_busy(
    operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
)

# Or get a specific backend:
# backend = service.backend("ibm_brisbane")

# Use a pass manager to transpile the circuit and observable for the specific backend being used:

pm = generate_preset_pass_manager(backend=backend, optimization_level=1)
isa_ansatz = pm.run(ansatz)
isa_hamiltonian = hamiltonian.apply_layout(layout=isa_ansatz.layout)

# Set estimator options
estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

# Open a Runtime session:

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)
    cost = cost_func_vqe(x0, isa_ansatz, isa_hamiltonian, estimator)

# Close session after done
session.close()
print(cost)

1.1120776913677988

เราจะกลับมาดูตัวอย่างนี้อีกครั้งใน Applications เพื่อสำรวจวิธีใช้ optimizer ในการวนซ้ำผ่าน search space โดยทั่วไปแล้ว ซึ่งรวมถึง:

• การใช้ optimizer เพื่อหาพารามิเตอร์ที่เหมาะสมที่สุด
• การ bind พารามิเตอร์ที่เหมาะสมลงใน ansatz เพื่อหา eigenvalue
• การแปลง eigenvalue กลับสู่นิยามของปัญหา

กลยุทธ์การวัด: ความเร็ว vs ความแม่นยำ

ดังที่กล่าวไว้ก่อนหน้า เราใช้คอมพิวเตอร์ควอนตัมที่มีสัญญาณรบกวนเป็น oracle แบบกล่องดำ ซึ่งสัญญาณรบกวนอาจทำให้ค่าที่ดึงมาไม่แน่นอน ก่อให้เกิดความผันผวนแบบสุ่มที่จะส่งผลเสีย — หรืออาจขัดขวางโดยสิ้นเชิง — ต่อการลู่เข้าของ optimizer บางตัวไปยังคำตอบที่เสนอ นี่คือปัญหาทั่วไปที่เราต้องแก้ไขขณะที่ค่อยๆ สำรวจประโยชน์ใช้สอยของควอนตัมและก้าวสู่ความได้เปรียบเชิงควอนตัม:

เราสามารถใช้ตัวเลือกการระงับข้อผิดพลาดและการบรรเทาข้อผิดพลาดของ Qiskit Runtime Primitive เพื่อจัดการกับสัญญาณรบกวนและเพิ่มประโยชน์ใช้สอยของคอมพิวเตอร์ควอนตัมในปัจจุบันให้สูงสุด

การระงับข้อผิดพลาด

การระงับข้อผิดพลาด หมายถึงเทคนิคที่ใช้ปรับปรุงและแปลง Circuit ในระหว่างการ compile เพื่อลดข้อผิดพลาดให้น้อยที่สุด นี่คือเทคนิคการจัดการข้อผิดพลาดพื้นฐานที่มักส่งผลให้มี overhead ในการประมวลผลล่วงหน้าเชิงคลาสสิกต่อรันไทม์โดยรวม overhead รวมถึงการ transpile Circuit เพื่อรันบนฮาร์ดแวร์ควอนตัมโดย:

• แสดง Circuit โดยใช้ Gate พื้นฐานที่ระบบควอนตัมรองรับ
• แมป virtual qubit ไปยัง physical qubit
• เพิ่ม SWAP ตามข้อกำหนดด้านการเชื่อมต่อ
• ปรับ Gate แบบ 1Q และ 2Q
• เพิ่ม dynamical decoupling ให้กับ qubit ที่ไม่ได้ใช้งานเพื่อป้องกันผลของ decoherence

Primitive ช่วยให้สามารถใช้เทคนิคการระงับข้อผิดพลาดได้โดยการตั้งค่าตัวเลือก optimization_level และเลือกตัวเลือกการ transpile ขั้นสูง ในหลักสูตรถัดไป เราจะเจาะลึกวิธีการสร้าง Circuit ต่างๆ เพื่อปรับปรุงผลลัพธ์ แต่สำหรับกรณีส่วนใหญ่ เราแนะนำให้ตั้ง optimization_level=3

เราจะแสดงให้เห็นคุณค่าของการเพิ่มการปรับปรุงในกระบวนการ transpile โดยดูตัวอย่าง Circuit ที่มีพฤติกรรมอุดมคติแบบง่าย

from qiskit.circuit import Parameter, QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

theta = Parameter("theta")

qc = QuantumCircuit(2)
qc.x(1)
qc.h(0)
qc.cp(theta, 0, 1)
qc.h(0)
observables = SparsePauliOp.from_list([("ZZ", 1)])

qc.draw("mpl")

Circuit ข้างต้นสามารถให้ค่าคาดหวังแบบไซน์ซอยด์ของ observable ที่กำหนด โดยมีเงื่อนไขว่าเราแทรก phase ที่ครอบคลุมช่วงที่เหมาะสม เช่น $[0,2\pi]$

## Setup phases
import numpy as np

phases = np.linspace(0, 2 * np.pi, 50)

# phases need to be expressed as a list of lists in order to work
individual_phases = [[phase] for phase in phases]

เราสามารถใช้ simulator เพื่อแสดงให้เห็นประโยชน์ของการ transpile ที่ปรับปรุงแล้ว เราจะกลับมาใช้ฮาร์ดแวร์จริงด้านล่างเพื่อแสดงให้เห็นประโยชน์ของการบรรเทาข้อผิดพลาด เราจะใช้ QiskitRuntimeService เพื่อรับ Backend จริง (ในกรณีนี้คือ ibm_brisbane) และใช้ AerSimulator เพื่อจำลอง Backend นั้น รวมถึงพฤติกรรมสัญญาณรบกวนของมัน

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_aer import AerSimulator

# get a real backend from the runtime service
service = QiskitRuntimeService()
backend = service.backend("ibm_brisbane")

# generate a simulator that mimics the real quantum system with the latest calibration results
backend_sim = AerSimulator.from_backend(backend)

ตอนนี้เราสามารถใช้ pass manager เพื่อ transpile Circuit ไปยัง "instruction set architecture" หรือ ISA ของ Backend นี่คือข้อกำหนดใหม่ใน Qiskit Runtime: Circuit ทั้งหมดที่ส่งไปยัง Backend ต้องเป็นไปตามข้อจำกัดของ Backend target นั่นคือต้องเขียนในรูปแบบของ ISA ของ Backend — กล่าวคือ ชุดคำสั่งที่อุปกรณ์สามารถเข้าใจและดำเนินการได้ ข้อจำกัด target เหล่านี้ถูกกำหนดโดยปัจจัยต่างๆ เช่น basis gate พื้นฐานของอุปกรณ์ การเชื่อมต่อ qubit และเมื่อเกี่ยวข้อง ข้อกำหนดด้านเวลาของ pulse และคำสั่งอื่นๆ

สังเกตว่าในกรณีนี้ เราจะทำสิ่งนี้สองครั้ง: ครั้งแรกด้วย optimization_level = 0 และครั้งที่สองโดยตั้งค่าเป็น 3 แต่ละครั้งเราจะใช้ Estimator primitive เพื่อประมาณค่าคาดหวังของ observable ที่ค่า phase ต่างๆ

# Import estimator and specify that we are using the simulated backend:

from qiskit_ibm_runtime import EstimatorV2 as Estimator

estimator = Estimator(mode=backend_sim)

circuit = qc

# Use a pass manager to transpile the circuit and observable for the backend being simulated.
# Start with no optimization:

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

pm = generate_preset_pass_manager(backend=backend_sim, optimization_level=0)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

noisy_exp_values = []
pub = (isa_circuit, isa_observables, [individual_phases])
cost = estimator.run([pub]).result()[0].data.evs
noisy_exp_values = cost[0]

# Repeat above steps, but now with optimization = 3:

exp_values_with_opt_es = []
pm = generate_preset_pass_manager(backend=backend_sim, optimization_level=3)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

pub = (isa_circuit, isa_observables, [individual_phases])
cost = estimator.run([pub]).result()[0].data.evs
exp_values_with_opt_es = cost[0]

สุดท้าย เราสามารถพล็อตผลลัพธ์ และเห็นว่าความแม่นยำของการคำนวณค่อนข้างดีแม้ไม่มีการปรับปรุง แต่มันแน่นอนว่าดีขึ้นเมื่อเพิ่มการปรับปรุงเป็นระดับ 3 สังเกตว่าใน Circuit ที่ลึกและซับซ้อนกว่า ความแตกต่างระหว่าง optimization level 0 และ 3 มักจะชัดเจนกว่ามาก นี่คือ Circuit ง่ายมากที่ใช้เป็น toy model

import matplotlib.pyplot as plt

plt.plot(phases, noisy_exp_values, "o", label="opt=0")
plt.plot(phases, exp_values_with_opt_es, "o", label="opt=3")
plt.plot(phases, 2 * np.sin(phases / 2) ** 2 - 1, label="ideal")
plt.ylabel("Expectation")
plt.legend()
plt.show()

การบรรเทาข้อผิดพลาด

การบรรเทาข้อผิดพลาด หมายถึงเทคนิคที่ช่วยให้ผู้ใช้ลดข้อผิดพลาดของ Circuit โดยการสร้างแบบจำลองสัญญาณรบกวนของอุปกรณ์ในขณะที่รัน โดยทั่วไปแล้ว สิ่งนี้ส่งผลให้มี overhead ในการประมวลผลล่วงหน้าเชิงควอนตัมที่เกี่ยวกับการฝึกโมเดล และ overhead ในการประมวลผลหลังเชิงคลาสสิกเพื่อบรรเทาข้อผิดพลาดในผลลัพธ์ดิบโดยใช้โมเดลที่สร้างขึ้น

ตัวเลือก resilience_level ของ Qiskit Runtime primitive ระบุปริมาณความยืดหยุ่นที่จะสร้างขึ้นต่อข้อผิดพลาด ระดับที่สูงขึ้นจะสร้างผลลัพธ์ที่แม่นยำกว่าแต่แลกมาด้วยเวลาประมวลผลที่นานขึ้นเนื่องจาก overhead ในการสุ่มตัวอย่างเชิงควอนตัม ระดับความยืดหยุ่นสามารถใช้กำหนดค่าการแลกเปลี่ยนระหว่างต้นทุนและความแม่นยำเมื่อนำการบรรเทาข้อผิดพลาดไปใช้กับ primitive query ของคุณ

เมื่อใช้เทคนิคการบรรเทาข้อผิดพลาดใดๆ เราคาดว่า bias ในผลลัพธ์ของเราจะลดลงเมื่อเทียบกับ bias ก่อนหน้าที่ไม่ได้รับการบรรเทา ในบางกรณี bias อาจหายไปเลย อย่างไรก็ตาม สิ่งนี้มาพร้อมกับต้นทุน เมื่อเราลด bias ในปริมาณที่ประมาณ ความแปรปรวนทางสถิติจะเพิ่มขึ้น (นั่นคือ variance) ซึ่งเราสามารถจัดการได้โดยเพิ่มจำนวน shot ต่อ Circuit ในกระบวนการสุ่มตัวอย่างของเรา สิ่งนี้จะเพิ่ม overhead เกินกว่าที่จำเป็นสำหรับการลด bias จึงไม่ได้ทำโดยค่าเริ่มต้น เราสามารถเลือกใช้พฤติกรรมนี้ได้ง่ายๆ โดยปรับจำนวน shot ต่อ Circuit ใน options.executions.shots ดังแสดงในตัวอย่างด้านล่าง

สำหรับหลักสูตรนี้ เราจะสำรวจโมเดลการบรรเทาข้อผิดพลาดเหล่านี้ในระดับสูงเพื่อแสดงให้เห็นการบรรเทาข้อผิดพลาดที่ Qiskit Runtime primitive สามารถทำได้โดยไม่ต้องอธิบายรายละเอียดการใช้งานทั้งหมด

การดับข้อผิดพลาดการอ่านค่าแบบ Twirled (T-REx)

Twirled readout error extinction (T-REx) ใช้เทคนิคที่เรียกว่า Pauli twirling เพื่อลดสัญญาณรบกวนที่เกิดขึ้นระหว่างกระบวนการวัดค่าเชิงควอนตัม เทคนิคนี้ไม่ได้สมมติรูปแบบเฉพาะของสัญญาณรบกวน ทำให้มีความทั่วไปและมีประสิทธิภาพสูง

ขั้นตอนการทำงานโดยรวม:

1. รวบรวมข้อมูลสำหรับ zero state ด้วย bit flip แบบสุ่ม (Pauli X ก่อนการวัด)
2. รวบรวมข้อมูลสำหรับ state ที่ต้องการ (มีสัญญาณรบกวน) ด้วย bit flip แบบสุ่ม (Pauli X ก่อนการวัด)
3. คำนวณฟังก์ชันพิเศษสำหรับชุดข้อมูลแต่ละชุด แล้วหาร

 

เราสามารถตั้งค่านี้ด้วย options.resilience_level = 1 ดังแสดงในตัวอย่างด้านล่าง

การประมาณค่าแบบ Zero noise

Zero noise extrapolation (ZNE) ทำงานโดยการขยายสัญญาณรบกวนใน Circuit ที่เตรียม quantum state ที่ต้องการก่อน รับการวัดที่ระดับสัญญาณรบกวนต่างๆ หลายระดับ และใช้การวัดเหล่านั้นเพื่ออนุมานผลลัพธ์ที่ไม่มีสัญญาณรบกวน

ขั้นตอนการทำงานโดยรวม:

1. ขยายสัญญาณรบกวนของ Circuit สำหรับ noise factor หลายค่า
2. รัน Circuit ที่มีการขยายสัญญาณรบกวนทุกตัว
3. ประมาณค่าย้อนกลับไปยังขีดจำกัด zero noise

 

เราสามารถตั้งค่านี้ด้วย options.resilience_level = 2 เราสามารถปรับปรุงเพิ่มเติมได้โดยสำรวจ noise_factors, noise_amplifiers, และ extrapolators หลากหลาย แต่สิ่งนี้อยู่นอกขอบเขตของหลักสูตรนี้ เราขอแนะนำให้ทดลองกับตัวเลือกเหล่านี้ตามที่อธิบายไว้ที่นี่

แต่ละวิธีมาพร้อมกับ overhead ที่เกี่ยวข้อง: การแลกเปลี่ยนระหว่างจำนวนการคำนวณควอนตัมที่จำเป็น (เวลา) และความแม่นยำของผลลัพธ์:

$$\begin{array}{c|c|c|c} \text{Methods} & R=1 \text{, T-REx} & R=2 \text{, ZNE} \\[1mm] \hline \text{Assumptions} & \text{None} & \text{Ability to scale noise} \\[1mm] \text{Qubit overhead} & 1 & 1 \\[1mm] \text{Sampling overhead} & 2 & N_{\text{noise-factors}} \\[1mm] \text{Bias} & 0 & \mathcal{O}(\lambda^{N_{\text{noise-factors}}}) \\[1mm] \end{array}$$

การใช้ตัวเลือกการบรรเทาและการระงับของ Qiskit Runtime

นี่คือวิธีคำนวณค่าคาดหวังขณะใช้การบรรเทาข้อผิดพลาดและการระงับข้อผิดพลาดใน Qiskit Runtime เราสามารถใช้ Circuit และ observable เดียวกันกับก่อนหน้าได้ แต่คราวนี้คงระดับการปรับปรุงไว้ที่ระดับ 2 และปรับ resilience หรือเทคนิคการบรรเทาข้อผิดพลาดที่ใช้อยู่ กระบวนการบรรเทาข้อผิดพลาดนี้เกิดขึ้นหลายครั้งตลอด optimization loop

เราดำเนินการส่วนนี้บนฮาร์ดแวร์จริง เนื่องจากการบรรเทาข้อผิดพลาดไม่สามารถใช้ได้กับ simulator

# Estimated usage: 8 minutes, benchmarked on an Eagle processor, 5-23-24

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import (
    Session,
    EstimatorOptions,
    EstimatorV2 as Estimator,
)

# We select the least busy backend

# Select the least busy backend
# backend = service.least_busy(
#    operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
# )

# Or use a specific backend
backend = service.backend("ibm_brisbane")

# Initialize some variables to save the results from different runs:

exp_values_with_em0_es = []
exp_values_with_em1_es = []
exp_values_with_em2_es = []

# Use a pass manager to optimize the circuit and observables for the backend chosen:

pm = generate_preset_pass_manager(backend=backend, optimization_level=2)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

# Open a session and run with no error mitigation:

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

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

    pub = (isa_circuit, isa_observables, [individual_phases])
    cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em0_es = cost[0]

# Open a session and run with resilience = 1:

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

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

    pub = (isa_circuit, isa_observables, [individual_phases])
    cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em1_es = cost[0]

# Open a session and run with resilience = 2:

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

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

    pub = (isa_circuit, isa_observables, [individual_phases])
    cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em2_es = cost[0]

เช่นเดียวกับก่อนหน้า เราสามารถพล็อตค่าคาดหวังที่ได้เป็นฟังก์ชันของมุม phase สำหรับระดับการบรรเทาข้อผิดพลาดสามระดับที่ใช้ จะเห็นได้ยากนิดหน่อยว่าการบรรเทาข้อผิดพลาดช่วยปรับปรุงผลลัพธ์เล็กน้อย อีกครั้ง ผลกระทบนี้จะเด่นชัดกว่ามากใน Circuit ที่ลึกและซับซ้อนกว่า

import matplotlib.pyplot as plt

plt.plot(phases, exp_values_with_em0_es, "o", label="unmitigated")
plt.plot(phases, exp_values_with_em1_es, "o", label="resil = 1")
plt.plot(phases, exp_values_with_em2_es, "o", label="resil = 2")
plt.plot(phases, 2 * np.sin(phases / 2) ** 2 - 1, label="ideal")
plt.ylabel("Expectation")
plt.legend()
plt.show()

สรุป

ในบทเรียนนี้ คุณได้เรียนรู้วิธีสร้าง cost function:

• สร้าง cost function
• วิธีใช้ Qiskit Runtime primitive เพื่อบรรเทาและระงับสัญญาณรบกวน
• วิธีกำหนดกลยุทธ์การวัดเพื่อปรับความเร็วและความแม่นยำให้เหมาะสม

นี่คือ variational workload ระดับสูงของเรา:

cost function ของเราทำงานในทุก iteration ของ optimization loop บทเรียนถัดไปจะสำรวจว่า classical optimizer ใช้การประเมิน cost function ของเราอย่างไรเพื่อเลือกพารามิเตอร์ใหม่

import qiskit
import qiskit_ibm_runtime

print(qiskit.version.get_version_info())
print(qiskit_ibm_runtime.version.get_version_info())

1.1.0
0.23.0