$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$

ボゾンサンプリングとは

例によって、主に公式文献のp.18に従って学んでいきます。
もう少し詳細な文献だとこれとか良さそうです。

ボゾンサンプリングはパッシブな線形光学素子だけで構成可能な、比較的シンプルな量子計算機です。
ただしそれゆえにユニバーサルな量子計算機にはなり得ないです。
では何が嬉しいかと言うと、ボゾンサンプリングはPermanent関数の結果を吐き出します。
ただし現実的に有用かというと実はそうでもないようです。

そこで今回は、このような流れで書いていこうと思います。

Permanent関数とは
ボゾンサンプリングの実装
ボゾンサンプリングの解析的計算
実際のところ役に立つの？

Permanent関数とは

Permanent関数は次のような数式で書けます。

$$
\mathrm{Perm(A)}= \sum_{\sigma \in S_n} \prod_{i=1}^{n} a_{i,\sigma (i)}
$$

$A$は$a_{i,j}$を要素にもつ$n \times n$行列で、$S_n$はn個の数字で可能な順列の集合です。
例えば、
$S_3 = $ { $(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1)$ }
です。

簡単に$A=2\times 2$行列の場合を考えると、
$\mathrm{Perm(A)}=a_{11} a_{22} + a_{12} a_{21}$
です。
$n=3$以降も計算するとわかるのですが、行列式と非常に似ています。
ただし行列式は各項の符号が偶置換or奇置換で変わりますが、Permanentは符号が変わりません。

応用例として、n人でランダムにプレゼント交換する場合に、誰も自分が持ってきたプレゼントを引き当てないような組み合わせの数を求めることができるようです。(リンク先の"derangement problem")

Permanent関数は古典計算だと最も効率的なアルゴリズムでも $O(2^{n-1} n^2)$ の計算量が必要だそうで、nについて指数関数的な計算リソースが必要です。よってもしこれを量子計算で効率的に求められれば量子計算の優位性を示すことができます。

ボゾンサンプリングの実装

ここはほとんど公式文献のまま、Jupiter notebook 実装ですが、入力光子数状態をランダム生成できるようにコードを追加しています(コメントアウトしてますが)。

import numpy as np
import strawberryfields as sf
import random as rd
from strawberryfields.ops import *
from numpy.linalg import multi_dot
from scipy.linalg import block_diag
import math

NodeNum = 4
SumOfPhoton = 3

# 入力光子数状態(FockList_ini)をランダムに生成
'''
FockList_ini = np.zeros(4, dtype=int)
for i in range(SumOfPhoton):
    tmp = rd.randint(0, NodeNum-1)
    FockList_ini[tmp] += int(1)
'''
    
# 1つのポートに2光子以上存在しない場合
'''
arr = np.arange(NodeNum)
np.random.shuffle(arr)
for i in range(SumOfPhoton):
    FockList_ini[arr[i]] += int(1)
'''

FockList_ini = [1,1,0,1]

eng, q = sf.Engine(4)

with eng:
    # prepare the input fock states
    Fock(FockList_ini[0]) | q[0]
    Fock(FockList_ini[1]) | q[1]
    Fock(FockList_ini[2]) | q[2] # vacuum state (optional)
    Fock(FockList_ini[3]) | q[3]

    # rotation gates
    Rgate(0.5719)
    Rgate(-1.9782)
    Rgate(2.0603)
    Rgate(0.0644)

    # beamsplitter array
    BSgate(0.7804, 0.8578)  | (q[0], q[1])
    BSgate(0.06406, 0.5165) | (q[2], q[3])
    BSgate(0.473, 0.1176)   | (q[1], q[2])
    BSgate(0.563, 0.1517)   | (q[0], q[1])
    BSgate(0.1323, 0.9946)  | (q[2], q[3])
    BSgate(0.311, 0.3231)   | (q[1], q[2])
    BSgate(0.4348, 0.0798)  | (q[0], q[1])
    BSgate(0.4368, 0.6157)  | (q[2], q[3])
    # end circuit
    # not performing measurement

# run the engine
state = eng.run('fock', cutoff_dim=8)

# extract the joint Fock probabilities
probs = state.all_fock_probs()

# print the joint Fock state probabilities
print(probs[1,1,0,1])

出力結果は、出力側光子検出器で検出されたフォトン数が{1,1,0,1}となる確率を示しています。
当然、光子数は入出力間で保存します。

0.17468916048563937

RゲートとBSゲートの結合は1つの$4 \times 4$ユニタリゲート$U$と等価で、入力ポート$i$のフォトンが出力ポート$j$から出てくる確率は$U_{ij}$で表せます。
この結果がPermanent関数と結びつく直感的な説明はこれのp.5がわかりやすいです。

例えば、入力ポート1,2に1つずつ存在したフォトンが、出力ポート3,4から1つずつ検出される確率は以下のように書けますね。
$$
U_{13} U_{24} + U_{14} U_{23}
$$
上の$2\times 2$行列のPermanent関数と見比べると、これは$4 \times 4$ユニタリゲート$U$を上手く整形してPermanent関数に渡してやれば出てきそうな気がしませんか？

ボゾンサンプリングの解析的計算

出力光子数状態${n_{1},n_{2},…,n_{N}}$が測定される確率は以下の通りです。
$$
\braket{ n_{1},n_{2},…,n_{N}}{ \psi^{'}} ^{2} =
\frac{|\mathrm{Perm}(U_{st})|^{2}}{m_{1}!m_{2}!,…,m_{N}!n_{1}!n_{2}!,…,n_{N}!}
$$

これを計算する実装例は以下の通りです。
元ネタとより詳細な説明はこちらになります。

def transform(photonNumList):
    resList = []
    for i in range(len(photonNumList)):
        for j in range(photonNumList[i]):
            resList.append(int(i))
    return resList

# Permanent関数を求める古典アルゴリズム
def perm(M):
     n = M.shape[0]
     d = np.ones(n)
     j =  0
     s = 1
     f = np.arange(n)
     v = M.sum(axis=0)
     p = np.prod(v)
     while (j < n-1):
         v -= 2*d[j]*M[j]
         d[j] = -d[j]
         s = -s
         prod = np.prod(v)
         p += s*prod
         f[0] = 0
         f[j] = f[j+1]
         f[j+1] = j+1
         j = f[0]
     return p/2**(n-1)

def Fact_Fock(FockList_ini, FockList_after):
    FockList_ini_Fact = [np.math.factorial(i) for i in FockList_ini]
    FockList_after_Fact = [np.math.factorial(i) for i in FockList_after]

    res = 1
    for i in range(len(FockList_ini)):
        res *= FockList_ini_Fact[i]

    for i in range(len(FockList_after)):
        res *= FockList_after_Fact[i]
    return res


# RゲートとBSゲートを結合し、1つのユニタリゲートで記述する
###################
Uphase = np.diag([np.exp(0.5719*1j),np.exp(-1.9782*1j),np.exp(2.0603*1j),np.exp(0.0644*1j)])

BSargs = [
(0.7804, 0.8578),
(0.06406, 0.5165),
(0.473, 0.1176),
(0.563, 0.1517),
(0.1323, 0.9946),
(0.311, 0.3231),
(0.4348, 0.0798),
(0.4368, 0.6157)
]

t_r_amplitudes = [(np.cos(q), np.exp(p*1j)*np.sin(q)) for q,p in BSargs]
BSunitaries = [np.array([[t, -np.conj(r)], [r, t]]) for t,r in t_r_amplitudes]
# list of 2x2 numpy array

UBS1 = block_diag(*BSunitaries[0:2])
UBS2 = block_diag([[1]], BSunitaries[2], [[1]])
UBS3 = block_diag(*BSunitaries[3:5])
UBS4 = block_diag([[1]], BSunitaries[5], [[1]])
UBS5 = block_diag(*BSunitaries[6:8])
U = multi_dot([UBS5, UBS4, UBS3, UBS2, UBS1, Uphase])
###################

# 出力側光子検出器の検出結果(FockList_after)をランダムに生成
###########
FockList_after = np.zeros(4, dtype=int)

for i in range(SumOfPhoton):
    tmp = rd.randint(0, NodeNum-1)
    FockList_after[tmp] += int(1)
###########

# 任意の値を使用する場合
FockList_after = [1,0,1,1]
  
print("Input photon number states")
print(list(FockList_ini))
print("Measured photon number states")
print(list(FockList_after), '\n')

# ボゾンサンプリングの出力を解析的に計算
div = Fact_Fock(FockList_ini, FockList_after)

ini_vec = transform(FockList_ini)
after_vec = transform(FockList_after)
clasic_calc = np.abs(perm(U[:,ini_vec][after_vec]))**2 / div
Boson_calc = probs[FockList_after[0],
                   FockList_after[1],
                   FockList_after[2],
                   FockList_after[3]]
print('Probability calculated with:')
print('Classical algorithm',clasic_calc)
print('Boson sampling',Boson_calc)

大まかな流れとして、

RゲートとBSゲートを結合し、1つのユニタリゲート$U$で記述する
光子が存在する入出力ポートに合わせて$U$を上式中$U_{st}$に整形し、Permanent関数を計算する。
解析的に得られた出力確率と、先ほど実装したボゾンサンプリング回路で得られた出力確率を比較する。
といったことを行っています。

出力結果はこのようになります。
比較の結果で言うと、非常によく一致します。
Strawberry Fieldsが非常によくできていることのデモンストレーションですね。

Input photon number states
[1, 1, 0, 1]
Measured photon number states
[1, 0, 1, 1] 

Probability calculated with:
Classical algorithm 0.046297023812036465
Boson sampling 0.046297023812036486

このように特定の入出力状態とそれが得られる確率がわかれば、$\mathrm{Perm}(U_{st})$を計算できます。

実際のところ役に立つの？

ここまで上手くいっていますが、ボゾンサンプリングには実用上の問題点があるようです。
まず、現実にボゾンサンプリングを実行しても上記コードのように直接出力確率が得られることはなく、あくまで1つの出力結果が得られるだけです。
よって、ある入出力状態について確率を知るには十分な回数のサンプリングを繰り返す必要があります。

ところが得られる出力光子数状態の組み合わせは光子数$n$に対してsuper exponentialです。
よって、サンプリングによって一定の精度を担保しようとすると結局指数関数的にリソースが増えかねないです。
また、そもそもどれだけサンプリングすればどれだけの精度を担保できるか自体が非常に難しい問題です。

入力光子数状態の生成にも課題があります。
各ポートには光子がほぼ同時に入力される必要がありますが、光子生成源があるtime bin $\delta t$ の間に光子を吐く確率が$p$であるような性質を持っている場合、$\delta t$の間にn個の光子が生成される確率は $p^n$ となってしまいます。
よって試行回数は光子数$n$の指数関数でスケールされてしまいます。

確率的でなく決定論的に単一光子を生成するのは技術的ハードルが高いようです。
あまり詳しくないですが。
ただし、入力光子数状態の問題は後に Gaussian Boson Sampling という手法が提案され回避可能になります。

過去の研究者の苦難の跡が見えますね。色々と勉強になります。
次は流れ的には Gaussian Boson Sampling ですが、似ているので飛ばすかもしれません。
読んでいただきありがとうございました！

参考にした文献

・公式ドキュメント
・https://en.wikipedia.org/wiki/Permanent_(mathematics)
・https://arxiv.org/abs/1406.6767