Go to Qiita Advent Calendar 2024 Top
LoginSignup
97
74

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

量子コンピュータAdvent Calendar 2017

Day 23
@openql(OpenQL Project)

SymPy で量子プログラミングを体験してみましょう

Last updated at Posted at 2017-12-23

はじめまして、こんにちは。OpenQLプロジェクト1の中のひと、です。
本記事の対象は、次のような方を想定しております。

  • 量子コンピューターや量子プログラミング2に興味がある方
  • 基礎的なPythonの知識がある方
  • 数式に極度のアレルギーのない方
  • 量子コンピューターの仕組みや計算方法をなんとなく知っている方
    （本記事の前に、@tsujishin さんの「量子情報科学序論 IBM Qを動かして学ぶ量子コンピュータ」を事前に読んで頂くと、より理解が深まると思います）

それでは、始めましょう。よろしくお願いします。

数学・科学計算には欠かせない！？「Mathematica」

皆さんは、Mathematica というソフトウェアをご存知でしょうか。
例えば、何か計算式をプログラムで扱おうとしたときに、一般的なプログラムが数値計算を得意としているのに対して、$x$が含まれる関数式のまま扱うことができるアプリケーションです。微分や積分ができたり、因数分解ができたり、グレブナー基底まで求めてくれたりする、理系関係者には涙が出るほど嬉しいソフトウェアです。最近では、RaspberryPi の標準 OS には（無償で）搭載されているというので、RaspberryPiは、安価な Mathematica マシンとして、一部の理系人には注目されたこともあります。
Mathematica は Wolfram さんが中心に開発されており、最近では、この Mathematica で使われているコマンド群が進化して、Wolfram言語 と呼ばれ、プログラミング言語の仲間入りもしています。

さて、そんな（熱狂的な一部の人には）著名な Mathematica ですが、有償であるため、利用者は研究開発の分野に限られていました。最近になって、オンラインでも使えるような環境3が提供されたりしていますが、それでも、購入前の評価目的で公開されている様です。基本的にはWebブラウザからの利用になり、プログラムとして動作させるには難しく、使い方は限定的となります。

数式を数式として扱える「SymPy」

理系人にとって涙が出るほど嬉しい Mathematica を無償で使えないものだろうか？そう考えるのが、OSS界隈のエンジニアの素晴らしいところです。そして登場したのが「SymPy」です。
**Mathematica のように、数式を数式のまま扱える処理系をめざして開発されたのが「SymPy」です。**SymPy は、完全に Python 上で動作し、非常に軽量で、かつ依存するソフトウェアが少ないライブラリとして「無償で」提供されています。
SymPy は、Mathematica と同じように、通常の Python で使われる「プログラム上の変数」を扱うのではなく、まさに数学や理科の教科書にでてくるような$x$とか$y$という「数学上の変数」をそのまま扱います。**変数を抽象的なシンボルのまま、代入や計算が行えるのがSymPyの特徴です。**この「Symbol」＋「Python」がこのライブラリの名前の由来なのでしょう。

Qiita にも SymPy の利用者は多く、カテゴリで検索すると、面白い使い方を投稿されている方がいらっしゃいます。

  • @hiroyuki827 さんが、英語のチュートリアルを日本語4にされています。
  • @mrrclb48z さんが、SymPy を使った色々な利用例を投稿されています。（中学校数学を sympy を使って解く5は面白い用例だと思います。）

本記事では、SymPyの数学系の使い方には触れません。ある程度の基本的な使い方は想定できていることを前提として進めます。ほかにも、Qiita には、SymPy の利用方法について投稿がありますので、基本的な使い方は、それらの記事を参考にしてください。

「SymPy」で量子プログラミング

さて、ここからが本題です。
この SymPy パッケージには、数学計算に便利な基礎的な利用用途とは別に、様々な応用計算が標準で含まれています。実は、その一つに、本記事の主題である「量子コンピューターの計算」ができます。

それでは「SymPyで量子プログラミング」を体験してみましょう。
まずは、SymPyをインストールします。Pythonの環境さえ整っていれば、pip でインストール6できます。

sympyのインストール
$ pip install sympy

量子計算に必要なモジュールは、SymPy の中でも極一部です。物理計算のなかの量子計算である sympy.physics.quantum にあります。

sympy.physics.quantumのインポート
>>> from sympy.physics.quantum import *
         # ↑基本機能のインポート
>>> from sympy.physics.quantum.qubit import Qubit,QubitBra
         # ↑ブラケットを使えるようにします
>>> from sympy.physics.quantum.gate import X,Y,Z,H,T,S,CNOT,SWAP
         # ↑基本的な量子ゲートを使えるようにします
>>> from sympy.physics.quantum.gate import IdentityGate as I
         # ↑恒等変換は省略されないため I と省略できるようにします

量子ゲートや、恒等変換の省略形は、単純なプログラムではそのまま使えるかもしれませんが、少し大きなプログラムでは、XGate, YGate, ZGate, HadamardGate, TGate, PhaseGate, CNotGate, SwapGate, IdentityGateを使いましょう。（本記事では、以降も省略形のゲートを使って説明します）

量子ビットの表現

SymPy で量子ビットは、ケット・ベクトルで表現し、計算するまではケット・ベクトルのままで処理されます。ケット・ベクトルは、Qubitクラスを使います。Qubitオブジェクトを生成するときに指定する引数は、

  • 各量子ビットをカンマ区切りのリストで指定する
  • 量子ビット列を文字列で指定する

の２通りの指定の仕方があります。いずれの指定も、0,1 以外を指定するとエラーになります。

Qubitの生成
>>> from sympy.physics.quantum.qubit import Qubit
>>> Qubit(1,0,0)
|100>
>>> Qubit('0101') # 4量子ビット. Qubit('101') は3量子ビットで意味が異なります.
|0101>
>>> q = Qubit('0'*100)    # 100桁の0並びでも処理できます. (*)
>>> q                     # 100量子ビットのケットベクトルです.
|000...00>

(*) SymPyでは、 シンボルのままで処理が進むため、例えば、q=Qubit('0'*100000) という式もリソースがあれば動作します。しかし、その領域の確保も同時に行うため、非常に時間がかかります。（ただし、実際に確保される領域は、本来の量子ビットが表せる$2^{100,000}$個の複素数ではなく、単に100,000文字列の領域でしかありません。）

Qubitクラスには、次のようなメソッドが提供されています。

Qubitが提供するメソッド
>>> Qubit('0101')
|0101>
>>> q.nqubits
4
>>> len(q)
4
>>> q.dimension
4
>>> q.qubit_values
(0, 1, 0, 1)
>>> q.flip(1)  # 引数で指定された位置(*)の右から2番目の量子ビットを反転ます.
|0111>

(*) 「Qubitのビット位置の指定」は、一番右の量子ビットが 0 で、左へ向かって順番に +1 ずつした数字で示します。$n$量子ビットをケット・ベクトルにすると、$ \lvert q_n q_{n-1} ... q_2 q_1 q_0 \rangle $ と表せますが、この添え字の番号がビット位置を示します。

複素共役を表すための dagger（$\dagger$） をQubitオプジェクトに作用させると、ブラ・ベクトルになります。内積も計算できます。シンボルのままで内積にするには、単にブラ・ベクトルQubitBraとケット・ベクトルQubitを掛ける（*）だけです。実際の値を計算するには、この内積オブジェクトのdoit()メソッドで計算します。

複素共役を作用してブラ・ベクトルに
>>> from sympy.physics.quantum.dagger import Dagger
>>> Dagger(q)
<0101|
>>> type(Dagger(q))
<class 'sympy.physics.quantum.qubit.QubitBra'>
>>> ip = Dagger(q)*q
>>> ip
<0101|0101>
>>> ip.doit()
1

ユニタリーの計算

次に、ユニタリー演算のための、ユニタリー行列をみてみましょう。
SymPy には、次の基本ゲートがあらかじめ定義されています。

  • １量子操作ゲート
    • Idゲート（IdentityGate）
    • パウリゲート（XGate,YGate,ZGate）
    • アダマールゲート（HadamardGate）
    • Ｔゲート（TGate）
    • 位相ゲート（PhaseGate）
  • ２量子操作ゲート
    • CNOTゲート（CNotGate）
    • SWAPゲート（SwapGate）

１量子操作ゲートは、引数を１つとり、２量子操作ゲートは、引数を２つとります。引数は、どの位置の量子ビットに操作するかを示す数字です。上の「Qubitのビット位置の指定」で説明した数字が引数になります。

量子ゲートの次元は計算時に確定
また、SymPy では、量子ゲートもシンボリックに扱われるため、ゲートを行列として考える際の次元に関する情報は、実際に計算するとき（後づけで）付加されます。
例えば、パウリＸであるXGateで説明します。通常１量子ビットに作用させることを考えると、XGate(0)は、$$
\displaystyle XGate(0) = \left( \begin{array}{cc} 0 & 1 \cr 1 & 0 \end{array}\right) $$ と、$2 \times 2 $行列です。同じXGate(0)でも、２量子ビットに作用させるときには、$$
\displaystyle XGate(0) = I \otimes X = \left( \begin{array}{cc} 0 & I \cr I & 0 \end{array}\right)
= \left( \begin{array}{cccc} 0 & 0 & 1 & 0 \cr 0 & 0 & 0 & 1 \cr 1 & 0 & 0 & 0 \cr 0 & 1 & 0 & 0 \end{array}\right)$$ と、$4 \times 4 $行列となります。さらに、XGate(1)は、$4 \times 4 $行列以上の次元をもち、例えば、２量子ビットに作用させるときには、$$
\displaystyle XGate(1) = X \otimes I = \left( \begin{array}{cc} X & 0 \cr 0 & X \end{array}\right)
= \left( \begin{array}{cccc} 0 & 1 & 0 & 0 \cr 1 & 0 & 0 & 0 \cr 0 & 0 & 0 & 1 \cr 0 & 0 & 1 & 0 \end{array}\right)$$ となります。同様に３量子ビットに作用させる場合は（行列表記は省きますが）、次のようになります。$$ XGate(0) = I \otimes I \otimes X \\
XGate(1) = I \otimes X \otimes I \\
XGate(2) = X \otimes I \otimes I \\ $$ 少し踏み込んで、XGate同士の積についても例を挙げておきましょう。$$
\displaystyle XGate(0)*XGate(1) = X \otimes X = \left( \begin{array}{cc} 0 & X \cr X & 0 \end{array}\right)
= \left( \begin{array}{cccc} 0 & 0 & 0 & 1 \cr 0 & 0 & 1 & 0 \cr 0 & 1 & 0 & 0 \cr 1 & 0 & 0 & 0 \end{array}\right)$$ となります。なお、引数が同じXGate同士の積は、恒等変換 $I$ になります。（ $ XGate(0) * XGate(0) = I $、$ XGate(1) * XGate(1) = I $ ...）

実際の計算は「qapply」
SymPy を使ったプログラムでは、計算がシンボルで進みます。そのため、 SymPy では、ユニタリーを量子ビットに作用するには、明示的に「qapply」メソッドを行って実際の行列計算を行います。次の例をみてください。qapplyを行うと、計算できるところまで計算されます。SymPy の特性を生かして、$ \sqrt{2} $ は、qapplyしても近似的な数値にならず、sqrt(2)のままで扱われます。

qapplyの用例
>>> from sympy import sqrt
>>> from sympy.physics.quantum.qubit import Qubit
>>> from sympy.physics.quantum.gate import HadamardGate
>>> from sympy.physics.quantum.qapply import qapply
>>> qapply(HadamardGate(0)*Qubit('1'))
sqrt(2)*|0>/2 - sqrt(2)*|1>/2
>>> # Hadamard on bell state, applied on 2 qubits.
>>> psi = 1/sqrt(2)*(Qubit('00')+Qubit('11'))
>>> qapply(HadamardGate(0)*HadamardGate(1)*psi)
sqrt(2)*|00>/2 + sqrt(2)*|11>/2

>>> from sympy.physics.quantum.gate import CNOT
>>> from sympy.physics.quantum.qapply import qapply
>>> from sympy.physics.quantum.qubit import Qubit
>>> c = CNOT(1,0)
>>> qapply(c*Qubit('10')) # note that qubits are indexed from right to left
|11>

観測

ここまでの量子計算では、量子状態のまま計算されます。シミュレーション上の理論的な計算では、qapplyまで行えば、目的が達成することが多いです。 SymPy には、観測のためのメソッドも備わっています。measure_all()では、引数で渡された量子ビット（ケット・ベクトル）のすべての状態の測定確率を計算で求めて出力します。そのほかにも幾つか観測のためのメソッドが提供されます。例えば、measure_partial()は、第一引数で指定された量子ビットに対して、第二引数で指定された部分的な振幅の測定を実行します。

観測(measure_all,measure_partial)の用例
>>> from sympy.physics.quantum.qubit import Qubit, measure_all
>>> from sympy.physics.quantum.gate import H, X, Y, Z
>>> from sympy.physics.quantum.qapply import qapply
>>> c = H(0)*H(1)*Qubit('00')
>>> c
H(0)*H(1)*|00>
>>> q = qapply(c)
>>> measure_all(q)
[(|00>, 1/4), (|01>, 1/4), (|10>, 1/4), (|11>, 1/4)]

>>> from sympy.physics.quantum.qubit import Qubit, measure_partial
>>> from sympy.physics.quantum.gate import H, X, Y, Z
>>> from sympy.physics.quantum.qapply import qapply
>>> c = H(0)*H(1)*Qubit('00')
>>> c
H(0)*H(1)*|00>
>>> q = qapply(c)
>>> measure_partial(qapply(q),(0,))
[(sqrt(2)*|00>/2 + sqrt(2)*|10>/2, 1/2), (sqrt(2)*|01>/2 + sqrt(2)*|11>/2, 1/2)]
>>> measure_partial(qapply(q),(1,))
[(sqrt(2)*|00>/2 + sqrt(2)*|01>/2, 1/2), (sqrt(2)*|10>/2 + sqrt(2)*|11>/2, 1/2)]
>>> measure_partial(qapply(q),(0,1))  # measure_all(q)と同じ
[(|00>, 1/4), (|01>, 1/4), (|10>, 1/4), (|11>, 1/4)]

量子アルゴリズムを試す

ここまで、基礎的な使い方を解説してきました。ここからは、SymPy を使った量子プログラミングを見ていきましょう。基本的には、最初に必要な量子ビットを準備して、その量子ビットに、量子ゲートを操作していくのが量子プログラムです。最後は、観測まで行うべきですが、ユニタリー操作の論理的な結果だけを求めるのであれば、qapplyで計算するところまでをプログラムすれば十分です。シミュレーションとして、理論的な計算が必要なケースが多く、ここで取り上げる量子アルゴリズムの例もqapplyまでに留めます。

Grover のアルゴリズム

Grover のアルゴリズムは、既に sympy.physics.quantum.grover にまとめられて、検証できるようになっています。しかし、ここでは、このモジュールは使わずに、「IBMQ」で作った量子回路の例を SymPy で書いてみましょう。
２量子ビットの重ね合わせの状態 $ \frac{1}{2}(\lvert00\rangle+\lvert01\rangle+\lvert10\rangle+\lvert11\rangle) $ から、$\lvert11\rangle$を探索する量子回路を考えます。
grover_search_3_on2qubit_2.png
この回路を SymPy で計算してみます。

Groverで|11>を探索する例
>>> from sympy.physics.quantum import *
>>> from sympy.physics.quantum.gate import X,Y,Z,H,T,S,CNOT
>>> from sympy.physics.quantum.gate import IdentityGate as I
>>> from sympy.physics.quantum.qubit import Qubit
>>> qi = H(0)*H(1)*Qubit('00')
>>> marker_3 = H(0)*CNOT(1,0)*H(0)
>>> amplifier = H(1)*CNOT(1,0)*(Z(1)*X(0))*H(1)
>>> grover = amplifier * marker_3
>>> qapply( grover * qi)
|11>

Shor のアルゴリズム

Shor のアルゴリズムも、既に sympy.physics.quantum.shor にまとまっています。ここでは、Shor のアルゴリズムの詳細は本題ではないので、詳しくは述べません。次のような位相発見のためのメソッドが定義されていることだけご紹介しておきます。

位相発見メソッド
def period_find(a, N):
    """Finds the period of a in modulo N arithmetic

    This is quantum part of Shor's algorithm.It takes two registers,
    puts first in superposition of states with Hadamards so: ``|k>|0>``
    with k being all possible choices. It then does a controlled mod and
    a QFT to determine the order of a.
    """

重ね合わせ状態を作って、処理時間を計測

さて、最後の用例です。$n$次のアダマール操作で、すべての量子状態の重ね合わせ状態を作る操作を考えます。果たして、SymPy はどれほどの次元まで扱えるのでしょうか。さらに、その計算にはどれくらいの時間がかかるのですか。実際に調べてみます。

アダマールで重ね合わせ状態を作る例
#-*- using:utf-8 -*-
import time
from sympy.physics.quantum import *
from sympy.physics.quantum.qubit import Qubit
from sympy.physics.quantum.gate import HadamardGate

if __name__ == '__main__':
  for n in range(1,19): # このまま実行すると非常に時間がかかります.
                        # 19 を 5 程度にしてお試しください.
    p = '0' * n
    q = Qubit(p)               # n-qubit( |0...0> )を準備します.
    h = HadamardGate(0) 
    for i in range(1,n):
      h = HadamardGate(i) * h  # H(n)*H(n-1)*...*H(0) の操作
    start = time.time()        # 計算時間計測をスタート
    r = qapply(h * q)          # (H(n)*H(n-1)*...*H(0))|0...0> を計算
    elapsed_time = time.time() - start   # 計算時間計測をストップ
    print ("{0},{1}".format(n,elapsed_time))

このプログラムを手元のノートPC7で試してみました。手元の限られたリソースでも、19までは処理が行われました。しかし、$n$が大きくなると、重ね合わせ状態を生成するだけなのに、かなりの時間（一晩放置しておくレベルの時間）がかかりました。量子ビット数と処理時間（実測値）の関係をグラフにしました。縦軸が処理時間の秒数、横軸が量子ビット数です。（$n=19$は一晩放っておいた結果で、4時間15分ほどかかって計算が終わっていたようです。）
hadamard19.png
ただし、19量子ビットの計算でも、時間はかかったものの、リソース不足などのエラーは起きませんでした。つまり、$ 2^{19} $ パターンの重ね合わせた状態を生成できたことになります。それぞれの$n$量子ビットに対するユニタリー行列の演算（操作）は試していませんが、貧弱なリソース環境でもある程度の量子ビット数は扱えそうです。

これからの「SymPyで量子プログラミング」

ここまで見てきたように、SymPy は量子計算の有用なツールです。
しかしながら、SymPy はあくまで古典コンピューター上でのシミュレーターにすぎず、上記で見たとおり、高速に動作するツールではありません。
そこで、OpenQLプロジェクトでは、SymPy のこの量子計算のライブラリ部分を改善したいと思い、一部を改修する活動を行っております。

最後までお読み頂きまして、ありがとうございます。SymPy を使えば、すぐに量子プログラミングが楽しめそうということが、感じていただけたでしょうか？

さて、終わりに少しだけ宣伝させてください。この記事をご覧の方で、PythonによるOSS開発に興味がある、もしくは、量子情報で何かプログラムを作ってみたい、という有志の方がいらっしゃいましたら、OpenQLコミュニティ1に参加8してみませんか。このコミュニティには「量子コンピューターが当たり前に使える時代を自分たちの手で創りたい」というエンジニアや学生さんも参加されています。皆さまのご参加をお待ちしております。

May your Christmas wishes come true!
May happy Quantum Computer world has come!

  1. OpenQLプロジェクトは、量子コンピューターを扱うための様々なライブラリを開発するためのオープンソースプロジェクトです。量子情報、量子コンピューターに興味がある人たちが集うコミュニティを運営しております。 ↩ ↩2

  2. 「量子プログラミング」という言葉そのものが新しい言葉です。現在、量子コンピューターが実用的なものではない状況で、量子コンピューター上で動作する「量子アプリケーション」がどのようなソフトウェアで、それを記述するためのプログラムがどういったものかは未成熟です ↩

  3. http://www.wolframalpha.com/ ↩

  4. SymPy チュートリアル　日本語ノート
    Githubにある日本語チュートリアル ↩

  5. 平成29年度学力テスト中学校数学B-問4をだいたいsympyを使って解く ↩

  6. SymPy は、Anaconda にも対応しております。ご自身の環境にあった方法でインストールしてください。 ↩

  7. 実行環境
    機種名: MacBook
    機種ID: MacBook8,1
    プロセッサ名: Intel Core M
    プロセッサ速度: 1.3 GHz
    プロセッサの個数: 1
    コアの総数: 2
    二次キャッシュ（コア単位）: 256 KB
    三次キャッシュ: 4 MB
    メモリ: 8 GB ↩

  8. 普段はオンラインでコミュニケーション（世間話？）をしながらライブラリ開発を進めています。世界のどこからでもご参加いただけます。年明けから東京で勉強会も開催することになりました。チェックしてみてください。https://openql.connpass.com/ ↩

97
74
1

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
97
74

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?