【四】如何用量子計算機解線性方程組

09-20

【四】如何用量子計算機解線性方程組

來自專欄 Coraline19 人贊了文章

本專欄著力簡單介紹量子信息，量子糾錯的有關原理，

並使用IBM Q實現一些基礎的量子計算。

參考文獻：

Quantum Information Meets Quantum Matter

Bei Zeng, Xie Chen, Duan-Lu Zhou,Xiao-Gang Wen，arXiv: 1508.02595v4.

Quantum Computation and Quantum Information

Michael A.Nielsen, Isaac L. Chuang.

Precise Creation, Characterization, and Manipulation of Single Optical Qubits

Nicholas Peters, Joseph Altepeter, Evan Jeffrey, David Branning, and Paul Kwiat.

Quantum Algorithm Implementations for Beginners

Patrick J. Coles, Stephan Eidenbenz, Scott Pakin, Adetokunbo Adedoyin, John Ambrosiano, Petr Anisimov, William

Casper, Gopinath Chennupati, Carleton Coffrin, Hristo Djidjev, David Gunter, Satish Karra, Nathan Lemons, Shizeng

Lin, Andrey Lokhov, Alexander Malyzhenkov, David Mascarenas, Susan Mniszewski, Balu Nadiga, Dan

O』Malley, Diane Oyen, Lakshman Prasad, Randy Roberts, Phil Romero, Nandakishore Santhi, Nikolai

Sinitsyn, Pieter Swart, Marc Vuffray, Jim Wendelberger, Boram Yoon, Richard Zamora, and Wei Zhu

Experimental quantum computing to solve systems of linear equations

X.-D. Cai1, C. Weedbrook2, Z.-E. Su1, M.-C. Chen1, Mile Gu3,4, M.-J. Zhu1, Li Li1, Nai-Le Liu1,Chao-Yang Lu1, Jian-Wei Pan

Wikipedia Phase Estimation

在本篇文章中，我們將集中介紹一種基礎演算法以及其應用：

量子Fourier變換
相位估計技術
用量子計算機求解線性方程組

量子Fourier變換

量子Fourier變換類似於離散Fourier變換(DFT)。

Definition 4.1 離散Fourier變換

對於 $N$ 點序列 ${displaystyle left{x[n] ight}_{0leq n<N}}$ ，它的離散傅里葉變換（DFT）為

$hat{x}[k]=dfrac{1}{sqrt{N}}sum_{n=0}^{N-1} e^{-ifrac{2pi}{N}nk}x[n] qquad k = 0,1,ldots,N-1.$

離散傅里葉變換的逆變換（IDFT）為

$xleft[n ight]={1 over sqrt{N}}sum_{k=0}^{N-1} e^{ ifrac{2pi}{N}nk}hat{x}[k] qquad n = 0,1,ldots,N-1.$

在數字系統中，由於採樣必須是有限長並且離散的點，所以我們只能採用離散Fourier變換。

實際上，這件事情還更有趣，如果學過固體物理的同學，應該覺得離散Fourier變換的式子很眼熟，實際上離散Fourier變換這套操作我們在處理一維晶格的時候一直在用。它來源於我們所加Born-von Karman邊界條件，它將無窮長的晶鏈截斷成 $N$ 個原胞的系統。每一種離散Fourier變換實際上對應了一種聲子的模式。

對於量子體系，我們實際上可以定義量子Fourier變換。

Definition 4.2 量子Fourier變換

對於 $N$ 個正交基( ${|0 angle,|1 angle,dots,|N-1 angle}$ )的量子系統，定義以下的量子操作

$|j angle=hat{mathcal{F}}{ |k angle} odfrac{1}{sqrt{N}}sum^{N-1}_{j=0}e^{2pi i j k/N}|k angle$

稱為量子Fourier變換。

實際上量子Fourier變換當中 $N=2^n$ ，其中 $n$ 為物理qubit數， $N$ 是計算基底的數量。

在這裡我們可以使用二進位數的表述，定義

$j=j_1j_2dots j_n=j_1 2^{n-1}+j_2 2^{n-2}+dots j_n 2^{0}$

那麼我們可以計算一下一個任意計算基底的Fourier變換

$egin{aligned} |j angle o& dfrac{1}{2^{n/2}}sum^{2^n-1}_{k=0}e^{2pi i j k/2^n}|k angle\ &= dfrac{1}{2^{n/2}}sum^1_{k_1=0}dots sum^1_{k_n=0}e^{2pi ij(sum^n_{l=1}k_l 2^{-l})}|k_1k_2dots k_n angle\ &=dfrac{1}{2^{n/2}}sum^1_{k_1=0}dots sum^1_{k_n=0}prod^n_{l=1}e^{2pi ijk_l 2^{-l}}|l_l angle\ &=dfrac{1}{2^{n/2}}prod^{n}_{l=1}[|0 angle+e^{2pi ij2^{-l}}1 angle]\ &=dfrac{(|0 angle+e^{2pi i0.j_n}|1 angle)(|0 angle+e^{2pi i0.j_{n-1}j_n}|1 angle)dots (|0 angle+e^{2pi i0.j_1dots j_n}|1 angle)}{2^{n/2}} end{aligned}$

其中符號 $0.j_l j_{j+1}dots j_m o j_l/2+j_{l+1}/4+dots j_{m}/2^{m-l+1}$

量子Fourier線路

這其中我們只需要使用一種特殊的量子門，就是 $R$ 門(rotation)

定義 $R_k=egin{pmatrix} 1&0\0 &e^{2pi i/2^k}end{pmatrix}$

則線路可以如上圖所示。

現在我們模擬一個toy model，一個兩qubit的Fourier變換，考慮初態為 $|psi angle=dfrac{1}{sqrt{2}}(|00 angle+|10 angle)$

相應的電路為

IBM-Q電路

一開始使用 $H$ 在第一個qubit上構造疊加態，後面使用了

$T^dagger CX T CX=ControlS$

這一由之前的分解定理可知。

對於 $|00 angle$ ， $j_1=1,j_2=0$ ，

$mathcal{F}{|00 angle}=(|0 angle+|1 angle)(|0 angle+|1 angle)=|00 angle+|01 angle+|10 angle+|11 angle$

$|10 angle$ ， $j_1=1,j_2=0$

$mathcal{F}{|10 angle}=(|0 angle-|1 angle)(|0 angle+|1 angle)=|00 angle-|01 angle+|10 angle-|11 angle$

故最終結果應該為

$mathcal{F}{|00 angle+|10 angle}=|00 angle+|01 angle$

而IBM Q x5反饋的結果為

註：IBMx5的qubit標記是從最後開始的，即最後一位是第一個qubit，以此類推

include "qelib1.inc";qreg q[5];creg c[5];h q[0];h q[0];tdg q[0];cx q[1],q[0];t q[0];cx q[1],q[0];h q[1];measure q[1] -> c[1];measure q[0] -> c[0];

相位估計(Phase Estimation)

相位估計是用來估計任意幺正操作 $U$ 的本徵值的方法，知道了本徵值，我們就可以實現對角化，這個意義是非常重大的。

相位估計電路

相位估計的問題被概括為

${displaystyle U|psi angle =e^{2pi i heta }left|psi ight angle ,0leq heta <1}.$

其中 $|psi angle$ 是 $U$ 的本徵矢量，求解 $heta$ 。相位估計的步驟如下：

在儲存器 $1$ 當中製備疊加態，儲存器的物理qubit數 $n$ 與要估計相位的精度有關。
從儲存器 $1$ 向儲存器 $2$ 施加 $U,U^2,dots,U^{2^n-1}$ 的受控幺正變換

${displaystyle U^{2^{j}}|psi angle =U^{2^{j}-1}U|psi angle =U^{2^{j}-1}e^{2pi i heta }|psi angle =cdots =e^{2pi i2^{j} heta }|psi angle }.$

此時，儲存器 $1$ 其中的態變為

${displaystyle {frac {1}{2^{frac {n}{2}}}}underbrace {left(|0 angle +e^{2pi i2^{n-1} heta }|1 angle ight)} _{1^{st} qubit}otimes cdots otimes underbrace {left(|0 angle +e^{2pi i2^{1} heta }|1 angle ight)} _{n-1^{th} qubit}otimes underbrace {left(|0 angle +e^{2pi i2^{0} heta }|1 angle ight)} _{n^{th} qubit}={frac {1}{2^{frac {n}{2}}}}sum _{k=0}^{2^{n}-1}e^{2pi i heta k}|k angle ,}$

很顯然，此時我們可以看出了，我們需要做逆Fourier變換，來獲取 $heta$

$mathcal{F}^{-1}{{frac {1}{2^{frac {n}{2}}}}sum _{k=0}^{2^{n}-1}e^{2pi i heta k}|k angle}={displaystyle {frac {1}{2^{n}}}sum _{x=0}^{2^{n}-1}sum _{k=0}^{2^{n}-1}e^{2pi i heta k}e^{frac {-2pi ikx}{2^{n}}}|x angle ={frac {1}{2^{n}}}sum _{x=0}^{2^{n}-1}sum _{k=0}^{2^{n}-1}e^{-{frac {2pi ik}{2^{n}}}left(x-2^{n} heta ight)}|x angle .}$

為了說明演算法結果的精度，考慮最大的 $a in[0,2^{t-1}]$ ，使得

$a/2^t=0.b_1dotsb_t$

最為接近實際的 $heta$ ，並考察它們的差值 $delta= heta-a/2^t$ 。

考慮儲存器 $1$ 的態為

${displaystyle {frac {1}{2^{n}}}sum _{x=0}^{2^{n}-1}sum _{k=0}^{2^{n}-1}e^{-{frac {2pi ik}{2^{n}}}left(x-a ight)}e^{2pi idelta k}|x angle otimes |psi angle .}$

測量之後系統處於 $|a angle$ 態上的概率為

${displaystyle Pr(a)=left|leftlangle aunderbrace {left|{frac {1}{2^{n}}}sum _{x=0}^{2^{n}-1}sum _{k=0}^{2^{n}-1}e^{{frac {-2pi ik}{2^{n}}}(x-a)}e^{2pi idelta k} ight|} _{ ext{State of the first register}}x ight angle ight|^{2}={frac {1}{2^{2n}}}left|sum _{k=0}^{2^{n}-1}e^{2pi idelta k} ight|^{2}={egin{cases}1&delta =0\&\{frac {1}{2^{2n}}}left|{frac {1-{e^{2pi i2^{n}delta }}}{1-{e^{2pi idelta }}}} ight|^{2}&delta eq 0end{cases}}}$

若測得的 $delta=0$ ，則 $a=2^n heta$ ，此時直接計算帶估計的 $heta$ .
若測得的 $delta eq 0$ ，則至少要求測得坍縮在對應的態 $|a angle$ 的概率為

$egin{align} Pr(a) &= frac{1}{2^{2n}} left | frac{1- {e^{2 pi i 2^n delta}}}{1-{e^{2 pi i delta}}} ight |^2 && ext{for } delta eq 0 \ &= frac{1}{2^{2n}} left | frac{2 sin left ( pi 2^n delta ight)}{ 2sin( pi delta)} ight |^2 && left| 1-e^{2ix} ight|^2 = 4left| sin (x) ight|^2 \ &= frac{1}{2^{2n}} frac {left | sinleft(pi 2^n delta ight) ight |^2}{| sin( pi delta) |^2} \ &geqslant frac{1}{2^{2n}} frac {left | sinleft(pi 2^n delta ight) ight |^2}{| pi delta |^2} && | sin(pi delta) | leqslant | pi delta | ext{ for } |delta| leqslant frac{1}{2^{n+1}} \ &geqslant frac{1}{2^{2n}} frac {|2 cdot 2^n delta|^2}{| pi delta |^2} && | 2cdot2^n delta | leqslant | sin(pi 2^ndelta) | ext{ for } |delta| leqslant frac{1}{2^{n+1}} \ &=frac {4}{pi^2} end{align}$

用量子演算法解線性方程

有了以上的工具，我們就可以考慮用量子演算法來接線性方程了。但是實際上這種演算法是沒有什麼顯著優勢的，在此我們僅僅做一討論。

對於線性方程組 $Avec{x}=vec{b}$ ，如果 $A$ 是Hermite的，我們可以直接求得

$|x angle=dfrac{A^{-1}|b angle}{||A^{-1}|b angle||}$

如果 $A$ 是非Hermite的，我們可以首先將其構造成Hermite矩陣，考慮新方程組

$egin{pmatrix}0 & A\A^dagger &0end{pmatrix} egin{pmatrix}0 \vec{x}end{pmatrix}=egin{pmatrix} vec{b} \ 0 end{pmatrix}$

滿足我們想要的性質。

對於算符 $A$ ，可以按其本徵態展開為

$A=sum^N_{i=1}lambda_i|lambda_i anglelangle lambda_i|$

同理

$A^{-1}=sum^N_{i=1}lambda_i^{-1}|lambda_i anglelangle lambda_i|$

同時，可以將向量 $|b angle$ 投影在 $A$ 的本徵態上

$|b angle =sum^N_{i=1}eta_i |lambda_i angle,eta_i=langle lambda_i|b angle$

因此有

$|x angle propto A^{-1}|b angle=(sum^N_{i=1}lambda_i^{-1}|lambda_i anglelangle lambda_i|)(sum^N_{j=1}eta_j |lambda_j angle)=sum^N_{i=1}dfrac{eta_i}{lambda_i}|lambda_i angle$

用量子演算法求解線性方程組的步驟如下：

在儲存器中製備 $n$ 個基態與矢量 $|b angle$ ，此時系統的態為

$|b angle |0 angle^{otimes n}=sum^N_{i=1}lambda_i |lambda_i angle|0 angle^{otimes n}$

利用相位估計方法獲取 $A$ 的本徵態，此時系統的態為

$sum^N_{i=1}eta_i|lambda_i angle |i angle$ ， $|i angle$ 是相位估計後儲存器的態

引入一個輔助位qubit，對其做旋轉操作，旋轉的角度與本徵值 $lambda_i$ 有關

$sum^N_{i=1}eta_i|lambda_i angle|i angle|0 angle o sum^N_{i=1}eta_i|lambda_i angle|i angle(sqrt{1-C^2/lambda^2_i}|0 angle+C/lambda_i |1 angle)$

其中 $C$ 與歸一化有關

做逆相位估計，將與相位估計的儲存器的關聯取消

$sum^N_{i=1}eta_i|lambda_i angle(sqrt{1-C^2/lambda^2_i}|0 angle+C/lambda_i |1 angle)$

做後選擇，在輔助儲存器 $|1 angle$ 的基上進行測量

$sum^N_{i=1}eta_i|lambda_i angledfrac{C}{lambda_i}|1 angle propto |x angle$

其係數即為解向量的各個分量。

相關線路如圖

量子線路

我們假定待求的方程為 $A|x angle=|b angle$

其中

$A=egin{pmatrix}1.5&0.5\0.5&1.5end{pmatrix}$

我們選擇這個 $A$ 的原因是

$Exp(Ipi A)=X$

而 $Exp(Ipi A/2)$ 由分解公式易得

當我們選擇

$|b angle=dfrac{1}{sqrt{2}}egin{pmatrix}1\1end{pmatrix}$

解得 $|x angle=dfrac{1}{sqrt{2}}egin{pmatrix}1\1end{pmatrix}$

註：測量是其模方

同理，對於其他的態也能得到正確答案

$|b angle=egin{pmatrix}cos(pi/6)\ sin(pi/6)end{pmatrix}$

$|b angle=egin{pmatrix}-1\ 0end{pmatrix}$ ， $|x angle=egin{pmatrix}3/4\1/4end{pmatrix}$

讀者可以一一驗證之，附qiskit代碼如下

import systry: sys.path.append("../../") # go to parent dir import Qconfig qx_config = { "APItoken": Qconfig.APItoken, "url": Qconfig.config[url]}except: qx_config = { "APItoken":"YOUR_TOKEN_HERE", "url":"https://quantumexperience.ng.bluemix.net/api"}# Useful additional packages import matplotlib.pyplot as plt%matplotlib inlineimport numpy as npimport mathfrom math import pifrom pprint import pprintimport lateximport qiskit.tools.qcvv.tomography as tomofrom qiskit.tools.visualization import plot_histogram,plot_statefrom qiskit.tools.qi.qi import state_fidelity, concurrence, purity, outerfrom qiskit import QuantumCircuit,QuantumProgram,ClassicalRegister,QuantumRegister,QISKitErrorfrom qiskit.tools.visualization import circuit_drawer#set QconfigQ_program = QuantumProgram()Q_program.set_api(Qconfig.APItoken, Qconfig.config[url]) # set the APIToken and API urlpprint(Q_program.available_backends())qr = Q_program.create_quantum_register("qr", 4)cr = Q_program.create_classical_register("cr", 4)backend = Q_program.get_backend_configuration(local_qasm_simulator)coupling=backend[coupling_map]qc=Q_program.create_circuit("qc", [qr], [cr])qc.h(qr[0])qc.h(qr[1])qc.ry(2*pi,qr[2])qc.cx(qr[0],qr[2])qc.ry(-pi/2,qr[2])qc.cx(qr[1],qr[2])qc.rz(pi/2,qr[2])qc.cx(qr[1],qr[2])qc.ry(pi/2,qr[2])qc.swap(qr[0], qr[1]) qc.h(qr[1])qc.t(qr[0])qc.cx(qr[1],qr[0])qc.tdg(qr[0])qc.cx(qr[1],qr[0])qc.h(qr[0])qc.swap(qr[0], qr[1])qc.ry(-pi/8, qr[3])qc.cx(qr[0],qr[3])qc.ry(pi/8, qr[3])qc.cx(qr[0],qr[3])qc.ry(-pi/4, qr[3])qc.cx(qr[1],qr[3])qc.ry(pi/4, qr[3])qc.cx(qr[1],qr[3])qc.swap(qr[0], qr[1])qc.h(qr[0])qc.tdg(qr[1])qc.cx(qr[0],qr[1])qc.t(qr[1])qc.cx(qr[0],qr[1])qc.h(qr[1])qc.swap(qr[0], qr[1])qc.ry(-pi/2,qr[2])qc.cx(qr[1],qr[2])qc.rz(-pi/2,qr[2])qc.cx(qr[1],qr[2])qc.ry(pi/2,qr[2])qc.cx(qr[0],qr[2])qc.h(qr[0])qc.h(qr[1])qc.measure(qr[2],cr[2])run1 = Q_program.execute([qc], backend=local_qasm_simulator)plot_histogram(run1.get_counts("qc"))

$秋山疊疊翠，夜月 extbf{圓圓}明$