Quantum counting by QAE

# Quantum counting

Quantum counting はGroverアルゴリズムを拡張するために使えるアルゴリズムです。
Groverアルゴリズムでは解の個数が通常1個だけであるときを扱います。
解の個数が$k$個であるときでも使えますが、適切な増幅回数は$k$に依存するため、$k$が既知でないといけません。
しかし問題の解がわからないのに個数はわかってるという状況はあまりなさそうです。

そこで、Quantum counting によって解の個数を推定します。[^1]

- qiskit, tutorial (https://qiskit.org/textbook/ch-algorithms/quantum-counting.html)
- wikipedia, Quantum counting algorithm (https://en.wikipedia.org/wiki/Quantum_counting_algorithm)

## 原理

解の個数の情報を、とある量子状態の確率振幅に埋め込みます。
その確率振幅をQAEで推定します。最後に確率振幅から解の個数を計算し、roundで丸め込めばOKです。

ぶっちゃけ、Quantum countingというのはQAEの別名みたいなもんです。

## 具体例

具体例でやってみます。
4 qubit状態(16個ある)の中に、解が5個含まれているとします。
解を全て重ね合わせたものを$\psi_{1} \rangle$, それ以外の状態の重ね合わせを$\psi_{0} \rangle$　とします。

$|\psi_{0} \rangle = \frac{1}{\sqrt{11}}(|0000\rangle + |0001\rangle + |0010\rangle + |0011\rangle + |0100\rangle + |0110\rangle + |1000\rangle + |1010\rangle + |1100\rangle + |1101\rangle + |1110\rangle) $

$|\psi_{1} \rangle = \frac{1}{\sqrt{5}}(|0101\rangle + |0111\rangle + |1001\rangle + |1011\rangle + |1111\rangle) $

としました。

ここで、4 qubit全てを$H$ゲートで重ね合わせた状態$|\psi \rangle$を考えます。それは、上記の状態を使って以下のように書けます。

$|\psi\rangle = \frac{\sqrt{11}}{\sqrt{16}}|\psi_{0} \rangle + \frac{\sqrt{5}}{\sqrt{16}}|\psi_{1} \rangle =  0.829|\psi_{0} \rangle + 0.5590|\psi_{1} \rangle$

解の個数が5個であるという情報が、$\sqrt{5}$のところに表現されています。
よって、$|\psi \rangle$に含まれている$|\psi_{1} \rangle$の確率振幅を推定すればいいことになります。
これはQAEで実現可能です。
[QAA,Grover,QAEを勉強する。](https://qiita.com/notori48/items/2082be87eea0985b01ba)

### 実装

qiskit.aquaのライブラリを活用していきます。

```python
from qiskit import QuantumCircuit
from qiskit.aqua.algorithms import Grover
from qiskit.tools.visualization import plot_histogram
from qiskit import Aer, execute
from qiskit.aqua import QuantumInstance
import numpy as np
```

```python
nqubit = 4
state_preparation = QuantumCircuit(nqubit, name='State preparation')
state_preparation.h(range(nqubit))

oracle = QuantumCircuit(nqubit)
# Oracle
oracle.h([2,3])
oracle.ccx(0,1,2)
oracle.h(2)
oracle.x(2)
oracle.ccx(0,2,3)
oracle.x(2)
oracle.h(3)
oracle.x([1,3])
oracle.h(2)
oracle.mct([0,1,3],2)
oracle.x([1,3])
oracle.h(2)
```
マーキングするためのオラクルは以下のような回路です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/d9b1f033-abb6-3d4d-c00c-4b8e7d16e923.png)

といってもよくわからないので、実際に4qubitの重ね合わせ状態を入れてオラクルを作用させた結果を示します。

`+(0.25-0.00j)*|0000>
+(0.25-0.00j)*|0001>
+(0.25+0.00j)*|0010>
+(0.25-0.00j)*|0011>
+(0.25-0.00j)*|0100>
+(-0.25+0.00j)*|0101>
+(0.25-0.00j)*|0110>
+(-0.25+0.00j)*|0111>
+(0.25-0.00j)*|1000>
+(-0.25+0.00j)*|1001>
+(0.25-0.00j)*|1010>
+(-0.25+0.00j)*|1011>
+(0.25-0.00j)*|1100>
+(0.25-0.00j)*|1101>
+(0.25-0.00j)*|1110>
+(-0.25+0.00j)*|1111>`

$ \frac{1}{\sqrt{16}}(|0000\rangle + |0001\rangle + |0010\rangle + |0011\rangle + |0100\rangle + |0110\rangle + |1000\rangle + |1010\rangle + |1100\rangle + |1101\rangle + |1110\rangle) - \frac{1}{\sqrt{16}}(|0101\rangle + |0111\rangle + |1001\rangle + |1011\rangle + |1111\rangle) =  \frac{\sqrt{11}}{\sqrt{16}}|\psi_{0} \rangle - \frac{\sqrt{5}}{\sqrt{16}}|\psi_{1} \rangle $

となり、解状態だけを反転させる操作が出来ています。

増幅演算子$Q$は以下のようになります。

```python
grover = Grover(oracle=oracle, good_state=['0000'], state_preparation=state_preparation)
grover.grover_operator.draw(output='mpl')
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/7a94950f-1afd-8f85-50a7-1a44abdcb2ed.png)

QAEの回路を立てます。要は増幅演算子$Q$へQPEをする回路です。

```python
from qiskit.aqua.algorithms import AmplitudeEstimation
num_eval_qubits = 4
QAE = AmplitudeEstimation(num_eval_qubits, state_preparation=state_preparation, grover_operator=grover.grover_operator, quantum_instance=Aer.get_backend('qasm_simulator'))
QAE.construct_circuit().draw('mpl')
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/6dedd25a-1ef6-a95c-69d7-3d1e5f853123.png)

QAEを実行し、確率振幅を得ます。確率振幅は戻り値$a$の平方根だそうです。

```python
result = QAE.run()
Amplitude = math.sqrt(result.a_estimation)
print(Amplitude)
```
`0.5555702475835077`
これが、$\frac{\sqrt{5}}{\sqrt{16}} = 0.5590$ の良い推定になっています。

解の個数に換算します。

```python
print((Amplitude*math.sqrt(16))**2)
```

`k ~ 4.938532800000001`

k = 5 ですので、ほぼあたっています。

これで終わりですが、QAEの挙動がブラックボックスすぎるので、少し中を見てみます。
増幅演算子$Q$へのQPEをした生の結果を取ります。

```python
circ = QuantumCircuit(4+num_eval_qubits,4)
circ.append(QAE.construct_circuit(),range(4+num_eval_qubits))
for i in range(num_eval_qubits):
  circ.measure(i,i)

result = execute(circ,backend=Aer.get_backend('qasm_simulator')).result()
plot_histogram(result.get_counts())
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/903b2298-07e3-6561-c39b-003e3f3569c4.png)

ここで、左端の量子ビットが$2^{-4}$の桁です。
最も頻度の多い1100は、$\theta = (1/16+1/8)*\pi$　であることを意味しています。確率振幅とは$\sin(\theta)$の関係にあります。
`math.sin((1/16+1/8)*math.pi) # 1100`
`0.5555702330196022`
これはQAEの関数が返した確率振幅と同じ値ですね。
ちなみに隣のピーク(1011)をとっても、sin関数の周期性から同じ確率振幅が得られます。
なので、どっちを取っても構いません。

# 結論

解の個数が知りたければ、QAEを使って振幅を推定すればすぐ出せる。

[^1]: 過去記事 [Distance calculation と Grover optimization に関する一考察](https://qiita.com/notori48/items/3f6b5f61fc9b710a2437)　で解の個数はQAEでわかると書いてますが、それがこの方法です。（伏線回収）