Quanvolutional Neural Network

#はじめに
前回量子GANに関する記事を書きましたが、今回は、量子回路とConvolutional Neural Network(CNN)を組み合わせたようなアルゴリズムに関する論文を見つけたのと、ネットワークの構造などの情報も書いてあったので、MDR社のblueqatとpytorchでコードを書いてみることにしました。

論文のリンクはこちらです。
https://arxiv.org/pdf/1904.04767.pdf

#本論文の要旨
ざっくり書くと以下です。

1. Convolutional Neural NetworkならぬQuanvolutional Neural Network(QNN)を提案した。
2. QNNは、適当なqubit数の量子回路の測定値から作ったフィルタの処理N個とCNNから構成されたアルゴリズムである。
3. 本論文ではCNNとQNNと、Randomな非線形変換処理をCNNの前に加えたRandom CNNを作成し、MNISTの手書きの数字の認識において性能を比較した。結果、QNNがCNNより早く収束した。

古典的なNNでは次元の呪いといって、扱うデータの次元数が増えるとともに計算量が増大するという問題点がありますが、量子NNでは重ね合わせの量子状態で大規模な並列計算をすることにより、計算時間の指数関数的スピードアップが見込めるとのことです。
QNNでやられている処理は、CNNの畳み込み層と同じようなはたらきをしますが、ランダムな量子回路を使って、局所領域ごとに入力画像を変換するところがCNNと異なる点であり、これにより学習精度を上げられると論文では仮定されています。

#本論文でのCNNとQNNとRandom CNN
CNNとQNNを図で表すと以下のようになります。QNNはCNNの前にquanv1というquanvolutional filterの処理が追加されています。Random CNNは入力画像の各画素に乱数による非線形変換を行ったもので、quanv1の部分を書き換えて作成されていました。

(作図に使った情報)https://towardsdatascience.com/a-comprehensive-guide-to-convolutional-neural-networks-the-eli5-way-3bd2b1164a53

#Quanvolutional filter (quanv1)の処理
これも図にまとめてみました。量子回路の部分は、10種類のゲートをランダムに選んだqubitに作用させ、ゲートの順番もランダムに配置させることで作成されました。

#blueqatとpytorchによるQNNの実装
以上で述べたようなCNN、QNN、Random CNNを模擬した実装を行いました。「模擬」としたのは、以下の理由からです。

• 全ての組み合わせのqubitではなく限定された組み合わせのqubitにCNot~ControlledUgateをあてはめた
• ゲートを作用させる順番や数を決め打ちにした(あとで述べますが、計算時間がすごかったので。。。)

コード書くにあたり、下記のサイトを参考にしました。

• 古典CNNの書き方 --> https://qiita.com/fukuit/items/215ef75113d97560e599
• QNNの量子回路の書き方 --> https://github.com/Blueqat/Blueqat-tutorials/blob/master/tutorial-ja/014_grover_ja.ipynb
• QNNの量子回路の結果を使った3x3の畳み込み処理 --> https://masamunetogetoge.com/deeplearning2

次から順番にコードを紹介します。

#コード
###古典CNN
フィルタサイズや全結合層の部分を論文の手法と同じにしました。

net.py

import torch.nn as nn
import torch.nn.functional as F

class Net(nn.Module):
    def __init__(self):
        super(Net, self).__init__()
        self.conv1 = nn.Conv2d(1, 50, 5)
        self.pool1 = nn.MaxPool2d(2, 2)
        self.conv2 = nn.Conv2d(50, 64, 5)
        self.pool2 = nn.MaxPool2d(2, 2)
        self.fc1 = nn.Linear(1024, 128)
        self.dropout2 = nn.Dropout2d()
        self.fc2 = nn.Linear(128, 10)

    def forward(self, x):
        x = F.relu(self.conv1(x))
        x = self.pool1(x)
        x = F.relu(self.conv2(x))
        x = self.pool2(x)
        x = x.view(-1, 1024)
        x = F.relu(self.fc1(x))
        x = self.dropout2(x)
        x = self.fc2(x)

        return x

###QNNの量子回路
制御Uゲートなど2qubit必要なものは、0~8のqubitに対し、どのqubitをつなげるかを0~8からランダムに選びました。元のqubitをl=[0,...8],つなげる先をl2=[0,...8]とリストにしてl2の中身をランダムに入れ替え、それぞれのゲートにl[0]とl2[0]、l[1]とl2[1]というようにqubitを割り当てました。このときl[0]=l2[0]=0のようなペアは除外しました。Xゲートなど1qubit必要なものはl[0]、l[1]などと順番に割り当てました。

quanv.py

import io
import sys
from blueqat import Circuit
import math
import random
import numpy as np

def QNet(col):
    l = [0, 1, 2, 3, 4, 5, 6, 7, 8]
    l2 =[0, 1, 2, 3, 4, 5, 6, 7, 8]
    while ((l[0] is l2[0]) or (l[1] is l2[1]) or (l[2] is l2[2]) or (l[3] is l2[3])
           or (l[4] is l2[4]) or (l[5] is l2[5])
           or (l[6] is l2[6]) or (l[7] is l2[7]) or (l[8] is l2[8])):
        random.shuffle(l)

    #######   ここからが9C2の組み合わせ回路用のゲート #################
    """
    #Control Ugate
    """

    id = np.random.randint(0, 15, size=15)

    id[0]=l[0]
    id[1] = l2[0]
    rn = np.random.rand()
    if (col[id[0]] > 0.):
        rcu = Circuit().i[8].x[id[0]].h[id[1]].cu1(rn * math.pi)[id[0], id[1]].m[:]
    else:
        rcu = Circuit().i[8].h[id[1]].cu1(rn * math.pi)[id[0], id[1]].m[:]

    """
    #CNOT(CX) gate
    """
    id[2] = l[1]
    id[3] = l2[1]
    if (col[id[2]]>0.):
        rcnot = Circuit().i[8].x[id[2]].cx[id[2], id[3]].m[:]   #control bit=1
    else:
        rcnot = Circuit().i[8].cx[id[2], id[3]].m[:]    #control bit=0

    """
    #Swap
    """
    id[4] = l[2]
    id[5] = l2[2]
    if (col[id[4]]>0.):
        rswap = Circuit().i[8].swap[id[4], id[5]].m[:]
    else:
        rswap = Circuit().i[8].x[id[4]].swap[id[4], id[5]].m[:]

    """
    ##sqrt swap
    """
    id[6] = l[3]
    id[7] = l2[3]
    if (col[id[6]] > 0.):
        sqrtswp = Circuit().i[8].x[id[6]].cx[id[6], id[7]].t[id[7]].h[id[6]].t[id[6]].cx[id[7], id[6]].tdg[id[6]].cx[id[7], id[6]].h[id[6]].cx[id[6], id[7]].m[:]
    else:
        sqrtswp = Circuit().i[8].cx[id[6], id[7]].t[id[7]].h[id[6]].t[id[6]].cx[id[7], id[6]].tdg[id[6]].cx[id[7], id[6]].h[id[6]].cx[id[6], id[7]].m[:]

    """
    ##RX(theta)
    """
    id[8] = l[4]
    rn = np.random.rand()

    if (col[id[8]] > 0.):
        xgt=Circuit().i[8].x[id[8]].rx(2*rn*math.pi)[id[8]].m[:]
    else:
        xgt = Circuit().i[8].rx(2*rn * math.pi)[id[8]].m[:]

    """
    ##RY(theta)
    """
    id[9] = l[5]
    rn = np.random.rand()
    if (col[id[9]] > 0.):
        ygt = Circuit().i[8].x[id[8]].ry(2*rn * math.pi)[id[9]].m[:]
    else:
        ygt = Circuit().i[8].ry(2*rn * math.pi)[id[9]].m[:]

    """
    ##RZ(theta)
    """
    id[10] = l[6]
    rn = np.random.rand()
    if (col[id[10]] > 0.):
        zgt = Circuit().i[8].x[id[10]].rz(2*rn * math.pi)[id[10]].m[:]
    else:
        zgt = Circuit().i[8].rz(2*rn * math.pi)[id[10]].m[:]

    """
    ##U(theta)
    """
    id[11] = l[7]
    rn = np.random.rand()
    if (col[id[11]] > 0.):
        U = Circuit().i[8].x[id[11]].h[id[11]].u1(2*rn * math.pi)[id[11]].m[:]
    else:
        U = Circuit().i[8].h[id[11]].u1(2*rn * math.pi)[id[11]].m[:]

    """
    ##T
    """
    id[12] = l[8]
    rn = np.random.rand()
    if (col[id[12]] > 0.):
        T = Circuit().i[8].x[id[12]].rz(math.pi/4)[id[12]].m[:]
    else:
        T = Circuit().i[8].rz(math.pi/4)[id[12]].m[:]

    """
    ##P   P=T^2
    """
    id[13] = l[0]
    rn = np.random.rand()
    if (col[id[13]] > 0.):
        P = Circuit().i[8].x[id[13]].rz(math.pi/2)[id[13]].m[:]
    else:
        P = Circuit().i[8].rz(math.pi/2)[id[13]].m[:]
    #aaa = (Circuit() + rcu + rcnot + rswap + sqrtswp + xgt + ygt + zgt + U + T + P).run(shots=1)
    aaa = (Circuit() + rcu).run(shots=1)  #上で定義した量子回路を組み合わせて実行し、観測結果を格納

    #######   ここまでが9C2の組み合わせ回路用のゲート #################

    ##観測結果を取り出す
    result = []

    with io.StringIO() as f:
        sys.stdout = f
        print(aaa)
        text = f.getvalue()
        sys.stdout = sys.__stdout__
        a = [float(text[10]),float(text[11]),float(text[12]),float(text[13]),
             float(text[14]),float(text[15]),float(text[16]),float(text[17]),float(text[18])]

        b=float(text[10]) + float(text[11]) + float(text[12]) + float(text[13])
        +float(text[14]) + float(text[15]) + float(text[16]) + float(text[17]) + float(text[18])

        if (b !=0):
            a = [float(text[10])/b, float(text[11])/b, float(text[12])/b, float(text[13])/b,
             float(text[14])/b, float(text[15])/b, float(text[16])/b, float(text[17])/b, float(text[18])/b]

        result.append(a)
    return result

###QNNの量子回路の結果を使った3x3の畳み込み
今回一番苦戦したのは上記の量子回路の測定結果を使って畳み込み処理をする部分でした。入力画像1枚の中で、座標によって異なる3x3のフィルタをかけるため、pytorchのnn.Conv2Dが使えなかったからです。
そこでim2colという畳み込みの手法を少し変えて使うことにしました。

im2colとはimage to columnの略で、画像をある種の行列に変換する処理です。この処理とほかの行列演算(内積含む)を組み合わせることで、普通の畳み込みで発生するfor loopの数を減らしながら計算を高速化します。
下図に例を示しますが、今回のim2colを使った畳み込みでは、入力画像を並べ替えた行列に、3x3の畳み込みフィルタを9行1列にしたものを内積しています。

このことと量子回路の測定結果を用いて以下の手順で畳み込みを実装しました。

1. im2colで作った行列と同じサイズの行列を作り、その中に量子回路の測定値を格納した。
2. 入力画像と量子回路由来の行列の、同じ要素の積をとった。
3. 上記2の結果を足し合わせるため、1を内積した。

※ちなみにrandom CNNは、1の部分を量子回路の測定値の代わりにrandomな値を格納して作成しました。

####メインのコードを次に示します。
学習画像を読み込んで、CNNやrandom CNN, QNNの量子回路の処理およびその畳み込み処理を呼び出す部分などが入っています。random CNNかQNNを実行するのにcnn_flagを立てました。

mnist.py

import argparse
import time
import matplotlib.pyplot as plt
import numpy as np
import torch
import torch.nn as nn
from torch.utils.data import DataLoader
import torchvision.transforms as transforms
from torchvision.datasets import MNIST
import torch.optim as optim
from net import Net
import quanv


def parser():
    '''
    argument
    '''
    parser = argparse.ArgumentParser(description='PyTorch MNIST')
    parser.add_argument('--epochs', '-e', type=int, default=1,
                        help='number of epochs to train (default: 2)')
    parser.add_argument('--lr', '-l', type=float, default=0.01,
                        help='learning rate (default: 0.01)')
    args = parser.parse_args()
    return args

def main():
    '''
    main
    '''
    args = parser()
    c=0

    #DatasetをDataLoader に渡してやるとシャッフルなどしてイテレーションしてくれるが、
    # MNISTのデータはPIL Imageで入っているのでそのまま渡すと
    # TypeError: batch must contain tensors, numbers, dicts or lists; found <class 'PIL.Image.Image'>
    # のエラーになる。 そこでtransforms.ToTesnor()でTensorに変換してやる必要がある。
    transform = transforms.Compose(
        [transforms.ToTensor(),
         transforms.Normalize((0.5, ), (0.5, ))])

    trainset = MNIST(root='./data',
                     train=True,
                     download=True,
                     transform=transform)
    testset = MNIST(root='./data',
                    train=False,
                    download=True,
                    transform=transform)

    trainloader = DataLoader(trainset,
                             batch_size=100,
                             shuffle=True,
                             num_workers=2)
    testloader = DataLoader(testset,
                            batch_size=100,
                            shuffle=False,
                            num_workers=2)

    for x, t in trainloader:
        break

    def im2col(input_data, filter_h, filter_w, stride=1, pad=1):
        """
        input_data : (データ数, チャンネル, 高さ, 幅)の4次元配列からなる入力データ
        filter_h : フィルターの高さ
        filter_w : フィルターの幅
        stride : ストライド
        pad : パディング

        Returns
        col : 2次元配列
        """
        N, C, H, W = input_data.shape  # データの形状
        out_h = (H + 2 * pad - filter_h) // stride + 1  # 出力の縦
        out_w = (W + 2 * pad - filter_w) // stride + 1  # 出力の横

        img = np.pad(input_data, [(0, 0), (0, 0), (pad, pad), (pad, pad)], 'constant')  # 必要ならパディング
        col = np.zeros((N, C, filter_h, filter_w, out_h, out_w))  # 答えが入ってくる行列の準備, 出力の大きさも入れておく

        for y in range(filter_h):
            y_max = y + stride * out_h
            for x in range(filter_w):
                x_max = x + stride * out_w
                col[:, :, y, x, :, :] = img[:, :, y:y_max:stride, x:x_max:stride]  # ストライド毎にデータを区切る

        col = col.transpose(0, 4, 5, 1, 2, 3).reshape(N * out_h * out_w, -1)  # 出力の大きさ*(入力データの数)の行列に変換
        return col

    def forward(x,W, b,cnn_flag, stride=1, pad=1):  # x:入力画像
            FN, C, FH, FW = W.shape
            N, C, H, W2 = x.shape
            out_h = 1 + int((H + 2 * pad - FH) / stride)  # 出力縦
            out_w = 1 + int((W2 + 2 * pad - FW) / stride)  # 出力横

            col = im2col(x, FH, FW, stride, pad)

            ##### ランダムCNNのとき。
            if cnn_flag==0:
                col_pre=np.random.random((N * out_h * out_w, 9))
                col_sum=np.sum(col_pre,axis=1)
                for i in range(N * out_h * out_w):
                    col_rep=col_pre[i]/col_sum[i]
                col_W = col_rep
                out = np.dot(col, col_W) + b

            #### 量子NNのとき
            if cnn_flag==1:

                net3=[]
                for i in range (N * out_h * out_w):
                    net2 = quanv.QNet(col[i])[0]
                    net3.append(net2)
                col_W2=np.array(net3)

                out = np.dot((col* col_W2),np.ones(9))

            out = out.reshape(N, out_h, out_w, -1).transpose(0, 3, 1, 2)  # (N,FN, OH,OW ) に整形
            return out

    ##### 学習に使われる画像のチェック
    # originalの学習データを出力
    inp=x[0].reshape(1,1,28,28)
    plt.imshow(inp.reshape(28,28),cmap="gray")
    plt.title("original")
    plt.savefig('out_original.png', bbox_inches='tight')
    W = np.array([1,1,1,1,1,1,1,1,1]).reshape(1, 1, 3, 3)  #適当な3x3のフィルタを作っておく
    b=0

    # QNN,random CNNの一層目で処理された学習データを出力
    #cnn_flag = 1
    #c=forward(x, W, b, cnn_flag,stride=1, pad=1)
    #plt.imshow(np.sum(c[0], axis=0), cmap="gray")
    #plt.title("qnn")
    #plt.savefig('out_qnn.png', bbox_inches='tight')

    #cnn_flag =0
    #c=forward(x, W, b, cnn_flag,stride=1, pad=1)
    #plt.imshow(np.sum(c[0], axis=0), cmap="gray")
    #plt.title("random_cnn")
    #plt.savefig('out_random_cnn.png', bbox_inches='tight')
    ###############################

    # model
    net = Net()

    # define loss function and optimier
    criterion = nn.CrossEntropyLoss()
    optimizer = optim.SGD(net.parameters(),
                          lr=args.lr, momentum=0.99, nesterov=True)

    # train
    for epoch in range(args.epochs):
        running_loss = 0.0
        num=0
        for i, (inputs, labels) in enumerate(trainloader, 0):
            # zero the parameter gradients
            optimizer.zero_grad()

            # forward + backward + optimize
            W = np.array([1/9, 1/9, 1/9, 1/9, 1/9, 1/9, 1/9, 1/9, 1/9]).reshape(1, 1, 3, 3)
            b = 0
            cnn_flag=1  #0ならrandom CNN, 1ならQNNを実行

            ## 古典CNN以外の入力データ処理。 ######
            num+=1
            print("\ri={0}".format(num),end="")
            #if (i%5==0):          #すべてのインプット画像に処理を行うと時間がかかるので、一部の画像にかけて計算を省略したい場合は入れる。
            inputs=forward(inputs,W, b, cnn_flag,stride=1, pad=1)
            inputs = torch.from_numpy(inputs).clone()
            ###################################

            outputs = net(inputs.float())

            loss = criterion(outputs, labels)
            loss.backward()
            optimizer.step()

            # print statistics
            running_loss += loss.item()
            if i % 100 == 99:
                print('[{:d}, {:5d}] loss: {:.3f}'
                      .format(epoch+1, i+1, running_loss/100))
                running_loss = 0.0
    print('Finished Training')

    # test
    correct = 0
    total = 0
    with torch.no_grad():
        for (images, labels) in testloader:
            outputs = net(images)
            _, predicted = torch.max(outputs.data, 1)
            total += labels.size(0)
            correct += (predicted == labels).sum().item()

    print('Accuracy: {:.2f} %'.format(100 * float(correct/total)))


if __name__ == '__main__':
    start_time = time.time()
    main()
    print('elapsed time: {:.3f} [sec]'.format(time.time() - start_time))

#結果
今回quanv.pyで10種類量子ゲートを作りました。CNN、QNN、Random CNNの実行結果を比較しようとしましたが、QNNをまともに実行しようとするとかなりの時間がかかってしまったので、

1. CNN
2. 制御UゲートからなるQNNを学習データすべてにかけた場合 (QNN)
3. 全学習データのうち20%に制御Uゲート,CNOTゲート, SqrtSwapゲートからなるQNNをかけた場合 (CNN+QNN)
4. random CNN

の結果を比較しました。60000枚の学習データを100枚ずつ学習させてepoch数を1としました。

####学習時のlossの変化を見てみます。
実装したコード内では、CNNが一番早くlossの値が落ちていました。最終的なlossは
random CNN < CNN < CNN+QNN < QNNとなりました。

論文で報告された結果と比較(※)すると、Random CNNに関しては同じような収束の仕方になっていました。CNNは論文の方が収束が遅かったです。

※iterationの軸が同じだと仮定しています。論文では60000枚の学習データを何枚ずつ学習させたかの情報がなかったので、自分の実装と一対一対応ではないかもしれません。また、論文ではquanv1の処理, randomCNNの変換を25回かけたとあるので、そこでも結果にずれが発生していると思います。

####検証時の精度を見てみます。
10000枚の検証用データを使った精度は以下の表のようになりました。Random CNNが一番高く、CNN+QNNがその次に高かったです。計算時間も載せましたが、QNNが一番計算時間が長かったです。。

(余談ですが、制御Uゲートを使ったQNNの計算時間が長すぎたので、3種のゲートを使ったQNNはあきらめてCNN+QNNをやりました。)

#考察
QNNではquanvolutional filter, random CNNではこれをrandomな値で置き換えたフィルタをCNN処理の前に追加しました。これらのフィルタ係数は、学習時に誤差逆伝播に従って更新されるものではありません。なので、結局は学習をうまく進めるために学習データに何らかの処理を行うことに相当すると考えられます。CNNにインプットする直前のCNN(要はもとの学習データ)、QNN、CNN+QNN、Random CNNの学習データを見てみると次のようになりました。

QNNでは制御Uゲートしか使っていないためか、画像が粗くなっていました。CNN+QNNでは、3種類のゲートを使ったためか、画像がQNNよりスムーズになっていました。どちらも白い数字の部分が薄くなってしまっていたため、もっとゲートを増やしたり、ランダム化したら、うまいこと元の学習データやrandom CNNの画像に近づくのかもしれません。

#まとめ
QNNを実装できました。ですが、時間がかかりすぎたため、もっと短縮できる方法がないかが今後の課題になりそうです。