IBM 與洛克希德·馬丁合作在量子電腦上模擬甲烯

# IBM 與洛克希德·馬丁合作在量子電腦上模擬甲烯今年五月份，IBM 與美國軍火商 Lockheed Martin 合作在量子電腦上使用 52 qubits 執行 SQD 演算法（sample-based quantum diagonalization）預測甲烯（methylene）能量，是首次將 SQD 演算法應用在 open shell 系統上，研究發表於 [《Journal of Chemical Theory and Computation》](https://pubs.acs.org/doi/10.1021/acs.jctc.5c00075)。

Open shell 就是指有未成對電子

## 前情提要在量子化學中，我們最常在鑽研的方程式就是以下知名的 Schrodinger 方程式： \begin{align}\tag{1} H|\psi\rangle=E|\psi\rangle \end{align} $H$ 會是一個矩陣，包含分子的所有能量資訊，$|\psi\rangle$ 也可以想成一個矩陣或向量，代表該分子的軌域（orbitals），兩者相乘後得到右邊一個數字 $E$ 乘上原來的向量 $|\psi\rangle$，其中 $E$ 就是該分子在這 orbitals 組合 $|\psi\rangle$ 下的能量值。然而當分子越來越大，結構越來越複雜，$H$ 矩陣會越來越大，計算（1）式所需的電腦資源會劇增，以知名的 FCI 方法（full configuration interaction）為例，想要在 22 個電子系統上獲得精確解已經逼近極限，想要再應用到更大的分子就需要用到許多近似方法，這也是長久以來計算化學在研究的工作之一。

由於這是新聞導讀，在很多專有名詞上為了科普會做很多簡化，以剛剛提到的「精確解」其實是指 exact diagonalization，22 個電子其實更明確說是 22 個電子搭配 22 個 orbitals 或是 26 個電子搭配 23 個 orbitals

而這時候，量子電腦給出了嶄新的一條路徑來攻克這困難，然而現在的量子電腦都不完美，造成以下兩個問題： 1. 為了得到更準確的結果就必須要對同的電路測量非常多次，測量次數越高，整體計算時間就拉長 2. 針對化學問題設計的 [ansatz](https://www.entangletech.tw/lesson/optim-08#toc-8)，其[電路深度](https://www.entangletech.tw/lesson/hardware-general-01)都很深，以知名的 coupled cluster 為例，對於 $M$ 個 (spin) orbitals 的分子，電路深度就達到 $M^4$，這麼深的電路意味著有非常多 gate，gate 越多，錯誤率也會隨之提高為了解決這問題，IBM 團隊提出 Sample-Based Quantum Diagonalization （SQD）演算法試圖克服量子電腦應用在化學上的瓶頸，在解釋 SQD 之前，榮我們來介紹 [VQE](https://www.entangletech.tw/lesson/optim-07) 演算法。 ## VQE 演算法在我們的[量子優化系列](https://www.entangletech.tw/lesson/optim-07)文章就有深入淺出介紹 VQE 演算法，有興趣的讀者可以點開連結，這麼我們就做簡短介紹。假設今天一個分子的 $H$ 可以寫成以下矩陣： \begin{align} H= \begin{bmatrix} -1.2 & 0 & 0 & 0.3 \\ 0 & -0.5 & 0.2 & 0 \\ 0 & 0.2 & -0.5 & 0 \\ 0.3 & 0 & 0 & 0.1 \\ \end{bmatrix} \end{align} VQE 的電路負責製備量子態 $|\psi\rangle$，接著經典電腦利用這量子態和 $H$ 矩陣計算能量值： \begin{align} \langle\psi|H|\psi\rangle=E \end{align} 透過此能量值和量子態提供新的參數傳回量子電路，如此循環直到 $E$ 沒什麼變化（簡稱收斂）。比方說我們這次的 VQE 電路長這樣：在某一次循環中，經典電腦給出的參數為$（-2\pi, 2\pi)$，那當下電路輸出的量子態為： \begin{align} RX(-2\pi)_0 \otimes RY(2\pi)_1|00\rangle= \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \\ 0 \\ 0 \\ \end{bmatrix}= \begin{bmatrix} 1 \\ 0 \\ 0 \\ 0 \\ \end{bmatrix} \end{align} 利用這結果計算能量值為： \begin{align} \langle \psi|H|\psi\rangle= \begin{bmatrix} 1 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} -1.2 & 0 & 0 & 0.3 \\ 0 & -0.5 & 0.2 & 0 \\ 0 & 0.2 & -0.5 & 0 \\ 0.3 & 0 & 0 & 0.1 \\ \end{bmatrix} \begin{bmatrix} 1 \\ 0 \\ 0 \\ 0 \\ \end{bmatrix}=-0.7 \end{align} 得到 $E=-0.7$，然後經典電腦再根據此提供新參數去調整電路的角度，理應最後會得到 $-1.26$。然而 VQE 的其中一個挑戰是，當 $H$ 超級大時，直接計算上式非常耗費電腦資源，而 SQD 演算法就是為此而生。 ## SQD 演算法 VQE 電路生成的量子態經過多次測量後便會進入下圖所繪的流程圖，而圖中繪製的流程圖都在經典電腦上運行：

SQD 流程圖

Picture comes from doi: 10.1126/sciadv.adu999

首先，因為現在的量子電腦都不完美，會有 noise 影響到製備的量子態，因此第一步是要檢查測量得到的數據有沒有受到 noise 干擾，找到被干擾的部分並做糾正，第二步從這大量的測量數據中隨機抽幾個樣本形成一個 batch $S$，共會抽 $K$ 個 batch，利用每個 batch 的數據對 $H$ 做投影和[對角化](https://www.entangletech.tw/lesson/math-08#toc-4)，得到比較小的 $H_S$，再利用這 $H_S$ 解出 $E$，$K$ 個 batch 就會得到 $K$ 個 $E$，從中挑出最低能量值作為這一次循環算出的能量值，然後再重複第一步驟直到 $E$ 收斂為止，按照經驗，通常三次循環就收斂了。 **有沒有發現最神奇的地方**，與 VQE 演算法不同，要再提供新的參數回傳給量子電路，如此循環，SQD 演算法僅需在**量子電路上執行一次**，接著全都在經典電腦上運行即可，不過在論文中也有提及，要達到這神奇的效果也要靠 ansatz 的設計，在這篇論文中團隊使用 LUCJ ansatz，ansatz 的細節不會在本文章提及，有興趣的讀者可以閱讀原論文，如果使用不好的 ansatz，那最後也得將 SQD 的結果給 optimizer 計算梯度，提供新參數給 VQE，如此循環。接下來我們會詳述 SQD 的每個步驟究竟怎麼運作： ### 第一步：recovery configuration 第一步就是要找到被 noise 干擾到的測量結果，並對之作修正：

多次測量得到測量結果，從中找到有錯的地方並做糾正

Picture comes from 2025 Qiskit Summer School

以下圖最左邊為例，在 IBM 的編碼規則中，$|0011\space 0011\rangle$ 的左邊四碼（$0011$）代表 spin down 的電子（藍色箭頭），右邊四碼代表 spin up 的電子（紅色箭頭），$0$ 代表該 orbital 上沒有電子，$1$ 則代表有電子，在閱讀 $0011$ 這編碼時，最右邊代表最低能量的 orbital，最左邊代表高能量的 orbital，因此 $0011$ 代表最低的兩個 orbitals 被兩個電子佔據。

Picture comes from 2025 Qiskit Summer School

然而有時候會因為 noise 的問題有個應該是 $0$ 的變成 $1$，像上圖的中間，使得 spin up 的電子原本應為 2 個變成 3 個，就必須從這 3 個 $1$ 中找出有錯的翻轉成 $0$，上圖最右也是一樣的道理，有個應該是 $1$ 變成 $0$，就要從 3 個 $0$ 中找到有錯的，然後翻轉成 $1$。而要知道究竟是哪個 qubit 要做糾正就得靠大量測量結果，累積統計每個 qubit 是 $1$ 的機率有多高，如果該 qubit 是 $1$ 的機率很大但當次測量得到是 $0$ 就對該 qubit 做糾正，在 SQD 中，第一次迭代並沒有這統計數據（這數據稱之為平均佔據數），所以在第一步會直接跳過此步驟。 ### 第二步：subsample 假設一個 VQE 電路測量了 100,000 次，做完糾正後接著從這十萬筆數據中隨機選 50 次測量結果組成一個批次（batch），共組合出 5 個 batch。

Picture comes from 2025 Qiskit Summer School

那為了方便演示 SQD 的運作，我們以前面提及的矩陣做範例，從大量測量結果中隨機抽 2 次測量結果組成一個 batch，分別是： \begin{align} |x_0\rangle = |01\rangle \\ |x_1\rangle = |10\rangle \end{align} ### 第三步：矩陣投影與對角化剛剛抽出的每個 batch 裡面的測量結果組合成投影矩陣 $P_S$，接著利用這投影矩陣將矩陣 $H$ 投影到 $H_S$，透過這樣的操作讓原本很龐大的矩陣縮小到很小的矩陣，對小矩陣做對角化得到能量值相比原本要對角化一個很龐大的矩陣來得容易。

以氮氣分子維利，其 H 矩陣為 40 億乘 40 億，對這樣的矩陣做對角化，十分不容易，但透過 SQD 的操作可以把這龐大的矩陣簡化到僅僅 50 乘 50

Picture comes from 2025 Qiskit Summer School

以我們剛剛的舉例，投影矩陣為： \begin{align} P_S&=|x_0\rangle\langle x_0 |+|x_1\rangle\langle x_1| \\&=|01\rangle\langle01|+|10\rangle\langle 10|\\ &= \begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ \end{bmatrix}+ \begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0\\ \end{bmatrix}\\ &= \begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0\\ \end{bmatrix} \end{align} 然後利用這投影矩陣去投影我們欲解的 $H$： \begin{align} H_S&=P_SHP_S\\ &=\begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0\\ \end{bmatrix} \begin{bmatrix} -1.2 & 0 & 0 & 0.3 \\ 0 & -0.5 & 0.2 & 0 \\ 0 & 0.2 & -0.5 & 0 \\ 0.3 & 0 & 0 & 0.1 \\ \end{bmatrix} \begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0\\ \end{bmatrix}= \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & -0.5 & 0.2 & 0 \\ 0 & 0.2 & -0.5 & 0 \\ 0 & 0 & 0 & 0 \\ \end{bmatrix}\\ &=\begin{bmatrix} -0.5 & 0.2 \\ 0.2 & -0.5 \\ \end{bmatrix} \end{align} 透過這樣的操作，原本 $4\times 4$ 的矩陣變成 $2\times 2$ 的矩陣。 ### 第四步：收集最低能量值，更新平均佔據數如果你閱讀過我們撰寫有關 [eigenvalue](https://www.entangletech.tw/lesson/math-08) 的文章，這麼簡單的矩陣你一定能解出他的 eigenvalue $E$ 為 $-0.3$ 與 $-0.7$，從中找到最低的能量 $-0.7$，然後統計這次選出的測量結果，得到每個 qubit 為 $1$ 的機率是多少，重複第一步驟，直到 $E$ 沒什麼變化。 ## 實驗結果在這篇論文之前，SQD 已經被廣泛應用在許多分子上，卻沒有用在 open shell 系統上過，因此這次 IBM 將其應用在計算 methylene 兩種構型之間的能量差異：

methylene 的兩種構型：singlet 與 triplet，兩種構型的能量差異為 14.4 mHa(9.03 kcal/mol)

在量子化學中，能量單位為 Ha（Hartree），也可以記做 Eh

當 methylene 的兩個孤對電子平行排列時，這狀態稱作 triplet，會被稱作 triplet 是因為當兩個 spin 平行排列時會有三個自旋量子數（這在圖中沒畫出來）；當兩個孤對電子反平行排列，這狀態稱作 singlet，此時這對孤對電子可能位在圖中灰色的 orbital 或是白色的 orbital（再一次，會被稱作 singlet 是因為只有一個自旋量子數，跟圖中有兩種不同的分布位置是兩回事）。這次的實驗在 *ibm nazca* 量子電腦上執行，模擬 6 個電子與 23 個 orbitals，因此需要 46 個 qubits 表示每個 orbitals 有沒有電子，以及其他 6 個 qubits（總共 52 個 qubits）和 33,242 個 two qubits gates，電路深度達到 1663，測量 100,000 次。 ### 能量預測與 SQD 一起做比較的為經典計算方法，包括 SCI、CCSD 與 SCF。

SCI：selected configuration interaction
CCSD：coupled cluster singles and doubles
SCF：self-consistent field

從下圖可以看到，對於 triplet 態，SQD 的結果與 SCI 非常接近（圖 (a)），與 SCI 的差異落在 $1\sim28$ mHa 範圍內（圖 (b)），當鍵長來到 2 Å 和 2.5 Å 時，兩者差異最大；對於 singlet，SQD 的結果與 SCI 都插在 $1\sim4$ mHa。

(a) 不同鍵長，不同方法計算兩態的能量值（b)不同鍵長，不同方法計算兩態的能量值與 SCI 差異（c）不同鍵長，不同方法計算兩態的能量差值

圖（c）展示每個鍵長下 singlet-triplet 兩態能量差異，在接近平衡鍵長（1.09 Å），三者預測的差值非常接近，在這裡，SQD 預測兩態能量差值 19 mHa，與 SCI 計算的 18 mHa 非常接近，實驗值為 14 mHa，而其他兩種經典方法 CCSD 和 SCF 則分別預測 24 和 40 mHa，體現出 SQD 的預測值非常接近 SCI 預測和實驗值。

不同方法預測在平衡鍵長下兩態的能量差異與實驗值

### 波函數振幅分析為了了解為何 SQD 在 triplet state 鍵長來到 2~2.5 Å 時就開始表現不佳，作者對 SQD 預測的量子態做分析，將每個量子態的振幅平方 $|\psi_x|^2$ 並做圖，並與 SCI 的預測做比較。

SCI（黃線）與 SQD（紅點）的波函數振幅分析。（a）為 singlet （b）為 triplet。從左到右分別是鍵長 1.11 Å, 1.60 Å, 2.30 Å, 與 2.50 Å

可以明顯發現，圖中越靠左邊，SQD 與 SCI 的預測越接近，越靠右邊，兩者預測越差越遠，另外當鍵長越長，兩者也越差越多，會造成這樣的現象，作者推測有兩種可能原因： 1. 當鍵長越大，相比平衡鍵長，波函數振幅越來越不集中 2. 越接近平衡鍵長，波函數展現出很強的平均場性質，而這樣的性質在鍵長越來越大時就隨之消失而這兩個效應影響到 SQD 計算出的平均佔據數以及如何做 configuration recovery，也就是說，當波函數被平均場主宰時，每個 qubit 的平均佔據數會非常接近 0 或 1，SQD 可以很好判斷每個 qubit 理應是 0 還是 1，導致 SQD 在平衡鍵長上的表現最好。

此文章僅作初步導讀，有更多內容在本文中沒有詳細提到，歡迎參看原論文

## 參考文獻與延伸閱讀 - [原論文](https://pubs.acs.org/doi/10.1021/acs.jctc.5c00075) - [IBM SQD 簡易教學](https://github.com/qiskit-community/qgss-2025/blob/main/solutions/lab-3/lab3-solution.ipynb) - [SQD 原始論文](https://www.science.org/doi/10.1126/sciadv.adu9991) - [SQD 原始碼](https://github.com/Qiskit/qiskit-addon-sqd/tree/main/qiskit_addon_sqd)