Qulacs向けJordan Wigner変換のC++実装

※ 時期的にはちょっと遅いですが、折角なので[量子コンピューター Advent Calendar 2023](https://qiita.com/advent-calendar/2023/quantum)にこの記事を追加させていただきました。

# 概要
- Qulacsのソースコードに手を加え、以下の内容の関数を実装してみました
```python3
OpenFermionの`InteractionOperator`クラスのインスタンスを受け取り,
Jordan-Wigner変換をしてQulacsの`Observable`クラスのインスタンスを返す関数.
```
- 実装した関数の中身はほとんどC++ですが、関数自体はPythonで呼べるようにしてあります。
- 基底関数にcc-pVDZを使い、実装した関数を使った場合と、そうでない既存の関数を場合で「Jordan-Wigner変換+`Observable`生成」にかかった時間を比較したら下の表のようになりました。ただし、表の「既存の関数の組み合わせ」は、データ変換の回数が多いので、表はJordan-Wigner変換そのものの厳密な比較にはなっていません。表は「Qulacsユーザーが高度なAPIを使ってJordan-Wigner変換とObservable生成を行う自然なコードを書いたときに実行時間」の比較になります。

| 分子 | Qubit数 | 既存の関数 (秒) | 実装した関数 (秒) |
| ----|---- | ----|---- | 
| H2 | 20 | 0.5088 | 0.0022 | 
| LiH | 38 | 9.6098 | 0.0991 | 
| H2O | 48 | 21.1508 | 0.2058 |
| CH4 | 68 | 129.9373 | 1.6031 |
| CO2 | 84 | 224.3241 | 2.2575 |

# 実装内容のリンク

Qulacsをforkして`jordan-wigner`という名前のブランチを作って、そこに新しいJordan-Wigner変換用の関数を追加しました。↓はそのブランチと`main`のdiffです。

https://github.com/qulacs/qulacs/compare/main...shim98a:qulacs:jordan-wigner



# 想定読者
- QulacsとOpenFermion/OpenFermion-PySCFを組み合わせて使っているときに、Jordan-Wigner変換が遅いと感じている方

そんなにいないと思いますが、少しでもいて、この記事が役に立てば嬉しいです…

# 背景
量子コンピュータを用いた化学計算の関連技術に取り組む際は、パウリ演算子で表されたハミルトニアンを生成するためにJordan-Wigner変換やBravyi-Kitaev変換を使うことがありますよね。こうした変換を自動で行う関数は、[OpenFermion](https://github.com/quantumlib/OpenFermion)や[PennyLane](https://github.com/PennyLaneAI/pennylane)、[Qiskit-Nature](https://github.com/qiskit-community/qiskit-nature)といったライブラリで実装されています[^1]。しかし、これらのライブラリで実装されている変換の関数は実行時間が長く、ハミルトニアン生成で待たされるのがネックでした。勿論、生成したハミルトニアンに対し、古典コンピュータでVQEを状態ベクトルでシミュレーションしたい場合はネックにならない[^2]と思いますが、そうでない稀なケース(ゲート数の見積もりなど)では変換が速いと結構嬉しかったりします。

ということで、今回は量子計算用ライブラリである[Qulacs](https://github.com/qulacs/qulacs)向けに、Jordan-Wigner変換を高速にやる関数を実装しようと思いました。Bravyi-Kitaevは今回は実装していません。

# 実装した関数の使い方
概要にある通り、Qulacs上に
```Python3
OpenFermionの`InteractionOperator`クラスのインスタンスを受け取り,
Jordan-Wigner変換をしてQulacsの`Observable`クラスのインスタンスを返す関数.
```
を実装しました。この関数の名前は`create_observable_from_openfermion_op`としました。開発した`jordan-wigner`ブランチをcloneして、[README](https://github.com/shim98a/qulacs/tree/jordan-wigner#install-python-library-from-source)の通りにrequirementsを満たした状態で`pip install`すれば使えます(少なくとも、筆者のUbuntu22.04LTS on WSL2環境では。他の環境は未実証です)。READMEの内容は元の`main`ブランチと同じです。


水素分子に対する使用例は次のような感じです：
```python3
from openfermionpyscf import generate_molecular_hamiltonian
from qulacs.observable import create_observable_from_openfermion_op

basis = "sto-3g"
multiplicity = 1
geometry = [("H", (0, 0, 0)), ("H", (0, 0, 0.45))]

# type(hamiltonian) is InteractionOperator
hamiltonian = generate_molecular_hamiltonian(geometry, basis, multiplicity)

# Jordan-Wigner変換
observable = create_observable_from_openfermion_op(hamiltonian)
```
上記同様のことを、既存のQulacs/OpenFermionの関数を使ってやろうとすると下のようなコードになると思います：
```python3
from openfermion import jordan_wigner
from openfermionpyscf import generate_molecular_hamiltonian
from qulacs.utils import convert_openfermion_op

basis = "sto-3g"
multiplicity = 1
geometry = [("H", (0, 0, 0)), ("H", (0, 0, 0.45))]

# type(hamiltonian) is InteractionOperator
hamiltonian = generate_molecular_hamiltonian(geometry, basis, multiplicity)

n_qubits = hamilitonian.n_qubits

# この変数observableは厳密には`Observable`クラスではなく`GeneralQuantumOperator`クラス
observable = convert_openfermion_op(jordan_wigner(hamiltonian), n_qubits)
```
importの部分と最後の1行が違いますね。既存のやり方の最後の行では、
```
OpenFermionのInteractionOperatorクラス
     --(OpenFermionでJordan-Wigner変換)--> OpenFermionのQubitOperator
     --(Qulacsのクラスに変換)--------> QulacsのGeneralQuantumOperator
```
とデータ構造の変換が2回行われます。この2回の変換はPythonベースのコードが大部分を占めます。特に、`OpenFermionのJordan-Wigner変換`はPure Pythonコードで、内部でPythonの標準の辞書型に対する操作が(大きなハミルトニアンの場合は)大量に行われます。小さなハミルトニアンでは変換にそんなに時間がかからないですが、ハミルトニアンを大きくすると結構遅く感じます[^3]。Pythonの標準コンテナに対する小さな操作を大量に行うコードは高速化の余地があります。今回は、`Qulacsのクラスに変換`の部分も含め、小さなサイズのデータを大量に動的に生成しているように見えたので、上記の流れをすべてC++化することにしました。なんでもかんでもC++にすれば速くなるわけではないことに注意ですが、今回は速くなりそうと判断しました。


※ 直前のコードの最後の1行はQulacsの別の関数`create_observable_from_openfermion_text`(中身がほとんどC++実装の関数)を使って
```python3
# この場合はobservableは`Observable`クラス
observable = 
create_observable_from_openfermion_text(str(jordan_wigner(hamiltonian)))
```
で置き換えることもできますが、これは`str(jordan_wigner(hamiltonian))`でOpenFermionの`QubitOperator`を文字列に変換する処理を挟むため、`create_observable_from_openfermion_text`を使えば劇的に速くなるというわけではないと思います。一応両方使ってみましたが`convert_openfermion_op`を使う場合と劇的な速度差は感じれらませんでした。

# なぜQulacsを実装対象にしたか

ずばり、パウリ演算子を管理するクラスがC++で書かれているからです。より厳密には、「パウリ演算子の直積の線形結合で表された演算子を管理するクラスがC++で書かれているから」ですかね。Qulacsにおけるそのクラスは`GeneralQuantumOperator`や`Observable`です。今回はハミルトニアン(エルミート演算子)を扱うの`Observable`クラスを生成する関数を作りました

Qulacs以外の有名量子ライブラリは、ハミルトニアンを管理するクラスはPythonで実装されています。Jordan-Wigner変換の部分だけ頑張ってC++で高速にしても
```
OpenFermionのInteractionOperator
    --(C++でJordan-Wigner変換)--> C++のオリジナルのハミルトニアンクラス
    --(C++の標準コンテナに変換)--> C++標準コンテナで書かれたハミルトニアンクラス
    --(Pythonの型に変換)--------> 既存のPythonのハミルトニアンクラス
```
といったような3回の変換を行わないといけなくなり、時間がかかってしまいます。これは、C++の実装をPythonで呼べるようにするときに使う有名ツール「[pybind11](https://github.com/pybind/pybind11)」が、C++の標準コンテナ(`std::vector`や`std::unordered_map`)とPythonの標準型(`list`や`dict`)の間の変換をサポートしていて、Pythonのハミルトニアンクラスに変換する際にいったんC++側で標準のコンテナに変換する必要があるからです。勿論
```
OpenFermionのInteractionOperator
        --(C++でJordan-Wigner変換)--> C++標準コンテナで表されたハミルトニアンクラス
        --(Pythonの型に変換)--------> Pythonのハミルトニアンクラス
```
と2回の変換で済ませる実装も可能です。しかし、C++標準の`std::unordered_map`が遅い[^4]ことが原因で、高速化具合に制限がかかってしまいます(パウリ演算子のデータを管理する際に`std::unoredered_map`をはじめとしたハッシュテーブルが使われます)。また、ライブラリのPythonハミルトニアンをすべて独自のC++クラスで書き直し高速化を狙うのも可能ですが、これは大幅変更になるので大変です。

一方、Qulacsはハミルトニアンのクラスが既にC++で実装されているため、わざわざ`std::unordered_map`を経由する必要もなく
```
OpenFermionのInteractionOperator
    --(C++でJordan-Wigner変換)---------> C++オリジナルのハミルトニアンクラス
    --(C++でQulacsのObservableに変換)---> Observable (Qulacs) 
```
という比較的スッキリした流れにできます。Python側には`Observable`クラスの参照だけが返るので、3回目のハミルトニアンデータの変換は起こりません。尚、理想的には`Observable`クラスの中身のデータ構造をいじれば
```
OpenFermionのInteractionOperator
    --(C++でJordan-Wigner変換)---> Observable(独自実装)
```
という1回の変換で済み、より速くなりますが、既存の`Observable`クラスや`GeneralQuantumOperator`クラスの中身まで変更する必要があります。そうするとやはり変更が多くなり大変なので今回はやっていません。より大きなqubit数の処理に対する需要がでてきたらやってもいいかもしれません。

# 既存の関数との比較
今回実装した関数を使うと、既存のAPIを使う場合と比べてどれくらい速くなるかの比較をしました。比較を行った環境は以下の通りです

- Ubuntu 22.04LTS on WSL2
- Python 3.10.10 / openfermion 1.6.0 / openfermionpyscf 0.5
- Qulacs: 開発したブランチのもの
- CPU: 11th Gen Intel(R) Core(TM) i7-11800H @ 2.30GHz (を1つ)
- メモリ: 32GB(16GB×2) DDR4 (ですがWSL2は全メモリの半分までしか使えないので実質16GBです)


比較実験には以下の2つのスクリプトを使いました
```bench.py
from argparse import ArgumentParser
from time import time

from openfermion import geometry_from_file, jordan_wigner
from openfermionpyscf import generate_molecular_hamiltonian

from qulacs.observable import (
    create_observable_from_openfermion_op,
)
from qulacs.utils import convert_openfermion_op

parser = ArgumentParser()
parser.add_argument("molecular_name", type=str)
args = parser.parse_args()

basis = "ccpvdz"  # basis set
multiplicity = 1  # spin multiplicity

molecular_name: str = args.molecular_name
geometry = geometry_from_file(molecular_name + ".xyz")

# `hamiltonian` is an instance of `InteractionOperator`
hamiltonian = generate_molecular_hamiltonian(geometry, basis, multiplicity)

n_qubits = hamiltonian.n_qubits

# Measure duration taken for the conventional approach
start = time()
obs_before = convert_openfermion_op(jordan_wigner(hamiltonian), n_qubits)
finish = time()
time_before = finish - start

# Measure duration taken for the new function
start = time()
obs_after = create_observable_from_openfermion_op(hamiltonian)
finish = time()
time_after = finish - start


print(
    "| {} | {} | {:d} | {:.4f} | {:.4f} | {:.1f} |".format(
        molecular_name.upper(),
        basis,
        hamiltonian.n_qubits,
        time_before,
        time_after,
        time_before / time_after,
    )
)
```
```run_bench.sh
#!/bin/bash

python bench.py h2
python bench.py lih
python bench.py h2o
python bench.py ch4
python bench.py co2
```
新しい関数を実装したQulacsをpip installした環境で`run_bench.sh`を実行しました。プログラム内で分子のxyzファイルを読み込まれています。今回は使ったのは以下の4つです。jordan-wigner変換の時間測定だけに使いたかったので、xyz座標の値は適当です。
```h2.xyz
H   0.00    0.00    0.00
H   0.00    0.00    0.45
```
```lih.xyz
Li  0.00    0.00    0.50
H   0.00    0.00    -1.50
```
```h2o.xyz
O   0.00    0.00    0.50
H   0.00    0.50    -0.60
H   0.00    -0.50   -0.60
```
```ch4.xyz
C   0.00    0.00    0.00
H   0.50    0.40    0.00
H   -0.50   0.40    0.00
H   0.00    -0.40   0.50
H   0.00    -0.40   -0.50
```
```co2.xyz
C   0.00    0.00    0.00
O   -0.50   0.50    0.00
O   0.50    -0.50   0.00
```

筆者の環境下では次のような結果となりました
| 分子 | Qubit数 | 既存の関数 (秒) | 実装した関数 (秒) |
| ----|---- | ----|---- | 
| H2 | 20 | 0.5088 | 0.0022 | 
| LiH | 38 | 9.6098 | 0.0991 | 
| H2O | 48 | 21.1508 | 0.2058 |
| CH4 | 68 | 129.9373 | 1.6031 |
| CO2 | 84 | 224.3241 | 2.2575 |

「既存の関数(秒)」は既存のQulacs/OpenFermionのAPIを使ってJordan-Wigner変換する方法、「実装した関数(秒)」が今回作ったJordan-Wigner変換の関数でやる方法ですが、明らかに速くなりました。H2の場合だと高速化の実感はあまりないかもですが、ハミルトニアンが大きくなるとかなり差が付きますね。CO2では100倍くらい速くなりました。

尚、この比較はJordan-Wigner変換そのもののプログラムの実行時間の比較ではありません。一度Jordan-Wigner変換をしたハミルトニアンのデータを、Qulacsの`Observable`のデータ構造に変換する処理も含まれています。そのため、この比較は、Jordan-Wigner変換とそのあとの`Observable`クラス生成にかかった合計時間にほぼ等しくなります。

# 今後
気が向けばBravyi-Kitaev変換もやろうかなと思います。

# 参考文献・レポジトリ
- Qulacs
  - Y. Suzuki et al., Qulacs: a fast and versatile quantum circuit simulator for research purpose, Quantum 5, 559 (2021).
- OpenFermion
  - J. R. McClean et al., OpenFermion: the electronic structure package for quantum computers, Quantum Sci. Technol. 5, 034014 (2020).
- PennyLane
  - Ville Bergholm et al., PennyLane: Automatic differentiation of hybrid quantum-classical computations, arXiv:1811.04968 (2018).
- Qiskit, Qiskit-Nature
  - Qiskit contributors, Qiskit: An Open-source Framework for Quantum Computing, doi:10.5281/zenodo.2573505 (2023)
  - Qiskit-Nature https://doi.org/10.5281/zenodo.7828767
- pybind11
  - GitHub repo: https://github.com/pybind/pybind11
  - Docs: https://pybind11.readthedocs.io/en/stable/index.html

[^1]: OpenFermionだと[この辺](https://github.com/quantumlib/OpenFermion/tree/master/src/openfermion/transforms/opconversions)、PennyLaneだと例えば[ここ](https://github.com/PennyLaneAI/pennylane/blob/master/pennylane/fermi/conversion.py)、Qiskit-Natureだと[この辺](https://github.com/qiskit-community/qiskit-nature/tree/main/qiskit_nature/second_q/mappers)で実装されています。
[^2]: VQEのシミュレーションをするにはqubit数に対し指数的な計算量を要する一方で、Jordan-Wigner/Bravyi-Kitaev変換には多項式的な計算量しかかからないからです。
[^3]: Jordan-Wigner変換に要する計算量がqubit数$N$に対して$\mathcal{O}(N^4)$以上だからです。$O(2^N)$程ではありませんが、$N$を2倍にすると変換が16倍になるので、「2倍にしただけなのに変換の時間が一気に長くなった…」と感じます。
[^4]: 「遅い」というのは主観的な感想ですが、例えば[色んなハッシュテーブル実装の比較ベンチマーク](https://martin.ankerl.com/2022/08/27/hashmap-bench-01/)では`std::unordered_map`は下位に入っています。