CUDA-Q Solversを使ったQAOAの実行

# CUDA-Q Solversとは
NVIDIAが開発している量子計算向け（特に量子-古典ハイブリッド計算向け？）のオープンソースSDKとして、**CUDA-Q**があります。CUDA-QはPythonとC++の両方から使うことができます。

https://developer.nvidia.com/cuda-q

さらに、CUDA-Qを基盤として構築された量子-古典ハイブリッド計算の実行用ライブラリが**CUDA-Q Solvers**です[^1]。こちらもPythonとC++から使うことができ、変分量子アルゴリズムを簡単に試すことができます。

https://nvidia.github.io/cudaqx/components/solvers/introduction.html

[^1]:CUDA-Qで書かれたライブラリ群として、[CUDA-QX](https://nvidia.github.io/cudaqx/index.html)があり、CUDA-Q Solversはその1つです。

変分量子アルゴリズムを簡単に試すことができるライブラリには、CUDA-Q Solversの他にも以下のようなライブラリがあります。
- [Qiskit Algorithms](https://qiskit-community.github.io/qiskit-algorithms/)
- [InQuanto](https://docs.quantinuum.com/inquanto/index.html)
- [Grove](https://grove-docs.readthedocs.io/en/latest/index.html)
- [OpenQAOA](https://el-openqaoa.readthedocs.io/en/latest/index.html)
- [scikit-qulacs](https://qulacs-osaka.github.io/scikit-qulacs/index.html)

2025年1月の時点では、以下3つの変分量子アルゴリズムがCUDA-Q Solversでサポートされているようです。
1. Variational Quantum Eigensolver (VQE)
1. Adaptive Derivative-Assembled Pseudo-Trotter VQE (ADAPT-VQE)
1. Quantum Approximate Optimization Algorithm (QAOA) 

本記事では、これらの中から3つ目のQAOAを取り上げ、実際に試してみた結果を紹介します。

# Quantum Approximate Optimization Algorithm (QAOA) 
量子近似最適化アルゴリズム (QAOA)は、組合せ最適化問題に対するヒューリスティックな量子-古典ハイブリッドアルゴリズムです。CUDA-Q Solversの公式ページでは、QAOAの次の4つの特徴が挙げられています。

1. **ハイブリッドアプローチ**: 量子と古典（従来型）のリソースを効率的に活用する
1. **反復的な最適化**: 古典オプティマイザーが回路のパラメーターを調整してエネルギーを最小化する
1. **NISQ互換性**: 現在のノイズが多い量子コンピューターで実行できるよう設計されている
1. **柔軟性**: 量子化学や幅広い最適化問題に適用可能である

https://nvidia.github.io/cudaqx/examples_rst/solvers/qaoa.html

QAOAについては既に多くの解説があるので、詳細はここでは割愛します。

# 最大カット問題をQAOAで解いてみる
CUDA-Q Solversの[公式Example](https://nvidia.github.io/cudaqx/examples_rst/solvers/examples.html)のコードを参考にして、CUDA-Q Solversに実装されているQAOAで小規模な組合せ最適化問題をシミュレーションにより解いてみました。問題にはJijの[Kapiyvaさん](https://zenn.dev/hydrochaeris)の[Goemans-Williamson Algorithmに関する記事](https://zenn.dev/jij_inc/articles/goemans_williamson)に記載されている最大カット問題 (MAX-CUT)を利用しました。

## 最大カット問題とは
最大カット問題は与えられたグラフに対して、エッジの総数が最大となるようなカットを見つける問題です。今回は以下の5つの頂点と6つの辺からなるグラフの最大カットを求めることを考えます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/317358/2e73660b-49f0-4742-130f-ae12231fddd6.png)



最大カット問題については、以下の記事も参考にしてください。

https://qiita.com/Kashalpha/items/43b2589229f5c31d337a

## 1. 最大カット問題の定式化
まずは`cudaq_solvers`を含め必要なライブラリをインポートします。
```python
import cudaq, cudaq_solvers as solvers
import networkx as nx, numpy as np
```
グラフは`networkx`を使って定義します。
```python
G = nx.Graph()
G.add_nodes_from([0, 1, 2, 3, 4])
G.add_edges_from([(0, 1), (0, 4), (1, 2), (1, 3), (2, 3), (3, 4)])
```
グラフが定義できると、対応する最大カット問題のハミルトニアンは`cudaq_solvers`に用意されている`cudaq_solvers.get_maxcut_hamiltonian`関数を利用して得ることができます。
```python
H = solvers.get_maxcut_hamiltonian(G)
```
自分でパウリ演算子に書き下す必要がないので、楽で良いですね。念のため、生成したハミルトニアンを確認してみましょう。
```python
print(H)
```
> [0.5+0j] ZIIIZ
[0.5+0j] IIIZZ
[-3+0j] IIIII
[0.5+0j] IIZZI
[0.5+0j] ZZIII
[0.5+0j] IZIZI
[0.5+0j] IZZII

5つのパウリ演算子の重み付きの和として表されています。CUDA-Qではビッグエンディアンの記法を採用しているので、一番右側の演算子が0番目のノードに作用する演算子となっていると思います。定数部分`[-3+0j] IIIII`に加えて、$ZZ$の形の項が辺の数と同数の6個なので、正しく最大カット問題（=2体のIsing模型）になっていそうです。

## 2. QAOAの実行
解きたい最適化問題のハミルトニアンが得られたので、QAOAを実行していきます。`cudaq_solvers`ではQAOAを実行するための関数として、`cudaq_solvers.qaoa`が用意されているのでこれを使います。`cudaq_solvers.qaoa`では単純なQAOAだけでなく、**counterdiabatic driving**と呼ばれる断熱時間発展の加速手法を利用した発展的なQAOAも利用できます。おそらく、以下の論文で提案された手法です：["Digitized-counterdiabatic quantum approximate optimization algorithm"](https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.4.013141)。今回は折角なので、counterdiabatic drivingありのQAOAを使うことにします。

QAOAではパラメータ付きの量子回路のパラメータを最適化アルゴリズムで最適化していきますが、まずはその回路パラメータの初期値を設定します。QAOAのパラメータの数は`cudaq_solvers.get_num_qaoa_parameters`で計算できます。QAOAの量子回路のレイヤー数を3にして計算してみましょう。
```python
num_layers = 3
parameter_count = solvers.get_num_qaoa_parameters(H,
                                                  num_layers,
                                                  full_parameterization=True,
                                                  counterdiabatic=True)
```
ここで`counterdiabatic=False`とすれば普通のQAOAです。今回の量子回路は5量子ビット、3レイヤーで`parameter_count`を確認すると51となりました。
```python
print(parameter_count)
```
> 51

パラメータの初期値を決める戦略にも色々あると思いますが、とりあえず一様な乱数で設定しました。
```python
init_params = np.random.uniform(-np.pi / 8, np.pi / 8, parameter_count)
```
これで準備ができたので、QAOAのシミュレーションを実行します。
```python
opt_value, opt_params, opt_config = solvers.qaoa(H,
                                                 num_layers,
                                                 init_params,
                                                 max_iterations=1000,
                                                 optimizer="lbfgs",
                                                 full_parameterization=True,
                                                 counterdiabatic=True,
                                                 verbose=True)
 ```
最適化アルゴリズムには[L-BFGS法](https://ja.wikipedia.org/wiki/L-BFGS%E6%B3%95)を利用しました。他にもCOBYLAが利用できるようです。使える最適化アルゴリズムの数が少ないのは課題かもしれません。`verbose`を`True`に設定しているので、途中経過が表示されて、以下のようになりました。
> \<H> = -2.878180228162
\<H> = -2.652994993099
\<H> = -3.128912030294
\<H> = -3.404621580222
\<H> = -3.404621580222
...
\<H> = -4.987770481448
\<H> = -4.987924819149
\<H> = -4.987905369308
\<H> = -4.988105580630

各最適化ステップでの出力状態のエネルギー期待値が表示されていて、ステップが経過していくことにエネルギーが低下していき、基底エネルギー$E_0=-5.0$に近づいていることがわかります。

## 3. 実行結果の確認
`cudaq_solvers.qaoa`では最適化後の量子回路により得られる量子状態のエネルギー、状態の測定結果、QAOAにより推定された最適解が出力されます。実際に今回の問題を解いた結果を見てみましょう。
```python
print('Optimal energy: ', opt_value)
print('Sampled states: ', opt_config)
print('Optimal Configuration: ', opt_config.most_probable())
```
> Optimal energy:  -4.988105580629855
Sampled states:  { 00010:1 01001:446 01010:3 01011:3 01101:46 10010:385 10100:6 10110:110 }
Optimal Configuration:  01001

QAOAは近似アルゴリズムのため、完全な基底状態を出力するわけではないので、Sampled statesには複数のビット列が現れています。ただ出現頻度には偏りがあり、$01001$と$10010$が多くなっています。このどちらも今回の問題では最適解なので、正しく最適化問題が解けていそうです。$01001$と$10010$の間にも頻度の差があり、$01001$がサンプリングされる確率が高かったので、Optimal Configurationは$01001$のほうになっています。

# まとめ
本記事では、NVIDIAの[CUDA-Q Solvers](https://nvidia.github.io/cudaqx/components/solvers/introduction.html)を利用して、最大カット問題をQAOAで解いてみました。今後、開発が進むと、より複雑な最適化問題も解けるようになっていくかもしれません。