文献情報
- タイトル：Hardware-efficient quantum annealing with error mitigation via classical shadow
- 著者：Takaharu Yoshida, Yuta Shingu, Chihaya Shimada, Tetsuro Nikuni, Hideaki Hakoshima, Yuichiro Matsuzaki
- 書誌情報(DOI)： https://doi.org/10.1103/stb4-wms8
- 出版年：2026年6月10日
概要
量子アニーリングは、問題ハミルトニアンの基底状態を求める手法です。しかし、実際のアニーリングマシンでは、環境との相互作用によるノイズや、有限時間の時間発展に伴う非断熱遷移が生じるため、必ずしも基底状態が得られるとは限りません。その結果、問題ハミルトニアンのエネルギー期待値が基底エネルギーからずれることがあります。本記事では、このずれをバイアスと呼びます。
このバイアスを小さくする方法の1つが、量子エラー抑制法の一種である仮想蒸留法です。仮想蒸留法では、ノイズを含む状態の期待値の代わりに、その状態の密度行列を累乗して規格化した状態(純粋化状態)に対する期待値を求めます。ただし、純粋化状態の期待値計算に必要な量を求める従来手法（アダマールテスト）は、複数の量子ビットにまたがるゲート操作を必要とするため、ゲート型量子コンピュータ向けの手法となっており、ゲート操作ができないアニーリングマシンでは実装が困難でした。
そこで、この論文では、アダマールテストの代わりに古典シャドウ法を用いることを提案しました。古典シャドウ法は、各量子ビットへのランダムな $X, Y, Z$ 測定と古典的な後処理のみで期待値を推定する手法です。この手法を用いることで、複数量子ビットにまたがる制御ゲートを必要とせず、$Z$測定に加え、$X, Y$ 測定が可能なアニーリングマシンでは、仮想蒸留法を実行できるようになります。
本記事では、この論文の内容を解説するとともに、記事の後半では原論文で行われた数値実験の結果を紹介し、仮想蒸留法の振る舞いを確認します。
背景
量子アニーリングの原理
量子アニーリングは、次の時間依存ハミルトニアンを持つ系を時間発展させることで、問題ハミルトニアン $H_p$ の基底状態を得る方法です[1]。
$$
\tag{1}
H(t) = A(t) H_p + B(t) H_d, \quad A(t) := \frac{t}{T}, \quad B(t) := 1 – \frac{t}{T}
$$
ここで、$H_p$ は基底状態を求めたい問題ハミルトニアン（例: QUBOに対応するイジングハミルトニアン）、$H_d$ は量子揺らぎを与えるドライバーハミルトニアン（例: 横磁場）です。$T$ はアニーリング時間を表し、量子アニーリングでは、ドライバーハミルトニアンの基底状態を初期状態として準備し、アニーリング時間まで状態を時間発展させます。
このような時間発展によって基底状態が得られることは、断熱定理によって説明されます。
断熱定理によれば、系が環境と相互作用しない孤立系であり、基底状態と励起状態の間に適切なエネルギーギャップが存在するなどの条件を満たす場合、十分ゆっくりとハミルトニアンを変化させることで、量子状態は瞬間的な基底状態に従って変化します。
したがって、理想的な断熱極限 \(T\to\infty\) では、初期状態として \(H_d\) の基底状態を準備すると、最終的に \(H_p\) の基底状態が得られます。
このとき、アニーリング後の状態を \(\rho\)、\(H_p\) の基底エネルギーを \(E_g\) とすると、$
\operatorname{Tr}(H_p\rho)=E_g$が成り立ちます。
アニーリングマシンでのノイズ
しかし、実際のアニーリングマシン（例: D-Waveが提供しているアニーリングマシン）は、環境との相互作用を完全に避けることは難しく、時間発展も有限時間で行われます。そのため、理想的な断熱時間発展とは異なる状態が得られる場合があります。
ここでは、量子アニーリングにおける誤差の原因として、インコヒーレントノイズと非断熱遷移の2つを考えます。
1つ目は、環境との相互作用などによって生じるインコヒーレントノイズです。このようなノイズによって、量子状態は一般に純粋状態から混合状態へと変化します。混合状態とは、複数の純粋状態が確率的に混ざった状態として表される量子状態です。
2つ目は、有限時間の時間発展に伴う非断熱遷移です。ハミルトニアンの変化が十分に遅くない場合、量子状態は瞬間的な基底状態に完全には追従できず、励起状態への遷移が生じることがあります。環境との相互作用を考えない孤立系では、時間発展はユニタリ演算子によって記述されます。そのため、非断熱遷移が生じても、初期状態が純粋状態であれば、最終状態も純粋状態のままです。
インコヒーレントノイズによってアニーリング後の状態は一般に混合状態となるため、以降では量子状態を密度行列 \(\rho\) で表すことにします。密度行列は純粋状態と混合状態の両方を表すことができ、純粋状態 \(|\psi\rangle\) に対しては \(\rho=|\psi\rangle\langle\psi|\) と書けます。
また、量子状態 \(\rho\) に対する物理量 \(O\) の期待値は、$
\langle O\rangle=\operatorname{Tr}(O\rho)$で表されます。ここで$\text{Tr}$はトレースと呼ばれる演算で、演算子の固有値の和を意味します。
問題設定
以降、問題ハミルトニアンは、QUBOに対応するハミルトニアン(Z演算子のみで構成されるイジングハミルトニアン)だとします。
アニーリング後の状態 $\rho$ に対して、問題ハミルトニアン $H_p$ の期待値$\text{Tr}(H_p \rho)$を求めることを考えます。ただし、測定回数は有限なので、実際に得られるのは期待値そのものではなく、$M$ 回の測定から求めた平均値です。
$$
\tag{2}
E^{(M)} = \frac{1}{M} \sum_{i=1}^{M}
\left(
\rho\text{ に対する }i\text{ 回目の }H_p\text{ 測定値}
\right)
$$
一方、測定回数$M$を増やすと、大数の法則により、平均値E^{(M)}は期待値$\mathrm{Tr}[H_p \rho]$に収束します。
$$
\tag{3}
\mathrm{Tr}[H_p \rho] = \lim_{M \to \infty} E^{(M)}
$$
これらを踏まえた上で、\(M\) 回の測定から得られる平均値 \(E^{(M)}\) がどのような分布に従うか考えます。
図1：$M$回測定で得られる の平均値の分布
理想的には \(H_p\) の基底状態が得られて平均値は基底エネルギー \(E_g\) に一致するはずですが、インコヒーレントノイズと非断熱遷移のため、実機では平均値の中心 \(\operatorname{Tr}[H_p\rho]\) が基底エネルギー \(E_g\) からずれる場合があります。このずれがバイアスです。この論文では、量子エラー抑制の一種である仮想蒸留法を用いて、特にノイズに由来するバイアスを小さくすることを目指します。
量子エラー抑制とは、ノイズの影響を受けた測定結果に古典的な後処理などを施すことで、ノイズがない場合の期待値をより正確に推定する方法の総称です。量子エラー抑制法には、外挿法や擬確率法など様々な方法が提案されています。(各量子エラー抑制法については、第69回物性若手夏の学校のテキスト[2]が詳しいです。)
続いて、典型的な量子エラー抑制を適用した場合の期待値の分布を見てみましょう。
図2：量子エラー抑制を適用した場合の平均値の分布
エラー抑制を行うことで、バイアスを小さくできる場合がありますが、その代償として推定値の分散が大きくなることがあります。分散の増加に対しては、測定回数（ショット数）を増やして統計誤差を小さくすることで対処します。
手法
仮想蒸留法（FVP）
仮想蒸留法（Full Virtual Purification, FVP）は、$\rho$ に対する期待値の代わりに、純粋化状態 $\rho^n / \mathrm{Tr}[\rho^n]$ に対する期待値を求めることでエラー抑制を行う手法です[3, 4]。ここで $n$ はコピー数と呼ばれる正の整数です。
なぜこれでエラー抑制ができるのかを見てみます。$\rho$ が正規直交する固有ベクトルの集合 $\{|\psi_k\rangle\}_{k=1}^{2^N}$ を用いて固有値展開できるとします（$N$ は量子ビット数）。
$$
\tag{4}
\rho = \lambda_1 |\psi_1\rangle\langle\psi_1| + \sum_{k=2}^{2^N} \lambda_k |\psi_k\rangle\langle\psi_k|, \quad \lambda_1 > \lambda_2 \geq \lambda_3 \geq \cdots
$$
このとき $\rho^n$ を計算すると、次のようになります。
$$
\tag{5}
\rho^n = \lambda_1^n \left( |\psi_1\rangle\langle\psi_1| + \sum_{k=2}^{2^N} \left( \frac{\lambda_k}{\lambda_1} \right)^n |\psi_k\rangle\langle\psi_k| \right)
$$
$\lambda_1$ が最大固有値なので $\lambda_k / \lambda_1 < 1$ であり、$n$ を大きくするほど第2項以降の係数が小さくなり、$\rho^n$ を規格化した$\rho^n / \mathrm{Tr}[\rho^n]$ は最大固有ベクトル$|\psi_1\rangle$に対応する純粋状態 $|\psi_1\rangle\langle\psi_1|$ に近づきます。もし最大固有ベクトルが理想的な基底状態とほぼ一致していれば、すなわち $|\langle\psi_1 | E_g \rangle| \approx 1$ であれば、純粋化状態$\rho^n / \mathrm{Tr}[\rho^n]$ に対するエネルギー期待値は 元の状態$\rho$ に対するエネルギー期待値よりも基底エネルギーに近づくことが期待できます。
ここで注意したいのは、仮想蒸留法が混合状態を純粋状態に近づけることでエラー抑制を行う点です。したがって、仮想蒸留法はインコヒーレントノイズによる期待値のずれを抑制できる可能性がありますが、純粋状態のまま生じる非断熱遷移によるずれは抑制できません。
- インコヒーレントノイズ（純粋状態 → 混合状態）：仮想蒸留法で抑制できる場合がある
- 非断熱遷移（孤立系では純粋状態のまま）：仮想蒸留法で抑制不可
孤立系での非断熱遷移はユニタリな時間発展によって生じ、純粋状態を混合状態にはしないため、仮想蒸留法では抑制できず、期待値に誤差として残ります。
純粋化状態の期待値をどう計算するか
以降では、コピー数 $n=2$ の場合を考えます。求めたいのは次の期待値です。
$$
\tag{6}
\mathrm{Tr}\left[ H_p \frac{\rho^2}{\mathrm{Tr}[\rho^2]} \right] = \frac{\mathrm{Tr}[H_p \rho^2]}{\mathrm{Tr}[\rho^2]}
$$
純粋化状態 $\rho^2 / \mathrm{Tr}[\rho^2]$ を直接準備するのは困難です。そこで、分母 $\mathrm{Tr}[\rho^2]$ と分子 $\mathrm{Tr}[H_p \rho^2]$ を別々に推定する方針をとります。ここで、2つの量子状態を入れ替えるSWAP演算子を用いた次の恒等式を使います。
$$
\tag{7}
\mathrm{Tr}[O \rho^2] = \mathrm{Tr}\left[ (O \otimes I)U_{\mathrm{SWAP}}(\rho \otimes \rho) \right]
$$
ここで $U_{\mathrm{SWAP}}$ は2つの量子状態を入れ替える演算子で、$U_{\mathrm{SWAP}} (|\phi_1\rangle \otimes |\phi_2\rangle) := |\phi_2\rangle \otimes |\phi_1\rangle$ を満たします。分母は $O = I$、分子は $O = H_p$ とすればよく、結局 $\mathrm{Tr}[(O \otimes I) U_{\mathrm{SWAP}} (\rho \otimes \rho)]$ をどう推定するかが問題になります。
(7)式の証明
まず、 $U_{\mathrm{SWAP}}$ の定義を確認します。$U_{\mathrm{SWAP}}$ は2つの量子状態を入れ替える演算子で、正規直交基底 ${|i\rangle}$ を用いると
$$U_{\mathrm{SWAP}} = \sum_{i,j} |i\rangle\langle j| \otimes |j\rangle\langle i|$$
と書けます。実際、これを $|a\rangle \otimes |b\rangle$ に作用させると
$$U_{\mathrm{SWAP}} (|a\rangle \otimes |b\rangle) = \sum_{i,j} \langle j|a\rangle \langle i|b\rangle \, |i\rangle \otimes |j\rangle = |b\rangle \otimes |a\rangle$$
となり、確かに2つの量子状態を入れ替える演算子になっていることが分かります。
この表式を使って、$(O \otimes I) U_{\mathrm{SWAP}} (\rho \otimes \rho)$ を計算します。
$$(O \otimes I) U_{\mathrm{SWAP}} (\rho \otimes \rho) = \sum_{i,j} \left( O|i\rangle\langle j| \right) \otimes \left( |j\rangle\langle i| \right) (\rho \otimes \rho) = \sum_{i,j} \left( O|i\rangle\langle j|\rho \right) \otimes \left( |j\rangle\langle i|\rho \right)$$
この両辺のトレースを取ります。テンソル積のトレースは、各系のトレースの積に分解できることを使うと、
$$\mathrm{Tr}\left[ (O \otimes I) U_{\mathrm{SWAP}} (\rho \otimes \rho) \right] = \sum_{i,j} \mathrm{Tr}\left[ O|i\rangle\langle j|\rho \right] \cdot \mathrm{Tr}\left[ |j\rangle\langle i|\rho \right]$$
となります。ここで、それぞれのトレースは
$$\mathrm{Tr}\left[ O|i\rangle\langle j|\rho \right] = \langle j|\rho O|i\rangle, \qquad \mathrm{Tr}\left[ |j\rangle\langle i|\rho \right] = \langle i|\rho|j\rangle$$
と書けるので、
$$\sum_{i,j} \langle j|\rho O|i\rangle \langle i|\rho|j\rangle = \sum_j \langle j| \rho O \left( \sum_i |i\rangle\langle i| \right) \rho |j\rangle$$
となります。ここで、$\sum_i |i\rangle\langle i| = I$であることを使うと、
$$\sum_j \langle j| \rho O \rho |j\rangle = \mathrm{Tr}[\rho O \rho]$$
となります。最後に、トレースの巡回性($\mathrm{Tr}[ABC] = \mathrm{Tr}[BCA]$)を使うと、
$$\mathrm{Tr}[\rho O \rho] = \mathrm{Tr}[O \rho^2]$$
となり、以上より
$$\mathrm{Tr}[O \rho^2] = \mathrm{Tr}\left[ (O \otimes I)U_{\mathrm{SWAP}}(\rho \otimes \rho) \right] $$
従来手法：アダマールテスト
従来手法では、この量をアダマールテストで求めていました。アダマールテストは、補助量子ビットを用意し、注目系に対する制御ユニタリゲートを作用させたのち、補助量子ビットを $Z$ 測定することで、注目系に対するユニタリ演算の期待値の実部を求める量子アルゴリズムです。
アダマールテストを実装する量子回路は次のとおりです。
図3：アダマールテストの量子回路
(画像出典：Quantum Native Dojo, 2-2. アダマールテスト, https://dojo.qulacs.org/ja/latest/notebooks/2.2_Hadamard_test.html)
最初に、補助系にアダマールゲートを作用させた後、制御Uゲートを作用させます(Uは注目系に対して期待値を求めたいユニタリ演算子)。その後、さらにアダマールゲートを作用させる量子回路です。この時点での(補助系+注目系の)状態を $|\psi_H\rangle$ とします。このとき、次の等式が成り立ちます。
$$\langle \psi_H | Z_0 | \psi_H \rangle = \text{Re}\langle \psi | U | \psi \rangle \tag{8}$$
左辺は、補助系+注目系の状態 $|\psi_H\rangle$ に対する補助系の $Z$ 期待値($Z_0$ 期待値)です。一方、右辺は、注目系の状態 $|\psi\rangle$ に対する $U$ の期待値の実部です。従って、上記の量子回路を通した後に、補助系の $Z$ 測定を繰り返し行って得た測定値から平均値を作ると、それは、注目系に対する $U$ の期待値の実部になります。
アダマールテストの量子回路で密度演算子の期待値 $\mathrm{Tr}\left[ (O \otimes I) U_{\mathrm{SWAP}} (\rho \otimes \rho) \right]$ を求める際は、注目系の初期状態を \(\rho\otimes\rho\) とし、\(U=(O\otimes I)U_{\mathrm{SWAP}}\) とします。ただし、\(O\) はユニタリ演算子である必要があります。 分母は \(O=I\) とすれば計算できます。一方、分子については、一般に \(H_p\) がユニタリではないため、\(H_p\) をパウリ演算子の線形結合に分解し、各項についてアダマールテストを行うことで計算できます。
(8)式の証明
まず、量子回路の各ステップを追って $|\psi_H\rangle$ を具体的に求めます。補助系の初期状態を $|0\rangle$、注目系の初期状態を $|\psi\rangle$ とします。
ステップ1:補助系にアダマールゲートを作用
$$|0\rangle|\psi\rangle \;\longrightarrow\; \frac{1}{\sqrt{2}}\left(|0\rangle+|1\rangle\right)|\psi\rangle$$
ステップ2:補助系を制御ビットとした制御$U$ゲートを作用
補助系が $|1\rangle$ のときだけ注目系に $U$ が作用するので、
$$\frac{1}{\sqrt{2}}\left(|0\rangle+|1\rangle\right)|\psi\rangle \;\longrightarrow\; \frac{1}{\sqrt{2}}\left(|0\rangle|\psi\rangle+|1\rangle U|\psi\rangle\right)$$
ステップ3:補助系に再びアダマールゲートを作用
$H|0\rangle=\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$、$H|1\rangle=\frac{1}{\sqrt{2}}(|0\rangle-|1\rangle)$ を使うと、
$$\frac{1}{\sqrt{2}}\left(|0\rangle|\psi\rangle+|1\rangle U|\psi\rangle\right) \;\longrightarrow\; \frac{1}{2}|0\rangle\left(|\psi\rangle+U|\psi\rangle\right)+\frac{1}{2}|1\rangle\left(|\psi\rangle-U|\psi\rangle\right)$$
これが $|\psi_H\rangle$ です。
$$|\psi_H\rangle = \frac{1}{2}|0\rangle\left(|\psi\rangle+U|\psi\rangle\right)+\frac{1}{2}|1\rangle\left(|\psi\rangle-U|\psi\rangle\right)$$
$Z_0$期待値の計算
$Z_0$ は補助系にのみ作用するパウリ$Z$演算子で、$Z|0\rangle=|0\rangle$、$Z|1\rangle=-|1\rangle$ です。したがって、
$$Z_0|\psi_H\rangle = \frac{1}{2}|0\rangle\left(|\psi\rangle+U|\psi\rangle\right)-\frac{1}{2}|1\rangle\left(|\psi\rangle-U|\psi\rangle\right)$$
これと $\langle \psi_H|$ の内積を取ります。補助系の $|0\rangle,|1\rangle$ は直交するので、
$$\langle \psi_H | Z_0 | \psi_H \rangle = \frac{1}{4}\Big(\langle\psi|+\langle\psi|U^\dagger\Big)\Big(|\psi\rangle+U|\psi\rangle\Big) – \frac{1}{4}\Big(\langle\psi|-\langle\psi|U^\dagger\Big)\Big(|\psi\rangle-U|\psi\rangle\Big)$$
右辺の各項を展開すると、
$$\frac{1}{4}\Big(\langle\psi|\psi\rangle+\langle\psi|U|\psi\rangle+\langle\psi|U^\dagger|\psi\rangle+\langle\psi|U^\dagger U|\psi\rangle\Big) – \frac{1}{4}\Big(\langle\psi|\psi\rangle-\langle\psi|U|\psi\rangle-\langle\psi|U^\dagger|\psi\rangle+\langle\psi|U^\dagger U|\psi\rangle\Big)$$
$\langle\psi|\psi\rangle$ の項と $\langle\psi|U^\dagger U|\psi\rangle$ の項は打ち消し合い、残るのは
$$\langle \psi_H | Z_0 | \psi_H \rangle = \frac{1}{2}\Big(\langle\psi|U|\psi\rangle+\langle\psi|U^\dagger|\psi\rangle\Big)$$
となります。ここで、$\langle\psi|U^\dagger|\psi\rangle
=\langle\psi|U|\psi\rangle^*$なので、 $\langle\psi|U^\dagger|\psi\rangle=\langle\psi|U|\psi\rangle$ となり、
$$
\begin{aligned}
\langle\psi_H|Z_0|\psi_H\rangle
&=
\frac{1}{2}
\left(
\langle\psi|U|\psi\rangle
+
\langle\psi|U|\psi\rangle^*
\right)\\
&=
\operatorname{Re}\langle\psi|U|\psi\rangle
\end{aligned}
$$
しかし、この方法には次の課題があります。
- 1回の回路実行で、同じ量子状態を2つ同時に準備する( $\rho \otimes \rho$ )必要がある
- 複数の量子ビットにまたがるゲート(複合量子ビットゲート)が必要である
- 分母と分子を別々の回路で推定する必要がある
特に複合量子ビットゲートは、そのようなゲート操作をすることができないアニーリングマシンでは実装が難しく、これがアニーリングマシンへの仮想蒸留法の適用を妨げていました。
提案手法：古典シャドウ法
そこでこの論文が提案するのが古典シャドウ法を用いた仮想蒸留法の実装です。古典シャドウ法は、各量子ビットに対してランダムに $X, Y, Z$ 測定を行い、その結果を古典的に後処理することで期待値を推定する手法です[5]。アダマールテストと比べて、次の利点があります。
- 1回の測定で $\rho$ のみを準備すればよい($\rho \otimes \rho$ は不要）
- $X, Y, Z$ 測定だけでよく、複合量子ビットゲートが不要である
- 分母と分子を同じ測定結果からそれぞれ推定できる
$Z$測定から$X, Y$ 測定を実現する方法
まず、$X, Y$ 測定を実現する方法を確認します。測定直前に適切なユニタリ演算子 $U$ を作用させてから $Z$ 測定を行えば、$X, Y$ 測定を実現できます。パウリ演算子 $P \in \{X, Y, Z\}$ に対して $P = U^\dagger Z U$ が成り立つように $U$ を選びます。
$$
\tag{9}
P = Z \to U = I, \quad P = X \to U = H, \quad P = Y \to U = H S^\dagger
$$
このとき、次が成り立ちます。
$$
\tag{10}
\mathrm{Tr}[P \rho] = \mathrm{Tr}[U^\dagger Z U \rho] = \mathrm{Tr}[Z U \rho U^\dagger]
$$
すなわち、$U \rho U^\dagger$ に対する $Z$ 測定は、元の状態に対する $P$ 測定と等価になります。ただし、実際のアニーリングマシンでこの方法を用いるには、測定前に適切な単一量子ビット操作を行えるか、各パウリ基底で直接測定できることが必要です。
1量子ビット系での期待値 $\mathrm{Tr}[O \rho]$ の推定
次に、1量子ビット系での期待値 $\mathrm{Tr}[O \rho]$ の推定手順を示します。
- $U \in {I, H, H S^\dagger}$ を一様ランダムに選ぶ
- $U \rho U^\dagger$ を $Z$ 測定し、測定結果 $\hat{b} \in \{0, 1\}$ を得る
- 古典シャドウ $\hat{\rho} = 3 U^\dagger |\hat{b}\rangle\langle\hat{b}| U – I$ を計算する
この1〜3を独立に $M$ 回繰り返して平均をとると、期待値が推定できます。
$$
\tag{11}
\mathrm{Tr}[O \rho] \approx \frac{1}{M} \sum_{i=1}^{M} \mathrm{Tr}[O \hat{\rho}^{(i)}], \quad \hat{\rho}^{(i)} := i\text{ 回目に計算した古典シャドウ}
$$
上記の手順で期待値を推定できる理由
古典シャドウ $\hat{\rho} = 3U^\dagger|\hat{b}\rangle\langle\hat{b}|U – I$ の期待値が元の状態 $\rho$ となることを示せば良いです。
測定チャンネルの導出
$U$ を $\{I, H, HS^\dagger\}$ から一様ランダムに選び（それぞれ$Z, X, Y$測定が対応する）、測定値 $\hat{b} \in \{0,1\}$ を得た後、$U^\dagger|\hat{b}\rangle\langle \hat{b}|U$ を計算します。この演算子の期待値を$\rho$ を入力とする量子チャネル $\mathcal{M}$ の出力として定義します。
$U$ を選んで $Z$ 測定で $\hat{b}$ が出る確率は $\langle \hat{b}|U\rho U^\dagger|\hat{b}\rangle$ で、あるため、
$$\mathcal{M}(\rho) := \frac{1}{3}\sum_{U \in \{I,H,HS^\dagger\}} \sum_{b=0}^{1} \langle b|U\rho U^\dagger|b\rangle \, U^\dagger|b\rangle\langle b|U$$
1量子ビットの状態をパウリ基底で $\rho = \frac{1}{2}\left(I + \sum_{P \in \{X,Y,Z\}} r_P P\right)$、$r_P := \mathrm{Tr}[P\rho]$ と展開します。上式に代入して計算すると、
$$\mathcal{M}(\rho) = \frac{1}{2}\left(I + \frac{1}{3}\sum_P r_P P\right) = \frac{1}{3}\rho + \frac{1}{3}\mathrm{Tr}[\rho]\, I$$
となります。
逆写像をかける
$\mathcal{M}$ は可逆で、逆写像は
$$\mathcal{M}^{-1}(X) = 3X – \mathrm{Tr}[X]\,I$$
と表すことができます。ここで $\mathcal{M}^{-1}$ は物理的なチャネルではありませんが、古典計算として適用する分には問題ありません。
そこで、測定で得た $U^\dagger|\hat{b}\rangle\langle \hat{b}|U$ に $\mathcal{M}^{-1}$ を作用させたものを古典シャドウと定義します。
$$\hat{\rho} := \mathcal{M}^{-1}\left(U^\dagger|\hat{b}\rangle\langle \hat{b}|U\right) = 3U^\dagger|\hat{b}\rangle\langle \hat{b}|U – I$$
不偏推定量になること
$U$ と $\hat{b}$ に関する期待値をとると、$\mathcal{M}^{-1}$ の線形性から
$$\mathbb{E}[\hat{\rho}] = \mathcal{M}^{-1}\left(\mathbb{E}\left[U^\dagger|\hat{b}\rangle\langle \hat{b}|U\right]\right) = \mathcal{M}^{-1}(\mathcal{M}(\rho)) = \rho$$
つまり古典シャドウ $\hat{\rho}$ は $\rho$ の不偏推定量です。よって任意の演算子 $O$ に対して
$$\mathbb{E}\left[\mathrm{Tr}[O\hat{\rho}]\right] = \mathrm{Tr}[O\,\mathbb{E}[\hat{\rho}]] = \mathrm{Tr}[O\rho]$$
が成り立ち、$M$ 回繰り返して標本平均をとれば大数の法則により (11) 式のように $\mathrm{Tr}[O\rho]$ が推定できます。
なお、$\hat{\rho}$ 自体は負の固有値を持ちうるため物理的な密度行列ではありませんが、期待値の推定には支障ありません。
N量子ビット系での期待値 $\mathrm{Tr}[O \rho]$ の推定
$N$ 量子ビット系へは、上の手続きを各量子ビットに独立に適用するだけで拡張できます。1回の測定から得られる全系の古典シャドウは、各量子ビットの古典シャドウのテンソル積になります。
$$
\tag{12}
\hat{\rho} = \bigotimes_{j=1}^{N} \hat{\rho}_j = \bigotimes_{j=1}^{N} \left( 3 U_j^\dagger |\hat{b}_j\rangle\langle\hat{b}_j| U_j – I \right)
$$
ここで $\hat{\rho}_j$ は第 $j$ 量子ビットの古典シャドウです。$\mathrm{Tr}[O \hat{\rho}]$ の計算は古典コンピュータで行いますが、古典シャドウは$\hat{\rho}= \bigotimes_{j=1}^{N} \hat{\rho}_j $と局所的な形で表せ、$\text{Tr}[\hat{\rho_j}]=1$なので、\(O\) が局所的であれば、その期待値の計算量は主に \(O\) が作用する量子ビット数に依存し、全系の量子ビット数 \(N\) には直接依存しません。
N量子ビット系での期待値 $\mathrm{Tr}[O(\rho \otimes \rho)]$ の推定
最後に、仮想蒸留法に必要な $\mathrm{Tr}[O(\rho \otimes \rho)]$ を推定します。これも同じ手順で独立な測定から得られた古典シャドウの集合を作り、異なる2つの古典シャドウを組み合わせて次のように平均をとります。
$$
\tag{13}
\mathrm{Tr}[O(\rho\otimes\rho)]
\approx
\frac{1}{M(M-1)}
\sum_{\substack{i,j=1\\i\neq j}}^{M}
\mathrm{Tr}\left[
O\left(
\hat{\rho}^{(i)}\otimes\hat{\rho}^{(j)}
\right)
\right]
$$
同じ古典シャドウ集合から、 $O = U_{\mathrm{SWAP}} $とすれば分母、 $O = (H_p \otimes I) U_{\mathrm{SWAP}}$ とすれば分子が推定できます。すなわち、\(Z\) 測定に加えて \(X,Y\) 測定が可能なアニーリングマシンを用いれば、その測定結果に古典シャドウ法を適用することで、コピー間の量子ゲート操作を行わずにアニーリングマシン上で仮想蒸留法を実行できることになります。
（参考）局所仮想蒸留法（LVP）
この論文では、仮想蒸留法をより効率的にする局所仮想蒸留法（Local Virtual Purification; LVP）も用いられています。ここでは参考として要点のみ紹介します。
FVPには、アダマールテストを用いた場合に全量子ビットにまたがるゲート操作が必要になることと、推定値の分散が大きくなりうることという2つの問題がありました。LVPは、期待値を求めたい物理量が局所的である場合に使える、測定コストを削減できる可能性がある仮想蒸留法です[6]。
手順は次のとおりです。
局所的な物理量 $P_a$ に対して、$P_a$ が作用する量子ビットの領域 $A$、それを囲む領域 $B$（大きさを調整できる）、その他の領域 $C$ を設定します。領域 $C$ を部分トレースで無視し、$\rho_{A+B} := \mathrm{Tr}_C[\rho]$ に対してのみ仮想蒸留法を適用します。
量子ビットのトポロジが二次元正方格子の場合を例にとります。量子ビットがノードで、相互作用する量子ビットの組をエッジが繋いでいます。
図4：LVPにおける領域の設定例（二次元正方格子）
オレンジの領域が領域A、青の領域が領域B、うす緑の領域が領域Cです。
領域を局所に限ることで、アダマールテストを用いる場合は大域的なゲートが不要になること、推定値の分散を小さくできる場合があることという利点が得られます。LVPとFVPのずれは、適切な条件下では、領域 \(B\) の幅に対して指数的に小さくなることが知られており、\(B\) を十分大きくとればLVPによる期待値がFVPによる期待値に近づきます。古典シャドウ法を用いてエラー抑制する場合でも、LVPの方がFVPよりサンプリングコストを小さくできる場合があります。
数値実験
ここでは、原論文で行われた数値シミュレーションの結果を紹介します。
数値実験の設定
問題ハミルトニアンには、周期境界条件下の1次元スピン1/2 XXZ模型に磁場項を加えたものを用います。
$$
H_p = J\sum_{i=1}^{N}\left(X_i X_{i+1} + Y_i Y_{i+1} + \Delta Z_i Z_{i+1}\right) + h\sum_{i=1}^{N} Y_{i}\tag{14}
$$
パラメータは $J = -1$, $h = 1$, $\Delta = -0.73$ です。ドライバーハミルトニアンには次のランダム横磁場を用います。
$$H_d = -\sum_i h_i \sigma_i^x \tag{15}$$
ここで $h_i$ は $[0, 1]$ の一様分布から生成した定数です。アニーリングスケジュールは $A(t) = t/T$, $B(t) = 1 – t/T$ とし、初期状態は $H_d$ の基底状態とします。
時間発展
開放系の時間発展は、次のGKSLマスター方程式に従います。
$$\frac{d\rho}{dt} = -i[H(t), \rho] – \frac{\lambda}{2} \sum_{i=1}^{N} \sum_{\{P_i \in X_i, Y_i, Z_i\}} [P_i, [P_i, \rho]], \qquad \lambda = 0.0025 \tag{16}$$
第2項が各サイトにかかる1量子ビットのデポラライジングノイズを表しています。この方程式に従って、アニーリング時間 $T$ まで時間発展させます。得られた密度行列に対して $n = 2$ コピーの仮想蒸留を適用します。時間発展の数値計算にはQuTiPが用いられています。
結果
以下は横軸をアニーリング時間$T$、縦軸を厳密な基底エネルギー $E_g$ に対するエネルギー期待値の相対誤差 $|\langle H_p \rangle – E_g| / |E_g|$としてプロットしたグラフです。w/o QEMが元の状態の期待値に対する結果で、FVPが純粋化状態に対する期待値の結果です。LVPについては、本記事では詳細を説明しないため、参考として併記するに留めます。白丸はLVPによる期待値を古典シャドウ法で推定した結果で、エラーバーは標準誤差を表しています。
図5：アニーリング時間に対する相対誤差
グラフから、次のことが読み取れます。
- \(T\) が小さいうちは、非断熱遷移の影響が大きいため、いずれの場合も誤差が大きい
- \(T\) を大きくすると誤差は一度減少するが、ある \(T\) を境に再び増加に転じる。これは \(T\) を大きくするとインコヒーレントノイズの影響が蓄積するためと考えられ、単純に $T\to\infty$とすればよいわけではない
- 図に示したすべての \(T\) の範囲で、純粋化状態の期待値は元の状態の期待値よりも基底エネルギーに近い
今回の数値実験では、純粋化状態の期待値は、元の状態の期待値よりも、基底エネルギーに近くなりました。誤差が最小となる最適な $T$ が存在し、$T$ を大きくすれば真の基底エネルギーに漸近する、というわけではない点が特徴です。
まとめ
- 古典シャドウ法により、\(X,Y,Z\) 測定が可能なアニーリングマシンでは、コピー間の量子ゲート操作を用いずに仮想蒸留法を実装できることが示された
- 原論文の数値実験で、純粋化状態に対するエネルギー期待値が元の状態に対するエネルギー期待値よりも基底エネルギーに近づくことを確認した
- 今後の課題として、実際のアニーリングマシンでの実装が挙げられる
※ 参考：元論文では、局所仮想蒸留法(LVP)と古典シャドウ法を組み合わせる手法が提案されており、図中のLVPの曲線および白丸はその結果です。本記事では仮想蒸留法と古典シャドウ法の基本的な仕組みの紹介に主眼を置くため、詳細には立ち入りません。
参考文献
[紹介論文] Yoshida Takaharu, et al. “Hardware-efficient quantum annealing with error mitigation via classical shadow.” (2026).
[1] Kadowaki Tadashi, and Hidetoshi Nishimori. “Quantum annealing in the transverse Ising model.” Physical Review E 58.5 (1998): 5355.
[2] 遠藤傑. “量子エラー抑制とその周辺分野.” 物性若手夏の学校テキスト 第 69 回物性若手夏の学校. 物性若手夏の学校準備局, 2025.
[3] Huggins William J., et al. “Virtual distillation for quantum error mitigation.” Physical Review X 11.4 (2021): 041036.
[4] Koczor Bálint. “Exponential error suppression for near-term quantum devices.” Physical Review X 11.3 (2021): 031057.
[5] Huang Hsin-Yuan, Richard Kueng, and John Preskill. “Predicting many properties of a quantum system from very few measurements.” Nature Physics 16.10 (2020): 1050-1057.
[6] Hakoshima Hideaki, et al. “Localized virtual purification.” Physical Review Letters 133.8 (2024): 080601.
あとがき
量子アニーリングに対して、従来は主にゲート型量子コンピュータでの実装が想定されていた仮想蒸留法を適用できるというのは興味深いです。この論文で仮定していた、量子アニーリングにおいて \(Z\) 測定のほかに \(X,Y\) 測定もできるということが現実になれば、今回の量子エラー抑制の適用以外にも、より様々なことが量子アニーリングで行えるようになると思います。
東北大学大学院 情報科学研究科 博士課程前期