1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@shintaro_minagawa(皆川 慎太朗)

量子情報における半正定値計画問題 2-2. 量子状態識別問題のSDPの数値計算

1
Last updated at Posted at 2026-01-17

はじめに

量子状態識別問題がSDPとして定式化されていく流れを前回の記事で議論しました。

本記事では、このSDPをPythonで数値的に解く方法を解説していきます。

参考文献

以下の本がCVXPYを使ったPythonでの凸最適化の勉強にとってもなります。

寒野善博，駒木文保（編），『最適化手法入門』，講談社サイエンティフィク，2019．

復習

与えられたアンサンブル$\mathcal E^\mathcal X_\mathrm{A}:=\{p_X(x);\rho^x_\mathrm{A}\}$に対して、最適な量子状態識別の成功確率は以下のようなSDPで与えられるのでした。

\begin{align}
\underset{E^x_\mathrm A}{\text{Maximize}}\quad & \sum_{x\in\mathcal X}p_x\mathrm{Tr}[E^x_\mathrm A\rho^x_\mathrm A]\\
\text{subject to} \quad & \sum_{x\in\mathcal X}E^x_\mathrm{A}=I_\mathrm{A}\\
\quad &E^x_\mathrm{A}\succeq 0_\mathrm{A}\quad(\forall x\in\mathcal X).
\end{align}

ここからはこの問題をPythonを使って解いていこうと思います。

数値計算

パッケージ

計算に使うパッケージはNumPyとCVXPYです。本来、凸最適化のソルバーは解きたい問題を決まった形式に落とし込んだ上で使うのですが、CVXPYがそこを橋渡ししてくれるので、自然な形で問題を記述できるようになります。CVXPYの使い方は後で見ていきます。

import numpy as np
import cvxpy as cp
import matplotlib.pyplot as plt

有益な関数の準備

次に、いくつか関数を用意しておきます。

# 複素共役転置
def dagger(A: np.ndarray) -> np.ndarray:
    return A.conj().T

# ベクトル |v> を演算子 |v><v| に（vは正規化されている想定）
def projector(v: np.ndarray) -> np.ndarray:
    v = np.asarray(v, dtype=np.complex128).reshape(-1, 1)
    return v @ dagger(v)

# 行列をエルミート化
def hermitize(A: np.ndarray) -> np.ndarray:
    A = np.asarray(A, dtype=np.complex128)
    return 0.5 * (A + dagger(A))

dagger関数は、A.conj().Tで、行列を複素共役を取って転置します。
projectorは、正規化されたベクトル|v>に<v|と行列にします。reshape(-1, 1)は、NumPy配列を列は1列に固定、行数は自動で調整することで縦ベクトルに変換しています。
hermitize関数は

    \left(\frac{A+A^\dagger}{2}\right)^\dagger=\frac{A^\dagger+(A^\dagger)^\dagger}{2}=\frac{A+A^\dagger}{2}

なので、任意の複素正方行列$A$について$\frac{A+A^\dagger}{2}$がエルミート行列になることを使っています。

基本的にこれらの関数はすべて複素NumPy配列を複素NumPy配列に変換するという想定ですが、念の為、np.asarray(v, dtype=np.complex128)と加えて、listやtupleが来てもNumPy配列に変換します。

状態の定義

ここでは、$\mathcal H_\mathrm{A}$を3次元複素ヒルベルト空間として、$\{|0\rangle_\mathrm{A},|1\rangle_\mathrm{A},|2\rangle_\mathrm{A}\}$を標準基底とします。次の3状態を識別する問題を考えます。

    \begin{align}
        |\psi_0\rangle_\mathrm{A} &= |0\rangle_\mathrm{A}\\
    	|\psi_1\rangle_\mathrm{A} &= \cos\theta\,|0\rangle_\mathrm{A} + \sin\theta\,|1\rangle_\mathrm{A}\\
    	|\psi_2\rangle_\mathrm{A} &= \cos\theta\,|0\rangle_\mathrm{A} + \sin\theta\,|2\rangle_\mathrm{A}\\
    \end{align}

$\theta$が0のときは全部同じ状態なので識別しようがなく、$\theta=\pi/2$のときに全状態が標準基底となって直交するので、

\{|0\rangle\langle 0|_\mathrm{A},|1\rangle\langle 1|_\mathrm{A},|2\rangle\langle 2|_\mathrm{A}\}

という測定で完全に識別できます。

それでは、この状態を準備するコードを書いていきましょう。

def states(theta, priors):
    #タプルのpriorsをNumPy配列に変換
    p = np.array(priors, dtype=float)

    # 標準基底
    e0 = np.array([1, 0, 0], dtype=np.complex128)
    e1 = np.array([0, 1, 0], dtype=np.complex128)
    e2 = np.array([0, 0, 1], dtype=np.complex128)

    #状態の生成(ベクトル)
    psi0 = e0
    psi1 = np.cos(theta) * e0 + np.sin(theta) * e1
    psi2 = np.cos(theta) * e0 + np.sin(theta) * e2

    #行列として状態を保存
    rhos = [projector(psi0), projector(psi1), projector(psi2)]

    return rhos, p

priorsはタプルで、3状態の出現確率です。states関数はpriorsをNumPy配列にしたpと、NumPy配列としての3状態が入ったリストrhosを返します。

SDPの部分

ここが本題です。当初の問題のSDPをCVXPYで解いていきます。

def SDP(p, rhos):
    # 恒等演算子の生成
    I = np.eye(3, dtype=np.complex128)
    # 変数の定義（POVM要素）
    Es = []
    for i in range(3):
        E = cp.Variable((3, 3), complex=True)
        Es.append(E)

    # 制約条件
    constraints = []

    for E in Es:
        constraints += [E == E.H]  # エルミート性
        constraints += [E >> 0]    # 半正定値性

    constraints += [Es[0] + Es[1] + Es[2] == I]  # 完備性（和がI）

    # 目的関数
    obj = 0
    for i in range(3):
        obj += p[i] * cp.real(cp.trace(rhos[i] @ Es[i])) #cp.realで実数の式と明示

    # 問題の定義
    problem = cp.Problem(cp.Maximize(obj), constraints)

    # 最適化
    problem.solve(solver="SCS", eps=1e-6, max_iters=20000)

    # 最適値を返す
    return problem.value

まず、POVM要素を変数としてもつ本SDPの変数リストを作ります。3状態識別なので、リストの要素も3つです。各要素は、$3\times 3$複素行列です。cp.VariableがCVXPYでの変数の定義方法で、

cp.Variable((3, 3), complex=True)

として変数を生成して、i=0から2まで回してEsにappendしていきます。

次に、CVXPYでは制約条件もリストとして扱います。

constraints += [条件]

で条件を追加していきます。Eがエルミートという条件をE==E.H、Eが半正定値という条件はE>>0、POVM要素の和が恒等演算子というのはEs[0] + Es[1] + Es[2] == Iというような書き方は、CVXPYの表記です。

目的関数$\sum_{x\in\mathcal X}p_x\mathrm{Tr}[E^x_\mathrm A\rho^x_\mathrm A]$はその通りに書いていくのですが、2点。CVXPYにはトレースがあるので、cp.traceを使えます。また、この式の値は絶対に実数なのですが、$E^x_\mathrm A$や$\rho^x_\mathrm A$は元々複素行列だったので、明示的に「これは実数の式」というためにcp.realが重要です。

問題の定義の部分

problem = cp.Problem(cp.Maximize(obj), constraints)

は、目的関数objをconstraintsの制約条件のもとでmaximizeするという意味です。

最適化を解くのはここの部分です。

problem.solve(solver="SCS", eps=1e-6, max_iters=20000)

SCSはSplitting Conic Solverのことで、無料で使え、Google Colabでも動くソルバーです。何度も反復して最適に近づいていくという設計を取っているみたいで、eps=1e-6はこれ以上改善できないとしてストップする精度、max_iters=20000は最大2万回反復計算という意味です。

Θを動かしてプロット

$\theta$を$0$から$\pi$まで動かして、各$\theta$で最適識別成功確率を計算し、描画します。3状態の出現する事前確率はそれぞれ$0.4, 0.3, 0.3$としています。

def plot(priors):
    thetas = np.linspace(0.0, np.pi, 181)
    p_guess = []

    for k in range(len(thetas)):
        theta = thetas[k]
        rhos, p = states(theta, priors=priors)
        val = SDP(p, rhos)
        p_guess.append(float(val))

    plt.figure()
    plt.plot(thetas, p_guess)
    plt.xlabel(r"$\theta$")
    plt.ylabel(r"optimal guessing probability")
    plt.grid(True)
    plt.show()

if __name__ == "__main__":
    plot(priors=(0.4, 0.3, 0.3))

実行

まとめると、コードが以下のようになります。

import numpy as np
import cvxpy as cp
import matplotlib.pyplot as plt

# ===============
# ユーティリティ
# ===============
# 複素共役転置
def dagger(A: np.ndarray) -> np.ndarray:
    return A.conj().T

# ベクトル |v> を演算子 |v><v| に（vは正規化されている想定）
def projector(v: np.ndarray) -> np.ndarray:
    v = np.asarray(v, dtype=np.complex128).reshape(-1, 1)
    return v @ dagger(v)

# 行列をエルミート化
def hermitize(A: np.ndarray) -> np.ndarray:
    A = np.asarray(A, dtype=np.complex128)
    return 0.5 * (A + dagger(A))

# ===============
# 状態の生成
# ===============
def states(theta, priors):
    #タプルのpriorsをNumPy配列に変換
    p = np.array(priors, dtype=float)

    # 標準基底
    e0 = np.array([1, 0, 0], dtype=np.complex128)
    e1 = np.array([0, 1, 0], dtype=np.complex128)
    e2 = np.array([0, 0, 1], dtype=np.complex128)

    #状態の生成(ベクトル)
    psi0 = e0
    psi1 = np.cos(theta) * e0 + np.sin(theta) * e1
    psi2 = np.cos(theta) * e0 + np.sin(theta) * e2

    #行列として状態を保存
    rhos = [projector(psi0), projector(psi1), projector(psi2)]

    return rhos, p

# ===============
# 状態識別のSDP
# ===============
def SDP(p, rhos):
    I = np.eye(3, dtype=np.complex128)
    # 変数の定義（POVM要素）
    Es = []
    for i in range(3):
        E = cp.Variable((3, 3), complex=True)
        Es.append(E)

    # 制約条件
    constraints = []

    for E in Es:
        constraints += [E == E.H]  # エルミート性
        constraints += [E >> 0]    # 半正定値性

    constraints += [Es[0] + Es[1] + Es[2] == I]  # 完備性（和がI）

    # 目的関数
    obj = 0
    for i in range(3):
        obj += p[i] * cp.real(cp.trace(rhos[i] @ Es[i]))

    # 問題の定義
    problem = cp.Problem(cp.Maximize(obj), constraints)

    # 最適化
    problem.solve(solver="SCS", eps=1e-6, max_iters=20000)

    # 最適値を返す
    return problem.value

# ===============
# Θごとにプロット
# ===============
def plot(priors):
    thetas = np.linspace(0.0, np.pi, 181)
    p_guess = []

    for k in range(len(thetas)):
        theta = thetas[k]
        rhos, p = states(theta, priors=priors)
        val = SDP(p, rhos)
        p_guess.append(float(val))
        
    plt.figure()
    plt.plot(thetas, p_guess)
    plt.xlabel(r"$\theta$")
    plt.ylabel(r"optimal guessing probability")
    plt.grid(True)
    plt.show()

if __name__ == "__main__":
    plot(priors=(0.4, 0.3, 0.3))

実行するとこんな感じになります。

download-1.png

やはり$\pi/2$で識別成功確率が1になっていますね。また、$\theta=0$の場合はもはや識別不可能なので、常に一番出現確率が高い状態を言うのが最適となって$0.4$となっています。

1
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?