Solovay-Kitaevの定理 〜Haskell実装例 gridsynth 〜


こんにちは｜こんばんは。カエルのアイコンで活動しております @kyamaz :frog: [^1] です。
[量子コンピュータ Advent Calendar 2022](https://qiita.com/advent-calendar/2022/quantum) の17日目（2022年12月17日）の投稿です。本稿では、ゲート方式の量子コンピューターについてのみ扱います。

## はじめに

### 本稿の対象者

本稿の対象は次のような方を想定しております。

- 量子コンピューターや量子プログラミングに興味がある方
- プログラミング言語Haskellの知識がある方、興味がある方

### 動作環境

Haskellの環境構築については　@mod_poppo さんの記事[Haskellの環境構築2023](https://zenn.dev/mod_poppo/articles/haskell-setup-2023)を参考にして整えてください。本記事の動作を確認している環境は以下の通りです：

- macOS Ventura (Apple Silicon) 
- GHC 9.2.5
- Cabal 3.8.1.0
- Stack 2.9.1
- Stackage: lts-20.3
- GHCup 0.1.18.0

## Solovay-Kitaevの定理

Solovay-Kitaevの定理は、量子計算の分野で最も重要な理論の１つ[^2]と言われており、任意の１量子ビット演算を基本的な量子演算に分解するときに役立つ定理です。１量子ビットに対する密なゲートセットが与えられたとすると、その演算を使って急速に任意の１量子ビットの演算全体を近似的にカバーできるようになるという性質を示した定理です。[Robert M. Solovay](https://en.wikipedia.org/wiki/Robert_M._Solovay)氏により1995年にメーリングリストで初めて提示され、[Alexei Kitaev](https://en.wikipedia.org/wiki/Alexei_Kitaev)氏が1997年に独自に証明の概要を発表しました。
この定理の存在が、限られた万能量子演算子だけで任意の量子演算を効率的に近似できるということの理論的な根拠となります。

『[The Solovay-Kitaev algorithm](https://arxiv.org/abs/quant-ph/0505030)』にこの定理は次のように書かれております。

> Let $G$ be an instruction set for $SU(d)$, and let a desired accuracy $\epsilon \gt 0$ be given. There is a constant $c$ such that for any $U \in SU(d)$ there exists a finite sequence $S$ of gates from $G$ of length $O(log ^{c}(1/\epsilon))$ and such that $d(U, S) \lt \epsilon$.

$d=2$の時が１量子ビットのユニタリ操作にあたります。
少しわかりづらいですが、翻訳すると「任意のゲート操作$U$に対して、計算誤差が$\epsilon$となり、長さが$O(log ^{c}(1/\epsilon))$の有限なゲート操作列$S$が存在し、その長さのパラメータ$c$は定数となる」となります。そして定義にある定数は $c \approx 4$ であることも示されています。

この定理の実装例を紹介するのが本稿の目的ですので、定理についてはここまでとして、この定理を扱っている投稿を紹介します。参考にしてみてください。

- @kyamaz :frog:の投稿『[量子コンピューターの計算精度とSolovay-Kitaevの定理](https://qiita.com/kyamaz/items/6d2c2708cb3351e57be1)』
（※過去記事のため表示が崩れてしまっています）
- @shnchrさんの投稿『[Solovay-Kitaevの定理を「見える化?」する](https://qiita.com/shnchr/items/ac8b6a0e2ee49ef6b383)』

## Haskell実装の例 gridsynth

Neil J. Ross氏とPeter Selinger氏による論文『[Optimal ancilla-free Clifford+T approximation of z-rotations (arXiv:1403.2975)](https://arxiv.org/abs/1403.2975)』では、アダマール(H)ゲートとTゲート、およびTゲートを２回連続で作用したSゲートの３個の基本ゲートセット$\\{H, T, S\\}$で、与えられたZ軸回転ゲートを任意の$\epsilon$の計算誤差で近似する研究がされています。この論文のauthorであるPeter Selinger氏のサイト『[Exact and approximate synthesis of quantum circuits](https://www.mathstat.dal.ca/~selinger/newsynth/)（以降、Quantum Synthesisのサイト）』には、その成果を実装したHaskellプログラム「gridsynth」が公開されております。ここからは、gridsynth を取り扱っていきます。

### ソースコードからビルドする

gridsynth は、最近のオープンソースらしくはなくGithubにはソースコードが公開されておりません。Quantum Synthesisのサイトにバイナリコードとソースコードのリンクがあります。本稿では、ソースコードからビルドする方法を紹介します。Hackageのページからソースコードのアーカイブをダウンロードしてください。

```
$ curl https://hackage.haskell.org/package/newsynth-0.4.0.0/newsynth-0.4.0.0.tar.gz --output newsynth-0.4.0.0.tar.gz
$ tar xvf newsynth-0.4.0.0.tar.gz
$ cd newsynth-0.4.0.0
```

stack.yaml を作成して、newsynth.cabal の依存ライブラリを編集し、ltsにない依存ライブラリはソースコードをダウンロードします。random パッケージは1.2になってから変更点が多く、その影響で`programs/gridsynth.hs`の変更も必要です。適切に変更できれば、次のビルドを実行すると正常終了します。

```
$ stack build
$ stack exec gridsynth -- -h    # 動作確認のためにコマンドヘルプを表示します。
```

gridsynth　は、Z軸回転の回転角度を引数にとって、誤差$\epsilon$で基本ゲートセット$\\{H, T, S\\}$の積を表示するプログラムです。コマンドヘルプの出力は次の通りです。

```
Usage: gridsynth [OPTION...] <theta>
Arguments:
 <theta>                   z-rotation angle. May be symbolic, e.g. pi/128
Options:
  -h        --help          print usage info and exit
  -d <n>    --digits=<n>    set precision in decimal digits (default: 10)
  -b <n>    --bits=<n>      set precision in bits
  -e <n>    --epsilon=<n>   set precision as epsilon (default: 1e-10)
  -p        --phase         decompose up to a global phase (default: no)
  -f "<n>"  --effort="<n>"  how hard to try to factor (default: 25)
  -x        --hex           output hexadecimal coding (default: ASCII)
  -s        --stats         output statistics
  -l        --latex         use LaTeX output format
  -t        --table         generate the table of results for the article
  -c <n>    --count=<n>     repeat count for --table mode (default: 50)
  -r "
"  --rseed="
"   set optional random seed (default: random)
```

### gridsynth を使ってみる

オプションの説明をします。重要なオプションは、誤差$\epsilon$を指定するオプションで`-d`,`-b`,`-e`になります。その中で`-d`が誤差のマイナス何乗かを指定するためのオプションとして使えます。また、`-p`はグローバル位相を加味して圧縮されるので指定するとよいです。
出力に関しては、`-s`が動作を知るには良い出力形式を提供しています。`-l`はLaTeX形式で出力されるので必要に応じて使ってください。なお、`-x`の出力や複数の結果を出力する`-t`（`-c`と同時に繰り返し回数を指定して利用することをお勧めします）は、あまり用途がないと思います。

引数のZ軸回転の回転角（ラジアン）については、数値で指定します。`pi`を用いて`pi/8`のような指定もできます。幾つか具体例を見てみましょう。（グローバル位相を加味した結果で見ます）

```
$ stack exec gridsynth -- -p pi/2
S
$ stack exec gridsynth -- -p pi/4
T
$ stack exec gridsynth -- -p pi/8
HTHTSHTHTSHTSHTHTHTHTHTSHTSHTSHTHTSHTSHTSHTSHTHTSHTSHTSHTHTHTSHT
SHTSHTHTHTHTHTHTSHTSHTHTHTSHTSHTSHTSHTSHTSHTSHTSHTHTHTHTSHTSHTHT
SHTSHTHTHTHTSHTSHTSHTHTSHTSHTHTHTHTHTSHTHTSHTHTHTSHTSHTSHTSHTHTH
TSHTHTSHTSHTHTHTHTSHTSHTHTSHTSHTHTSHTSHTSHTSHTHTSHTHTHTSHTH
```

Z軸回転が$\pi/2, \pi/4$のときは、それぞれSゲート、Tゲートです。Z軸回転が$\pi/8$のときは、上記のような近似計算になります。このときの近似値を確認するには、`-s`オプションをつけます（実行結果は次のとおり）。

```
$ stack exec gridsynth -- -p -s pi/8
HTH...(上記と同じ結果のため割愛)
Random seed: StdGen {unStdGen = SMGen 9620950084922317958 7233975811747058025}
T-count: 98
Lower bound on T-count: 98
Theta: pi/8
Epsilon: 10^(-10) 
Matrix: roothalf^50 * matrix [[Omega (-4095549) (-8122310) 6324583 25541547,Omega 10250618 (-1318839) (-8385498) 13177724],[Omega (-8385498) (-1318839) 10250618 (-13177724),Omega (-6324583) 8122310 4095549 25541547]]
Actual error: 0.5016660126*10^(-10)
Runtime: 0.024829s
Candidates tried: 1 (0 failed, 0 timed out, 1 succeeded)
Time/candidate: 0.024829s
```

`Epsilon:`で表示されているのが、目標とする精度（誤差）で、実際の誤差は`Actual error:`に出力されます。また、量子回路では、Tゲートの個数が実回路の難易度に影響があります。そこで`T-count:`でその数が出力されます。次の例では、精度を粗くして、`Epsilon: 10^(-4) `の例を見てみましょう。

```
$ stack exec gridsynth -- -p -d 4 -s pi/8
HTHTSHTSHTHTSHTHTSHTHTSHTHTSHTHTHTSHTSHTSHTHTSHTHTHTSHTHTHTSHTHT
SHTHTHTSHTSHTSHTSHTSHTSHTHTSHTHTHTHSS
Random seed: StdGen {unStdGen = SMGen 7456393667667392629 13655418611081506763}
T-count: 39
Lower bound on T-count: 39
Theta: pi/8
Epsilon: 10^(-4)
Matrix: roothalf^21 * matrix [[Omega (-861) (-504) 436 287,Omega 517 (-94) (-384) 637],[Omega 637 (-384) (-94) 517,Omega (-287) (-436) 504 861]]
Actual error: 0.8202046019*10^(-4)
Runtime: 0.005127s
Candidates tried: 1 (0 failed, 0 timed out, 1 succeeded)
Time/candidate: 0.005127s
```

当然ながら、出力される回路も短くなりましたし、`T-count:`も減りました。
皆さまも色々と試してみてはいかがでしょう。

（※今回検証している環境では、精度を高めると計算リソースの不足で計算できない事象を確認しております。ただし、同じ精度でも`-s`オプションを指定しなければ計算できたりするケースもありました。`Epsilon: 10^(-50)`あたりまでは計算できそうです。）

## むすび
本稿では、Solovay-Kitaevの定理のHaskell実装例である gridsynth の利用方法を中心にご紹介しました。ソースコードを見てみると論文の理解もでき、Haskellの勉強にもなります。コードリーディングもお薦めします。

:frog:が力をつけて、どこかで解説できたらと思います


（●）（●）
/""　＿＿""＼ @kyamaz :frog:は、オープンソース・コミュニティ[^1]を通じて、皆さんと共に量子情報・量子コンピューティングの分野で活動しております。今後とも引き続き、どうぞ宜しくお願い致します。

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


[^1]: @kyamaz :frog: が運営する OpenQLプロジェクトは、量子コンピューターを扱うための様々なライブラリを開発するためのオープンソースプロジェクトです。量子情報、量子コンピューターに興味がある人たちが集うコミュニティを運営しております。詳しくは[connpassのサイト](https://openql.connpass.com/)をご覧ください。

[^2]:Christopher M. Dawson氏、Michael Nielsen氏による論文『[The Solovay-Kitaev algorithm](https://arxiv.org/abs/quant-ph/0505030)』の冒頭には "The Solovay-Kitaev (SK) theorem is one of the most important fundamental results in the theory of quantum computation."という記述があります。