量子コンピュータのシミュレータをRustで作る

目的

最近、量子コンピュータを勉強しています。一般的な意味でいう「なんもわからん」状態ですので、エンジニアがいう「完全に理解した」状態になるために、量子コンピュータのシミュレータを作ってみる事にしました。
せっかくなので、Rustで実装することにしました。量子コンピュータもRustもわからない事だらけなので、以下を参考にしていみました。

• qulacs/qulacs: Variational Quantum Circuit Simulator for Quantum Computation Research
• hajifkd/rusq: Quantum computing simulator in Rust

基本、qulacsの簡易版をRustで移植して、とりあえずqulacsのチュートリアルを実装できる程度をめざします。

今回実装したコードは、こちらになります。

量子コンピュータのシュミレータは、量子状態を表す複素数の状態ベクトルに対して、行列で表現される量子ゲートを適用して計算していく処理になります。
状態ベクトルから取り出して、行列積を取るところのロジックがよくわからなくて、苦戦したので、わかったことをまとめてみました。

前提

ある程度の量子コンピュータの知識があることを前提としています。
状態ベクトル、計算基底、量子ビット、量子ビットゲート、量子回路等を知っている前提です。

量子計算

１量子ビットゲート(Hゲート)

3量子ビットの回路に対して、１量子ビットゲートのHゲートを適用することを考えます。

まず1番目の量子ビットにHゲートを適用します。回路図にすると以下になります。

シュミレータで実行させると以下のようになります。

• main.rs

use qustV2::circuit::QuantumCircuit;
use qustV2::gates::single::SingleGateApplicator;

fn main() {
    let n = 3;
    let mut circuit = QuantumCircuit::new(n);
    circuit.H(0);
    circuit.update_quantum_state();
    println!("{}", circuit);
}

❯ cargo run
* Qubit Count  : 3
* Dimension    : 8
* State vector :
000> 0.7071067811865475+0i
001> 0.7071067811865475+0i
010> 0+0i
011> 0+0i
100> 0+0i
101> 0+0i
110> 0+0i
111> 0+0i

Hゲートは、0と1の重ね合わせを作るため適用すると 000>と001>が$\frac{1}{2}$の確率で測定されるようになります。
※001>は、右から1番目の値をとります。この場合、1量子ビット目が1です。
量子コンピュータでは、振幅確率となるため、この場合は、$\sqrt{2}$ $\fallingdotseq$ (0.7071067811865475+0i)となります。
このように量子コンピュータは、量子状態ベクトルに対して、量子ゲートを適用していくことで計算を行います。

2番目の量子ビットにHゲートを適用した場合は、以下となります。

• main.rs

use qustV2::circuit::QuantumCircuit;
use qustV2::gates::single::SingleGateApplicator;

fn main() {
    let n = 3;
    let mut circuit = QuantumCircuit::new(n);
    circuit.H(1);
    circuit.update_quantum_state();
    println!("{}", circuit);
}

Hゲート適用するビットを0から1に変更します。
circuit.H(0)をcircuit.H(1)に変更

* Qubit Count  : 3
* Dimension    : 8
* State vector :
000> 0.7071067811865475+0i
001> 0+0i
010> 0.7071067811865475+0i
011> 0+0i
100> 0+0i
101> 0+0i
110> 0+0i
111> 0+0i

000>と010>が$\frac{1}{2}$の確率で測定されるようになります。

手計算

次に1番目と2番目にHゲートを適用した結果をそれぞれ手計算してみます。

1番目の量子ビットにHゲートを適用する場合

初期状態は、000>です。
つまり、

$$\ket{000} = \ket{0} \otimes \ket{0} \otimes \ket{0}$$

であり、行列表記にすると以下のようになる。

$$\ket{0} = \begin{pmatrix} 1 \\ 0 \\ \end{pmatrix} のため、 \ket{000} = \begin{pmatrix} 1 \\ 0 \\ \end{pmatrix} \otimes \begin{pmatrix} 1 \\ 0 \\ \end{pmatrix} \otimes \begin{pmatrix} 1 \\ 0 \\ \end{pmatrix}$$

$(1 0 0 0 0 0 0 0)^\top$ となります。

この初期状態ベクトルに対して、Hゲートを適用します。

Hゲートの行列は以下となります。

$${H} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \\ \end{pmatrix}$$

1番目に適用するため、状態ベクトルと次元をあわせるため2番目と3番目に単位行列（I）を適用します。
つまり、状態ベクトル全体に適用する行列は、以下のようになります。

$$I \otimes I \otimes H$$

※行列演算は、右から掛けていく形になるため、先頭が一番右となります。

$${I} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix}$$

のため、計算すると

$$I \otimes I \otimes H = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix} \otimes \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix} \otimes \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \\ \end{pmatrix} \\ = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \otimes \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \\ \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & -1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 \\ \end{pmatrix}$$

となります。
この行列を状態ベクトルに適用します。

$$\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & -1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 \\ \end{pmatrix} \begin{pmatrix} 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 \\ 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix}$$

2番目の量子ビットにHゲートを適用する場合

適用するHゲートは以下となります。

$$I \otimes H \otimes I$$

計算すると

$$I \otimes H \otimes I = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix} \otimes \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \\ \end{pmatrix} \otimes \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix} \\ = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 \\ 0 & 0 & 0 & 0 & 1 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & -1 \end{pmatrix}$$

となります。

この行列を状態ベクトルに適用します。

$$\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 \\ 0 & 0 & 0 & 0 & 1 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & -1 \end{pmatrix} \begin{pmatrix} 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 \\ 0 \\ 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix}$$

プログラムの結果と一致しました。

計算ロジック

上記、手計算を一般化して計算ロジックにしてみます。
このロジックは、京都大学数理解析研究所の講究録 1120巻の量子計算機シミュレーションシステムを参考にさせて頂きました。
以下、P113を引用します。
k ビット目に U ゲート（$2^{n} \times 2^{n}$ユニタリ行列U）をかける
このゲートは $2^{n} \times 2^{n}$のユニタリ行列であらわすと、

$$\bigotimes_{i=k+1}^{n-1} I \otimes U \otimes \bigotimes_{i=0}^{k-1} I$$

となる。この行列は以下のような規則的な構造をした非常に疎な行列である。

$$U \bigotimes_{i=0}^{k-1} I = X = \\ \begin{rcases} \begin{pmatrix} u_{11} & & u_{12} & & & & & & \\ & u_{11} & & u_{12} & & & & & \\ & & \ddots & & \ddots & & & \\ & & & u_{11} & & & u_{12} \\ u_{21} & & u_{22} & & & & & & \\ & u_{21} & & u_{22} & & & & & \\ & & \ddots & & \ddots & & & \\ & & & u_{21} & & & u_{22} \end{pmatrix} \end{rcases} 2^{k} \\ \bigotimes_{i=k+1}^{n-1} I \otimes U \bigotimes_{i=0}^{k-1} I = \bigotimes_{i=k+1}^{n-1} I \otimes X = \\ \begin{rcases} \begin{pmatrix} X_(1) \\ & X_(2) \\ & & \ddots & \\ &. & & X_(2^{n-k}) \end{pmatrix} \end{rcases} 2^{n}$$

※筆者注: ここで出てくるXは、Xゲートのことではなく、$U \bigotimes_{i=0}^{k-1} I$をXとしておいている。
1行におけるU行列の要素は2つでそれ以外の要素はすべて0となる。（行列Uが $2^{n-1}$ 個埋め込まれた形をしている。）よってこの $2^{n} \times 2^{n}$ 行列を生成する必要はなく、記憶要領には、 $2 \times 2$ 行列のみを保存し、規則性を利用して、U行列の要素があたる部分のみを繰り返し、計算すれば計算時間は、$O(2^{n})$ となる。

引用終わり。

つまり、1量子ビットゲートの場合は、$2^{n} \times 2^{n}$の行列演算となるため、量子状態ベクトルから2個ずつ取り出して行列演算して更新することを繰り返せばよい。

ソースコード

ソースコードは以下のようになります。
qustV2/state.rs at master · forest1040/qustV2から1量子ビットゲート部分を抽出しています。
（※ロジックをもう少し一般化して、複数量子ビットゲートや複数制御ビットに対応したいところです。。）

apply

pub fn apply(
    &mut self,
    qubits_ctrl: &Vec<usize>,
    qubits_target: &Vec<usize>,
    matrix: &Array2<Complex<f64>>,
) {
    let qsize = qubits_ctrl.len() + qubits_target.len();
    let mut qubits = qubits_ctrl.to_owned();
    let mut qubits_tgt = qubits_target.to_owned();
    qubits.append(&mut qubits_tgt);

    let masks = mask_vec(&qubits);
    for i in 0..self.dim >> qsize {
        let indices = indices_vec(i, &qubits_target, &masks);
        let values = indices.iter().map(|&i| self.states[i]).collect::<Vec<_>>();
        let new_values = matrix.dot(&arr1(&values));
        for (&i, nv) in indices.iter().zip(new_values.to_vec()) {
            self.states[i] = nv;
        }
    }
}

量子ゲートを適用する処理です。

for i in 0..self.dim >> qsize {
}

で量子ゲートを構成する行列の次元で割った分、ループします。
事前にmask_vec()で取得したマスクを適用して、
indices_vec()で、量子状態ベクトルの更新する配列のインデックスを取得します。
1量子ビットゲートの場合は、2個ずつ取り出します。後述の2量子ビットゲートの場合は、4個ずつ取り出します。
ゲートが何番目の量子ビットに適用されるかによって、インデックスを取る位置が変わってきます。
あとはmatrix.dot(&arr1(&values))で行列積を計算して更新します。

量子状態のインデックス

mask_vec()とindices_vec()がロジックの肝となりますがその前にどのような感じでインデックスを取得するか一覧を見てみます。

1量子ビットゲートのインデックス

idx	bi0	bi1
0	0	1
	2	3
	4	5
	6	7
1	0	2
	1	3
	4	6
	5	7
2	0	4
	1	5
	2	6
	3	7

• idx: 量子ビットゲートの開始インデックス
• bi0: 0基底として取得する状態ベクトルのインデックス
• bi1: 1基底として取得する状態ベクトルのインデックス

量子ビットゲートの開始インデックスによって、量子ビットゲートの行列が広がっていく形になります。そのため、状態ベクトルのインデックスも開けて取るようにします。（幅を取るイメージ）
例えば、2番目（idx:1）の量子ビットに対して、Uゲートを適用する場合は、最初のループで、状態ベクトルの0と2の値を取って、行列積を取り0と2をそれぞれ更新します。
次のループで1と3を取ります。

mask_vec

pub fn mask_vec(qubits: &Vec<usize>) -> Vec<usize> {
    let min_qubit_index = qubits.iter().min().unwrap();
    let max_qubit_index = qubits.iter().max().unwrap();
    let min_qubit_mask = 1usize << min_qubit_index;
    let max_qubit_mask = 1usize
        << if qubits.len() > 1 {
            *max_qubit_index - 1
        } else {
            *max_qubit_index
        };
    let mask_low = min_qubit_mask - 1;
    let mask_high = !(max_qubit_mask - 1);
    let mut res = Vec::with_capacity(3);
    res.push(max_qubit_mask);
    res.push(mask_low);
    res.push(mask_high);
    res
}

状態ベクトルのインデックスを取得するためのマスクを作成します。
関数のパラメータのqubitsには、量子ビットゲートのインデックスが指定されます。
複数量子ビットゲートに対応するため、Vec指定可能にしています。
指定された量子ビットゲートのminとmaxを取り、maxがマスクとなります。（幅を取るイメージ）
そして、apply()のループのインデックスで2個ずつ（1量子ビットゲートの場合）取るためのローマスクとハイマスクを取ります。

indices_vec

pub fn indices_vec(
    index: usize,
    qubits_tgt: &Vec<usize>,
    masks: &[usize],
) -> Vec<usize> {
    let mut qubits = qubits_tgt.to_owned();
    qubits.sort();
    let mut res = Vec::with_capacity(qubits.len());
    let mask = masks[0];
    let mask_low = masks[1];
    let mask_high = masks[2];
    let basis_0 = (index & mask_low) + ((index & mask_high) << qubits.len());
    let basis_1 = basis_0 + mask;
    res.push(basis_0);
    res.push(basis_1);

    res
}

let basis_0 = (index & mask_low) + ((index & mask_high) << qubits.len())
が状態ベクトルから0基底のインデックスを取る箇所です。
<< qubits.len())の左シフトは、量子ビットゲートのビット数に比例しています。
状態ベクトルから何個ずつ取るかというのと対応します。
(index & mask_low) + ((index & mask_high)で、apply()のループのインデックスを見て、取っていきます。
let basis_1 = basis_0 + maskが、マスクをかけて幅を取るイメージに対応しています。

参考コード

rusqの場合、indices_vec()辺りが該当コードになります。
しかし、2量子ビットゲートのときのインデックスの取り方が自分が思っているのとは違ったため、このロジックは使用しませんでした。。
コード自体は、Rustっぽくて、すごく綺麗です。
しかし、自分のRust力がまだまだなので、コード自体も理解できていないです。。修行せねば。

qulacsの場合、update_ops_matrix_dense_single.cここら辺のコードになります。
今回は、このqulacsのコードを参考にしています。

長くなったので、一旦1量子ビットゲートまでで公開します。
次は、2量子ビットゲートをアップします。
https://zenn.dev/forest1040/articles/67a7f49701847e