関数オラクル・QRAMを無理やり作る方法

# はじめに
量子アルゴリズムの中で関数の値をエンタングル状態で量子ビットに格納するという操作がよく現れる. つまり, 次のような操作である:
```math
\begin{align*}
    U_f\colon |j\rangle|0\rangle \mapsto |j\rangle|f(j)\rangle, \qquad j=0,1,\dots, 2^n-1, \tag{1}
\end{align*}
```
ただし, $X_n:=\\{0,1,\dots, 2^n-1\\}$ に対し, $f\colon X_n\to X_m$とする. この操作 $U_f$ は [QRAM](https://dojo.qulacs.org/ja/latest/notebooks/7.2c_quantum_random_access_memory.html) (_Quantum Random Access Memory_) と呼ばれている.

この $U_f$ は通常最初から存在が仮定されていることが多い. $U_f$ はオラクルであるので普通は構成する必要はないが, 少しシミュレーターでテストしてみたい等といった場合には無理やりにでも構成してみる必要が生ずるかも知れない[^1]. 

算術計算 (四則計算) については量子アルゴリズムが知られているので, 多項式計算を量子コンピュータ上で計算できることになる. そのため, 関数 $f$ を多項式で近似することにより, $f$ の手続きや定義式からオラクルとして $U_f$ を与えることは可能である. しかし, それを行うのは非常に大変だろうと思うので, 直接 $U_f$ を与えてしまおうというのが今回の目的である.

断っておくと, 量子アルゴリズムの中で $U_f$ がオラクルとして与えられているのは自然である. $U_f$ を実際に無理やり構成するにはすべての $j$ に関する $(j,f(j))$ の情報 (テーブル) を知る必要がある. 例えば, Grover の検索アルゴリズムで考えるならば, $U_f$ を構成しながらこの情報を同時に確認できるわけで, $U_f$ を作り終えた段階ですでに解は得られており, $U_f$ を使ってわざわざ検索アルゴリズムを実行する必要がそもそもない. したがって, 今回説明する内容については, 何かのアルゴリズム等のテストにおいて関数のオラクル $U_f$ が直接必要なときにしか実用上は役に立たないだろう[^2].

しかし, 一方で関数 $f$ のオラクル $U_f$ の存在をいつでも仮定できるかどうかを考えることには意味があるだろう. 結論を言えば, いかなる関数 $f\colon X_n\to X_m$ に対しても, (1) を満たすユニタリー $U_f$ は常に存在する[^3]. 本記事ではこの構成方法とその実装について述べる.


# 背景
数学的には単に Hilbert 空間の基底の取換えの話である. $\mathcal{H}$ を有限次元複素 Hilbert 空間とし, 正規直交基底 $\\{e_j\\}_j$, $\\{{f_j\\}}_j$ としよう. ただし内積 $\langle\cdot,\cdot\rangle$ は右線型・左反線型とする.

まず, 次の記法を導入する. $x,y\in\mathcal{H}$ に対して, 写像 $\mathcal{H}\times\mathcal{H}\ni z \mapsto \langle y,z\rangle x \in\mathcal{H}$ をここでは $x\hat{\otimes}y$ とかくことにする[^4]. $x\hat{\otimes}y$ は線型写像である. このとき, 次が成立する.

__命題.__ $\mathcal{H}$ 上の線型写像
```math
\begin{align*}
    U = \sum_j f_j\hat{\otimes}e_j
\end{align*}
```
はユニタリーである. また, すべての $j$ に対し, $Ue_j=f_j$ である.

__証明.__ $x,y,z,w\in\mathcal{H}$ に対し, 
```math
\begin{align*}
    & (x\hat{\otimes}y)^* = y\hat{\otimes}x, \\[.5em]
    & (x\hat{\otimes}y)(z\hat{\otimes}w)=\langle y,z\rangle x\hat{\otimes}w
\end{align*}
```
であることに注意すれば明らか. $\Box$


# 関数オラクルの構成[^5]
$n=1,2,\dots$ に対し, $\mathcal{H}_{n}:=\textrm{span}\\{|j\rangle\mid j\in X_n\\}$ とする. また, $f\colon X_n\to X_m$ に対し, $2^m\geq \max\_{j\in X_n} f(j)$ とする. テンソル積空間 $\mathcal{H}_n\otimes\mathcal{H}_m$ の正規直交基底は $\\{|j\rangle|k\rangle\\}\_{j\in X_n,k\in X_m}$ である.

前節を踏まえると, $j\in X_n$ に対して $|j\rangle|0\rangle$ と $|j\rangle|f(j)\rangle$ を入れ換え, その他はそのままとするとよい. つまり, 次のようになる:
```math
\begin{align*}
    & |0\rangle|0\rangle \overset{\text{入換え}}{\longleftrightarrow} |0\rangle|f(0)\rangle, \\
    & |1\rangle|0\rangle \overset{\text{入換え}}{\longleftrightarrow} |1\rangle|f(1)\rangle, \\
    & \hspace{2em}\vdots \\
    & |2^n-1\rangle|0\rangle \overset{\text{入換え}}{\longleftrightarrow} |2^n-1\rangle|f(2^n-1)\rangle, \\
    & |j\rangle|k\rangle \overset{\text{そのまま}}{\longleftrightarrow} |j\rangle|k\rangle, \qquad \text{その他}. \\
\end{align*}
```

実際は, $f(j)=0$ である $j$ については入れ換える必要はない. そこで, $\Gamma_f:=\\{(j,f(j))\in X_n\times X_m\mid j\in X_n\\}$ とする. すると $j\in X_n$,  $k\in X_m$ に対して次のように場合分けできる.
1. $(j,0)\in\Gamma_f$ のとき:
  $|j\rangle|k\rangle$ はそのまま.
1. $(j,0)\notin\Gamma_f$ かつ $(j,k)\in\Gamma_f$ のとき:
  $|j\rangle|k\rangle$ と $|j\rangle|0\rangle$ を入れ換える.
1. $(j,0)\notin\Gamma_f$ かつ $(j,k)\notin\Gamma_f$ のとき:
  $|j\rangle|k\rangle$ はそのまま. しかし, $|j\rangle|0\rangle$ は 2. で使われるので, カウントしない ($\Rightarrow k>0$).

以上をまとめると, $U_f$ が得られる.
```math
\begin{align*}
    U_f 
    &= \sum_{(j,0)\in\Gamma_f} |j\rangle\langle j|\otimes |k\rangle\langle k| \\
    &\hspace{2em} + \sum_{\substack{(j,0)\notin\Gamma_f\\ (j,k)\in\Gamma_f}} \Big\{ |j\rangle\langle j|\otimes |k\rangle\langle 0| + |j\rangle\langle j|\otimes |0\rangle\langle k| \Big\} + \sum_{\substack{(j,0)\notin\Gamma_f\\ (j,k)\notin\Gamma_f\\ k>0}} |j\rangle\langle j|\otimes |k\rangle\langle k|. \tag{2}
\end{align*}
```

__注意.__ 入れ換える必要がない $|j\rangle|k\rangle$ で関数 $f$ の値を計算するのに関係しないものは, 同じく入れ換える必要がない $|j\rangle|k\rangle$ と入れ換えてもよい. そのため, (1) を満たす $U_f$ は 1 つに定まらない.

__注意.__ $f(x,y)$ のような多変数関数の場合も考え方は変わらない. 単純に番号付けの問題であり, 2 変数関数の場合は $\mathcal{H}\_{n_1}\otimes\mathcal{H}\_{n_2}$ を $\mathcal{H}\_{n}$ と見ればよい. この場合については次節の実装において取り上げる.

### コード例
(2) の python コードは次のようになる. 多変数の場合にも対応した. 

$U_f$ を構成することはグラフ $\Gamma_f$ を構成することと同値であり, したがって, 引数はグラフ $\Gamma_f$ と変数の数を量子ビットに割り当てるリストである.

```python
import numpy as np

# U_f を作成する関数
def create_Uf(qubits_numbers, G_f):
    # qubits_numbers: |x1>|x2>…|xk>|f(x1,x2,…,xk)> に使う量子ビット数のサイズ k のリスト
    # G_f: f のグラフ Γ_f
    
    # |0>, |1>
    e = np.array([[1., 0.], [0., 1.]])
    
    H_n = [[i for i in range(2 ** n)] for n in qubits_numbers]
    H_sp = list(itertools.product(*H_n))

    def bit_base_vector(bit_string):
        # ビット列のベクトルを作成
        for i, bit in enumerate(bit_string):
            if i == 0:
                a = e[int(bit)]
            else:
                a = np.kron(a, e[int(bit)])
        return a

    Uf = np.zeros((2 ** (sum(qubits_numbers)), 2 ** (sum(qubits_numbers))))

    for v in H_sp:
        n = sum([v[i] * 2 ** (sum(qubits_numbers[i + 1:])) for i in range(len(v))])
        m = sum([v[i] * 2 ** (sum(qubits_numbers[i + 1:])) for i in range(len(v[:-1]))])
        bit_string_a = format(n, '0' + str(sum(qubits_numbers)) + 'b')
        bit_string_b = format(m, '0' + str(sum(qubits_numbers)) + 'b')

        # ビット列のベクトルを作成
        a = bit_base_vector(bit_string_a)
        b = bit_base_vector(bit_string_b)

        v_ = list(v[:-1])
        v_.append(0)
        v_ = tuple(v_)
        if v_ in G_f:
            # (j,0)∈Γ_f のとき
            Uf += np.outer(a, a)

        else:
            # (j,0)∈Γ_f でないとき
            if v not in G_f:
                # (j,k)∈Γ_f でないとき
                if v[-1] > 0:
                    # k>0 のときのみ
                    Uf += np.outer(a, a)
            else:
                # (j,k)∈Γ_f のとき
                Uf += np.outer(a, b) + np.outer(b, a)
                
    return Uf
```

# 数値例
ここでは $f(x,y)=[xy/2]$ を実装してシミュレーターで実行してみよう. ただし, $x\in\mathbb{R}$ に対して $[x]$ は $x$ を超えない最大の整数を表す. 量子ビット系は $\mathcal{H}_3\otimes\mathcal{H}_3\otimes\mathcal{H}_5$ とする.

今回シミュレーターは [Qulacs](https://dojo.qulacs.org/ja/latest/index.html) を用いる.

### 計算に必要なグラフを作成
$\Gamma_f=\\{(i,j,f(i,j))\mid (i,j)\in X_3\times X_3\\}$ の作成.
```python
import itertools

qubit_numbers = [3, 3, 5]    # 3+3+5量子ビット
X = [i for i in range(2 ** qubit_numbers[0])]
Y = [i for i in range(2 ** qubit_numbers[1])]

# X,Y のすべての組合せのリスト
X_Y = list(itertools.product(X, Y))
# print(X_Y)

# グラフ Γ_f の作成
G_f = {(i[0], i[1], (i[0] * i[1]) // 2) for i in X_Y}
print(G_f)
```
結果
> {(3, 4, 6), (6, 7, 21), (1, 0, 0), (7, 1, 3), (1, 5, 2), (1, 3, 1), (7, 5, 17), (6, 0, 0), (0, 5, 0), (1, 4, 2), (6, 2, 6), (1, 7, 3), (3, 6, 9), (7, 4, 14), (6, 5, 15), (0, 2, 0), (4, 0, 0), (2, 6, 6), (7, 7, 24), (2, 3, 3), (2, 7, 7), (4, 6, 12), (7, 6, 21), (6, 1, 3), (5, 7, 17), (3, 3, 4), (3, 2, 3), (6, 3, 9), (2, 0, 0), (0, 6, 0), (2, 4, 4), (1, 2, 1), (5, 1, 2), (3, 1, 1), (7, 3, 10), (4, 7, 14), (7, 0, 0), (0, 3, 0), (0, 0, 0), (4, 4, 8), (4, 2, 4), (6, 4, 12), (5, 6, 15), (4, 1, 2), (5, 0, 0), (1, 1, 0), (2, 2, 2), (6, 6, 18), (0, 7, 0), (5, 4, 10), (1, 6, 3), (2, 1, 1), (2, 5, 5), (3, 5, 7), (0, 4, 0), (5, 3, 7), (5, 5, 12), (3, 0, 0), (4, 3, 6), (7, 2, 7), (4, 5, 10), (5, 2, 5), (0, 1, 0), (3, 7, 10)}

### オラクル作成
$U_f$ を上述の関数 create_Uf を使って作成, ユニタリーの確認.
```python
U_f = create_Uf(qubit_numbers, G_f)
print(U_f)

# ユニタリーの確認
print(np.array_equal(U_f @ U_f.T, np.eye(2 ** sum(qubit_numbers))), np.array_equal(U_f.T @ U_f, np.eye(2 ** sum(qubit_numbers))))
```
結果
> [[1. 0. 0. ... 0. 0. 0.]
 [0. 1. 0. ... 0. 0. 0.]
 [0. 0. 1. ... 0. 0. 0.]
 ...
 [0. 0. 0. ... 1. 0. 0.]
 [0. 0. 0. ... 0. 1. 0.]
 [0. 0. 0. ... 0. 0. 1.]]
True True

### 計算
関数 $f(x,y)$ の計算.
```python
# 計算テスト
from qulacs import QuantumState
from qulacs import QuantumCircuit
from qulacs.gate import DenseMatrix

n_qubits = sum(qubit_numbers)
state = QuantumState(n_qubits)

# 入力 (x,y) のリスト
x = [4, 7]
print('[x,y]:', x)

# 入力 (x,y) と出力用のビット |0> を合わせたビット列に変換
x = sum([x[i] * 2 ** (sum(qubit_numbers[i + 1:])) for i in range(len(x))])
print('input: x= ', format(x, '011b'))
state.set_computational_basis(x)

circuit = QuantumCircuit(n_qubits)

gate = DenseMatrix([i for i in range(n_qubits)], U_f)
circuit.add_gate(gate)

circuit.update_quantum_state(state)
results = state.sampling(1)
f_value = results[0] % (2 ** qubit_numbers[2])
results = [format(result, '011b') for result in results]
print('result:', results)
print('f_value:', f_value)
```
結果
> [x,y]: [4, 7]
input: x=  10011100000
result: ['10011101110']
f_value: 14

"x= 10011100000", "result: ['10011101110']" の見方はそれぞれ $|100\rangle|111\rangle|00000\rangle$, $|100\rangle|111\rangle|01110\rangle$ である. すなわち, 
```math
\begin{align*}
    U_f|100\rangle|111\rangle|00000\rangle=|100\rangle|111\rangle|01110\rangle
\end{align*}
```
であることを示している. いろいろ $x,y$ の値を変えて試してもちゃんと計算されているようである.

# おわりに
この記事では, オラクル $U_f$ を実際に作って, シミュレーター上で計算させてみた. その結果, ちゃんと量子ビット上に計算結果が得られる様子を観察することができた.

さて最後に, 量子コンピュータといえば大人気アニメ「神様になった日」が思い起こされる. エンディングテーマ「Goodbye Seven Seas」は神曲なので是非聞いてみて欲しい.

# 参考文献
[1] Mitarai, K., [京大基研量子情報スクール用ノート](https://www2.yukawa.kyoto-u.ac.jp/~qischool2019/mitaraiCTO.pdf)
[2] Quantum Native Dojo, 7章, [コラム：量子ランダムアクセスメモリー (qRAM)](https://dojo.qulacs.org/ja/latest/notebooks/7.2c_quantum_random_access_memory.html)



[^1]: $|j\rangle|0\rangle$ にしか作用させないのであれば, シミュレーターでは $U_f$ は必ずしもユニタリーである必要はないかもしれない. QulacsだとDenseMatrixでユニタリーでない行列を回路に放り込めるようである (多分).
[^2]: このように, 量子コンピュータの優位性は必ずしも古典から直ちに得られるわけではない. 検索の話であれば, データベースが"量子"の枠組みの中で管理されている必要がある.
[^3]: 書籍等を見てもこの事実は語られていないような気がする. まあ当たり前と言われれば当たり前ではあるので, わざわざ言及していないだけかもしれない. また, このことと量子コンピュータ上に実際に実装できるかは別の話であり, ここでは単に数学的には必ず存在し, 仮定しても論理的に問題がないということである. この辺については参考文献 [1,2] によると, 現在研究が行われているテーマであるとのこと.
[^4]: ブラケット記法を用いれば, $|x\rangle\langle y|$ を表す.
[^5]: 構成してしまったものをはたしてオラクルと呼べるだろうか？