任意の論理関数を量子回路で表現してみる

$$
\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}}
$$

## はじめに

量子アルゴリズムの本を読むと、ドイッチュ・ジョザとかグローバーのあたりで「量子オラクル」なるものがでてきて、天下り的にこんな動作をする量子回路があったとする、という前提のもとで説明がなされたりしますが、具体的な関数があったときに、どうやったらそれを実現する量子回路が作れるのでしょうか。なかなかそれをシンプルに解説した記事が見当たらなかったのですが、以下の記事などを見ながら、自分なりにこうやれば良いということがわかってきたので、その手順を示してみたいと思います。そして、自作の量子計算シミュレータ[qlazy](https://github.com/samn33/qlazy)で、動作を確認します。

* [山下、松尾：「量子回路と最適化」](http://www.orsj.or.jp/archive2/or63-6/or63_6_342.pdf)
* [Multi-controlled NOTゲートを補助ビットなしに作ってみる(@gyu-don)](https://qiita.com/gyu-don/items/7c52c0d20f9998022803)

体系的にちゃんと勉強されている人にとっては、当たり前すぎる話かもしれませんし、もしかすると間違って理解しているところもあるかもしれませんが、そのあたりは大目に見ていただきながら(ご指摘いただけるとありがたいです)、以下読んでいただければと思います。

## 問題設定

ある論理関数$y=f(x_0,x_1,\cdots ,x_{n-1})$が与えられたとき、これを実現する量子回路を求めよ、という問題を考えます。すなわち、量子状態を以下のように変換するユニタリ演算子を基本量子ゲートの組み合わせだけで表現できれば、それが答えというわけです。

```math
\ket{x_0,x_1,\cdots ,x_{n-1}} \otimes \ket{y} \rightarrow \ket{x_0,x_1,\cdots ,x_{n-1}} \otimes \ket{y \oplus f(x_0,x_1,\cdots ,x_{n-1})}
```

論理関数は、真理値表で表すのがわかりやすいです。具体的に$n=3$として、例えば、

|$x_2$|$x_1$|$x_0$|$y$|
|-----|-----|-----|---|
|$0$  |$0$  |$0$  |$0$|
|$0$  |$0$  |$1$  |$1$|
|$0$  |$1$  |$0$  |$0$|
|$0$  |$1$  |$1$  |$0$|
|$1$  |$0$  |$0$  |$0$|
|$1$  |$0$  |$1$  |$1$|
|$1$  |$1$  |$0$  |$0$|
|$1$  |$1$  |$1$  |$0$|

のような真理値表(Trueを$1$、Falseを$0$としました)を題材にして、その設計手順を示してみます。ここで、$x_2$を最上位ビットとしたビット演算との対応のため、表の一番左の列を$x_2$としました。


## 設計手順

### 真理値表から論理式へ

まず、真理値表から論理式の形にします。どうするかというと、真理値表の中の出力がTrueになっている行に注目して、そのときの$x_0,x_1,x_2$の各値を見ます。Trueの場合はそのまま、Falseの場合は否定形にしてその論理積をとります。それを出力がTrueになっている行すべてについて繰り返し、(排他的)論理和をとります。以上です。とても簡単ですね。上の真理値表の場合、2行目と6行目がTrueになっているので、2つの論理積の(排他的)論理和として、

```math
f(x_0,x_1,x_2) = x_0 \bar{x_1} \bar{x_2} \oplus x_0 \bar{x_1} x_2  \tag{1}
```

のように書けます。本当にこれで良いかは簡単に確認できます。すべての{$x_0,x_1,x_2$}の組み合わせを上の式に代入すると、2行目、6行目以外は、全部$0$になります。


### 論理式からMPMCTへ

論理式から基本量子ゲートに行く前に、MPMCT(Mixed Polarity Multiple-Controll Toffoli)ゲートという表現を使って、一旦、中間的な論理設計という段階を踏みます。MPMCTゲートというのは、複数の正または負の極性をもった制御ビットと1つの標的ビットを持ち、正極性のビットがすべて1、かつ、負極性のビットがすべて0のときに、標的ビットの値を反転させるような演算を行う回路です。例えば、以下のように描かれたりします(ここで、塗り潰されている丸は「正極性」、塗り潰されていない丸を「負極性」と見てください)。

```
c0 --◯--
c1 --●--
c2 --●--
c3 --○--
t  --CX--
```

この例の場合、制御ビット$c_1$と制御ビット$c_2$が$1$、制御ビット$c_0$と制御ビット$c_3$が$0$のときのみ、標的ビット$t$が反転します。この表現を使うと、式(1)の論理関数は、以下のように表せます($y$の初期値は$0$とします)。

```
x0 --●--●--
x1 --○--○--
x2 --○--●--
y  --CX--CX--
```

### MPMCTからMulti-Controlled NOTを含む量子回路へ

MPMCTの●や○は、基本量子ゲートで置き換えられます。正極性の●は特に何もせず、通常の制御ビットを表すドットに置き換えます。負極性の○の方は、通常の制御ビットと動作を逆にしないといけないので、両側にXゲートを置きます。ということで、以下でOKです。

```
x0 -----*---------*-----
x1 --X--*--X---X--*--X--
x2 --X--*--X------*-----
y  -----CX--------CX----
```

(真中のX-Xは合わせるとキャンセルするので、本当は不要なのですが、とりあえず置いておいて、説明を先に進めます)

### Multi-Controlled NOTを基本量子ゲートで表現

ここで制御ビットが複数ある特殊な制御NOTゲートがでてきました。Multi-Controlled NOTとか、一般化制御NOTとか言われているようです。これを基本量子ゲートでどう表現するか？補助量子ビットを使わない優れたやり方を[@gyu-donさん](https://qiita.com/gyu-don/items/7c52c0d20f9998022803)が紹介されていますので、それを参考にしたいと思います。

基本方針は、Multi-Controlled P(位相シフトゲート)を作っておいて、Pの前後にアダマールをおきます。このときの回転角がちょうど$\pi$になるようにしておけば、HP($\pi$)H=HZH=Xなので、これでMulti-Controlled NOTが実現できます。ちなみに、Pの定義は、以下です。

```math
P(\theta) =
\begin{pmatrix}
1 & 0 \\
0 & e^{i \theta}
\end{pmatrix}
```

Mult-Controlled Pをどう作るかですが、[グレイコード](https://ja.m.wikipedia.org/wiki/%E3%82%B0%E3%83%AC%E3%82%A4%E3%82%B3%E3%83%BC%E3%83%89)を使ったとても巧妙なやり方になっています。ここではその手順だけをざっくり記載しておきます。

* [STEP.0] nビットの制御ビットと1ビットの標的ビットを設定し、グレイコードを0に初期化します。
* [STEP.1] グレイコードを1つ進めます。
* [STEP.2] グレイコードが0のときは[STEP.1]に戻ります。それ以外のときは[STEP.3]に進みます。
* [STEP.3] グレイコードを1つ進めたことで変化したビット番号を制御ビット、グレイコードの1が立っている最上位ビット番号を標的ビットとしたCNOTゲートをかけます。ただし、制御と標的が同じになってしまった場合、制御ビット番号の方を1つ下げます。また、グレイコードが1の場合、何もしません。
* [STEP.4] グレイコードの1が立っている最上位ビット番号を制御ビットにして、全体の標的ビットを標的にする制御P(位相シフト)をかけます。このときの回転角は、奇数番目のグレイコードの場合、$\phi = +\pi / 2^{n-1}$、偶数番目のグレイコードの場合、$\phi = -\pi / 2^{n-1}$とします。
* [STEP.5] グレイコードがnビットの範囲内にある場合、[STEP.1]に戻ります。それ以外となった場合、終了します。

一見したところ、とても不思議な手順なのですが、実際にこの通り回路を組んで手計算で確認すると、確かにMulti-Contorolled Pが実現されることが確認できます。が、なぜこれでうまくいくのか、正直きちんと理解できていません(汗)。わかったら補足・追記したいと思います。とりあえず、制御ビットが3つの場合の回路を示して、説明した(お茶をにごす)ことにします。


```
    (000) (001) (011)      (010)      (110)      (111)      (101)      (100)
x0 --------*-----*----------*---------------------*---------------------*----------
x1 --------|-----CX---*-----CX---*-----*----------|----------*----------|----------
x2 --------|----------|----------|-----CX---*-----CX---*-----CX---*-----CX---*-----
y  -------P(+)-------P(-)-------P(+)-------P(-)-------P(+)-------P(-)-------P(+)---
```

ここで、()内の3桁のビット列は3ビットのグレイコードです。左から右に１つずつ進めています。P(+)は位相$\phi = +\pi / 2^{n-1}$の制御Pゲート、P(-)は位相$\phi = -\pi / 2^{n-1}$の制御Pゲートを表しているものとします。これをMulti-Controlled NOTにするためには、標的ビットyの一番最初と最後にアダマールをかければ良いです。入力がnビットの場合も同様の議論で設計することが可能です。

以上で、与えられた論理関数に対する、量子回路が実現できたことになります。


### 複数出力をもつ論理関数への拡張

さて、ここまで説明してきた「問題設定」は、出力側のビット数が1ビットの場合でした。2ビット以上の場合どうなるの？ということが気になるかもしれませんが、簡単に拡張できます。例えば、

```
x0 --●--●----●--●--
x1 --○--○----○--○--
x2 --○--●----○--●--
y0 --CX--CX-----------
y1 -----------CX--CX--

```

こんな具合に、$y_0$に対する回路と$y_1$に対する回路をつなげれば良いです。


## シミュレータで確認

それではシミュレータで動作を確認します。上で示した真理値表を実装します。

### 実装

全体のPythonコードを以下に示します。

【2021.9.5追記】[qlazy](https://github.com/samn33/qlazy)最新版でのソースコードは[ここ](https://github.com/samn33/qlazy/blob/master/example/py/LogicalFunction)に置いてあります。

```python
from qlazypy import QState

# custom gates

def hadamard(self,id):

    for i in range(len(id)):
        self.h(id[i])
    return self

def multi_cx(self,id_in,id_out):

    #
    # hadamard
    #
    self.h(id_out[0])

    #
    # controlled-P(psi), psi=pi/(2**(bitnum-1))
    #
    bitnum = len(id_in)
    psi = 1.0/(2**(bitnum-1)) # unit=pi(radian)
    gray_pre = 0
    for gray in gray_code(bitnum):
        if gray == 0:
            continue
        msb = len(str(bin(gray)))-3
        chb = len(str(bin(gray^gray_pre)))-3
        if gray != 1:
            if chb == msb:
                chb -= 1
            self.cx(id_in[chb],id_in[msb])
        self.cp(id_in[msb],id_out[0],phase=psi)
        psi = -psi
        gray_pre = gray
    
    #
    # hadamard
    #
    self.h(id_out[0])

    return self
    
def mpmct(self,id_in,id_out,x):

    bitnum = len(id_in)
    
    for i in range(bitnum):
        bit = (x>>i)%2
        if bit == 0:
            self.x(id_in[i])
            
    self.multi_cx(id_in,id_out)

    for i in range(bitnum):
        bit = (x>>i)%2
        if bit == 0:
            self.x(id_in[i])

    return self

def q_oracle(self,id_in,id_out,func):

    num = 2**len(id_in)
    for x in range(num):
        y = func(x)
        if y == 1:
            self.mpmct(id_in,id_out,x)

    return self

# functions

def gray_code(n):

    for k in range(2**n):
        yield k^(k>>1)
        
def func(x):

    if x == 1 or x == 5:
        y = 1
    else:
        y = 0
    return y

def create_register(num):

    return [0]*num

def init_register(*args):

    idx = 0
    for i in range(len(args)):
        for j in range(len(args[i])):
            args[i][j] = idx
            idx += 1
    return idx

def result(self,id_a,id_b):

    # measurement
    id_ab = id_a + id_b
    iid_ab = id_ab[::-1]
    shots = (2**len(id_a))*5
    freq = self.m(id=iid_ab,shots=shots).frq
    
    # set results
    a_list = []
    r_list = []
    for i in range(len(freq)):
        if freq[i] > 0:
            a_list.append(i%(2**len(id_a)))
            r_list.append(i>>len(id_a))
    return (a_list,r_list)

if __name__ == '__main__':

    # add custom gates
    QState.hadamard = hadamard
    QState.multi_cx = multi_cx
    QState.mpmct = mpmct
    QState.q_oracle = q_oracle
    QState.result = result

    # set register
    id_in = create_register(3)
    id_out = create_register(1)
    qnum = init_register(id_in,id_out)

    # initial state
    qs = QState(qnum)
    qs.hadamard(id_in)

    # quantum oracle
    qs.q_oracle(id_in,id_out,func)

    # print results
    in_list,out_list = qs.result(id_in,id_out)
    for i in range(len(in_list)):
        print("{0:03b} -> {1:d}".format(in_list[i],out_list[i]))

    qs.free()
```

main部の処理を見ながら、何をやっているか簡単に説明します。


#### カスタムゲートの追加

```python
# add custom gates
QState.hadamard = hadamard
QState.multi_cx = multi_cx
QState.mpmct = mpmct
QState.q_oracle = q_oracle
QState.result = result
```

qlazypyのQStateクラスに自分で定義した量子ゲートをメソッドとして追加しています。関数定義はコードの上の方にあります。hadamardは、レジスタを構成する全量子ビットにアダマールをかけます。multi_cxは、Multi-Controlled NOTを実行します。mpmctは、与えられた{$x_0,x_1,\cdots,x_{n-1}$}の組に対応したMPMCTを実行します。q_oracleは、与えられた論理関数に対応した量子オラクルを実行します。resultは、最終的な量子状態を測定して、各レジスタの値のリストを返します。

#### レジスタの設定

```python
# set register
id_in = create_register(3)
id_out = create_register(1)
qnum = init_register(id_in,id_out)
```

制御ビットのレジスタ(入力)を3ビット、標的ビットのレジスタ(出力)を1ビットに設定して初期化し、全体で必要になる量子ビット数を取得します。id_in,id_outの実体は量子ビット番号を要素とする単なるリストです。create_registerですべての要素が0になっているリストを一旦作っておいて、init_registerで各レジスタに重なりがないように0から順番に量子ビット番号を入れています(関数定義を参照)。

#### 量子状態の初期化

```python
# initial state
qs = QState(qnum)
qs.hadamard(id_in)
```

量子状態を初期化して、入力レジスタ全体にアダマールをかけます。これで、3ビットの全状態の均等重ね合わせが実現できます。式で書くと、

```math
(\ket{000}+\ket{001}+\ket{010}+\ket{011}+\ket{100}+\ket{101}+\ket{110}+\ket{111}) \otimes \ket{0}
```

です。

#### 量子オラクルの実行

```python
# quantum oracle
qs.q_oracle(id_in,id_out,func)
```

論理関数funcに対する量子オラクルを実行します。今想定している真理値表を参照し、以下のようにしました。x=001,101のときのみTrueを返す関数です。各レジスタとこの関数を引数にしています。

```python
def func(x):

    if x == 1 or x == 5:
        y = 1
    else:
        y = 0
    return y
```

q_oracleの中身は以下です。

```python
def q_oracle(self,id_in,id_out,func):

    num = 2**len(id_in)
    for x in range(num):
        y = func(x)
        if y == 1:
            self.mpmct(id_in,id_out,x)

    return self
```

真理値yが1になったときのみ、そのときのxに対応したMPMCTを実行します。

mpmctの中身は以下です。

```python
def mpmct(self,id_in,id_out,x):

    bitnum = len(id_in)
    
    for i in range(bitnum):
        bit = (x>>i)%2
        if bit == 0:
            self.x(id_in[i])
            
    self.multi_cx(id_in,id_out)

    for i in range(bitnum):
        bit = (x>>i)%2
        if bit == 0:
            self.x(id_in[i])

    return self
```

Multi-Controlled NOT(multi_cx)をかけるのですが、その前後で、入力値xのビット列の中の負極性のビット番号にXゲートをかけます。

multi_cxの中身は以下です。

```python
def multi_cx(self,id_in,id_out):

    #
    # hadamard
    #
    self.h(id_out[0])

    #
    # controlled-P(psi), psi=pi/(2**(bitnum-1))
    #
    bitnum = len(id_in)
    psi = 1.0/(2**(bitnum-1)) # unit=pi(radian)
    gray_pre = 0
    for gray in gray_code(bitnum):
        if gray == 0:
            continue
        msb = len(str(bin(gray)))-3
        chb = len(str(bin(gray^gray_pre)))-3
        if gray != 1:
            if chb == msb:
                chb -= 1
            self.cx(id_in[chb],id_in[msb])
        self.cp(id_in[msb],id_out[0],phase=psi)
        psi = -psi
        gray_pre = gray
    
    #
    # hadamard
    #
    self.h(id_out[0])

    return self
```

先程説明したグレイコードを使った巧妙な手順を実行しています。アダマールを前後にかけるとして、その間でMulti_controlled Pを実行します。ちょっと複雑なので、言葉で説明するよりコードを見ていただいた方が良いと思います。中でグレイコードを1つずつ前に進めながらforループを回しています。

k番目のグレイコードは、以下のように、kをビットシフトしてkとのXORを実行することで生成できます。

```python
def gray_code(n):

    for k in range(2**n):
        yield k^(k>>1)
```

このq_oracleを実行すると、量子状態は、以下のような形になってるはずです。

```math
\ket{000}\otimes\ket{0}+\ket{001}\otimes\ket{1}+\ket{010}\otimes\ket{0}+\ket{011}\otimes\ket{0}+\ket{100}\otimes\ket{0}+\ket{101}\otimes\ket{1}+\ket{110}\otimes\ket{0}+\ket{111}\otimes\ket{0}
```

#### 測定と出力

最後に測定して、結果を得ます。

```python
# print results
in_list,out_list = qs.result(id_in,id_out)
for i in range(len(in_list)):
    print("{0:03b} -> {1:d}".format(in_list[i],out_list[i]))
```

resultの中で、全体の量子ビットを何度も測定して、制御ビット(入力)に相当するビット列と標的ビット(出力)に相当するビット列に分けたリストを返すようにしています(resultの関数定義を参照)。そして、そのリストを順に表示します。


### 結果

実行結果は、以下の通りです。

```
000 -> 0
010 -> 0
011 -> 0
100 -> 0
110 -> 0
111 -> 0
001 -> 1
101 -> 1
```

001または101の場合、1(True)で、それ以外の場合、0(False)が返る論理関数が、量子回路で実現できました。


## おわりに

あることを実現する量子回路は1つに決まるわけではありません。今回、このような形で論理関数を実現しましたが、実用を考えると、なるべく効率の良いもの(例えば、使用する量子ゲート数が少ないものとか、量子ビット数が少ないもの)が求められます。量子回路をいかに最適化するかというのは、各所で研究が進められている興味深いテーマの1つです。参考にさせていただいた[「量子回路と最適化」](http://www.orsj.or.jp/archive2/or63-6/or63_6_342.pdf)は、まさにその話題に関する解説記事です。これ以外にもいろいろな方向で取り組まれているようなので、今後のお勉強ネタの１つに「量子回路の最適化」も加えて行きたいと思います(やることがどんどん増えていく、、悩)。

## 訂正(2019.7.30)

Multi-Controlled NOTの説明で、Z軸周りの回転ゲート$Rz$と記載していた箇所は、すべて位相シフトゲート$P$の間違いでした。各々の定義は、

```math
Rz(\theta) = 
\begin{pmatrix}
e^{-i \theta /2} & 0 \\
0 & e^{i \theta /2}
\end{pmatrix},

P(\theta) =
\begin{pmatrix}
1 & 0 \\
0 & e^{i \theta}
\end{pmatrix}

```

です。違いは全体にかかっている位相項のみなので、今回の場合どっちで説明しても結果は同じでした。が、厳密には違います。

```math
H Rz(\theta) H = Rx(\theta)
```

という公式から$\theta=\pi$とすれば$HRz(\pi)X=HZH=X$だろ、と思い込んでしまっていましたが、$HRz(\pi)H = Rx(\pi) = iX$になるのでした。つまり、$Rz$を使って上の手順を踏むと、位相項iだけ違うMulti-Controlled NOTができあがります。$Rz$を$P$に変えると$HP(\pi)H=HZH=X$なので、位相項は現れません。というわけで、訂正します。


以上