※ 時期的にはちょっと遅いですが、折角なので量子コンピューター Advent Calendar 2023にこの記事を追加させていただきました。

概要

Qulacsのソースコードに手を加え、以下の内容の関数を実装してみました

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です。

想定読者

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

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

背景

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

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

実装した関数の使い方

概要にある通り、Qulacs上に

OpenFermionの`InteractionOperator`クラスのインスタンスを受け取り,
Jordan-Wigner変換をしてQulacsの`Observable`クラスのインスタンスを返す関数.

を実装しました。この関数の名前はcreate_observable_from_openfermion_opとしました。開発したjordan-wignerブランチをcloneして、READMEの通りにrequirementsを満たした状態でpip installすれば使えます(少なくとも、筆者のUbuntu22.04LTS on WSL2環境では。他の環境は未実証です)。READMEの内容は元のmainブランチと同じです。

水素分子に対する使用例は次のような感じです：

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の関数を使ってやろうとすると下のようなコードになると思います：

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の標準の辞書型に対する操作が(大きなハミルトニアンの場合は)大量に行われます。小さなハミルトニアンでは変換にそんなに時間がかからないですが、ハミルトニアンを大きくすると結構遅く感じます³。Pythonの標準コンテナに対する小さな操作を大量に行うコードは高速化の余地があります。今回は、Qulacsのクラスに変換の部分も含め、小さなサイズのデータを大量に動的に生成しているように見えたので、上記の流れをすべてC++化することにしました。なんでもかんでもC++にすれば速くなるわけではないことに注意ですが、今回は速くなりそうと判断しました。

※ 直前のコードの最後の1行はQulacsの別の関数create_observable_from_openfermion_text(中身がほとんどC++実装の関数)を使って

# この場合は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」が、C++の標準コンテナ(std::vectorやstd::unordered_map)とPythonの標準型(listやdict)の間の変換をサポートしていて、Pythonのハミルトニアンクラスに変換する際にいったんC++側で標準のコンテナに変換する必要があるからです。勿論

OpenFermionのInteractionOperator
        --(C++でJordan-Wigner変換)--> C++標準コンテナで表されたハミルトニアンクラス
        --(Pythonの型に変換)--------> Pythonのハミルトニアンクラス

と2回の変換で済ませる実装も可能です。しかし、C++標準のstd::unordered_mapが遅い⁴ことが原因で、高速化具合に制限がかかってしまいます(パウリ演算子のデータを管理する際に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

OpenFermionだとこの辺、PennyLaneだと例えばここ、Qiskit-Natureだとこの辺で実装されています。 ↩
VQEのシミュレーションをするにはqubit数に対し指数的な計算量を要する一方で、Jordan-Wigner/Bravyi-Kitaev変換には多項式的な計算量しかかからないからです。 ↩
Jordan-Wigner変換に要する計算量がqubit数$N$に対して$\mathcal{O}(N^4)$以上だからです。$O(2^N)$程ではありませんが、$N$を2倍にすると変換が16倍になるので、「2倍にしただけなのに変換の時間が一気に長くなった…」と感じます。 ↩
「遅い」というのは主観的な感想ですが、例えば色んなハッシュテーブル実装の比較ベンチマークではstd::unordered_mapは下位に入っています。 ↩