叢 爽,丁 嬌,張 坤

(中國(guó)科學(xué)技術(shù)大學(xué)自動(dòng)化系,安徽合肥 230027)

1 引言

量子層析(quantum tomography)是最常用的量子狀態(tài)估計(jì)方法,最早由Stokes于1851年提出[1].量子層析是對(duì)量子狀態(tài)的大量全同副本進(jìn)行多次測(cè)量,作為量子狀態(tài)在某個(gè)方向上的投影(坍縮)的概率,通過對(duì)所獲的概率進(jìn)行逆變換,反解出量子狀態(tài)的密度矩陣[2].量子層析在包括量子計(jì)算機(jī)應(yīng)用在內(nèi)的量子態(tài)制備,以及量子系統(tǒng)控制中的狀態(tài)反饋控制都具有重要意義.一個(gè)n比特量子系統(tǒng)可用密度矩陣ρ ∈Cd×d(d2n)描述其狀態(tài),并且密度矩陣需滿足半正定、單位跡的厄米矩陣約束.該系統(tǒng)的完備測(cè)量次數(shù)為d2,它是隨量子位數(shù)n呈指數(shù)增長(zhǎng).壓縮傳感理論是由Donoho,Candes等人于2004年提出[3–4],它為降低量子狀態(tài)估計(jì)中的測(cè)量次數(shù)問題提供了新的思路和解決問題的新手段.該理論指出:如果信號(hào)本身是稀疏的,或者矩陣信號(hào)是低秩的,那么,通過構(gòu)造一個(gè)滿足一定條件的測(cè)量矩陣,就可以遠(yuǎn)少于完備測(cè)量的次數(shù),將該信號(hào)無損地壓縮到低維空間,再通過求解一個(gè)優(yōu)化問題,精確地恢復(fù)出原始信號(hào)[5].現(xiàn)實(shí)中,人們感興趣的量子系統(tǒng)往往處于純態(tài)或近似純態(tài)[6],這意味量子系統(tǒng)的密度矩陣ρ為低秩r ?d.此時(shí),可以采用壓縮感知理論,只要采用較少的測(cè)量次數(shù),就可以精確地重構(gòu)出ρ.Gross證明采用泡利矩陣進(jìn)行觀測(cè)時(shí),由測(cè)量結(jié)果構(gòu)造的測(cè)量矩陣滿足限制等距特性(restricted isometry property,RIP)條件,最少需要mO(rd log d)?d2個(gè)測(cè)量值,就可以通過求解一個(gè)最優(yōu)化問題,精確地恢復(fù)出密度矩陣[7],并定義采樣率為ηm/d2,m為測(cè)量次數(shù).

在實(shí)際量子測(cè)量過程中干擾是不可避免的,當(dāng)在密度矩陣某些位置的元素中存在稀疏干擾,耦合在測(cè)量b ∈Rm中,這種稀疏干擾使得密度矩陣的重構(gòu)變得困難[8].對(duì)于考慮含有稀疏干擾的量子態(tài)估計(jì)問題,可以將問題轉(zhuǎn)化為:在考慮量子態(tài)ρ半正定、單位跡的厄米矩陣的約束條件的情況下,分解求解密度矩陣ρ的核范數(shù),以及稀疏干擾S的l1范數(shù)的兩個(gè)子問題的優(yōu)化問題[9].2014年,對(duì)于量子狀態(tài)本身存在稀疏干擾的問題,Li首次將交替方向乘子法(alternating direction method of multipliers,ADMM)用于基于壓縮傳感的量子態(tài)重構(gòu)中,解決了同時(shí)包含核范數(shù)和l1范數(shù)的優(yōu)化問題[8].該算法求解過程需要計(jì)算一個(gè)m×d2矩陣的偽逆,計(jì)算復(fù)雜度為O(d6).2016年,Zheng等人將不動(dòng)點(diǎn)算法與ADMM結(jié)合,提出了不動(dòng)點(diǎn)方程的ADMM 算法(fixed point ADMM,FP–ADMM)[10].該算法避免了高維矩陣偽逆的求解,將計(jì)算復(fù)雜度下降到O(md4),提高了密度矩陣的重構(gòu)精度.2017年,Zhang等人提出結(jié)合迭代收縮閾值的算法(iterative shrinkage thresholding algorithm,ISTA)[11],有效地解決帶有量子態(tài)約束的優(yōu)化問題,將計(jì)算復(fù)雜度進(jìn)一步下降到O(md2),并且能夠保證最終的優(yōu)化結(jié)果為全局最優(yōu)解.不過,ISTA算法每次迭代的步長(zhǎng)是固定的,所以算法迭代的收斂速度受到限制.2018年,Zhang等人提出非精確的ADMM 算法(imprecise ADMM,I–ADMM)[12],該算法通過采用近鄰梯度法,來近似求解密度矩陣和稀疏干擾的子問題,算法的估計(jì)精度得到進(jìn)一步的提高,并通過仿真實(shí)驗(yàn)證明I–ADMM算法優(yōu)于ISTA算法.非精確的ADMM算法計(jì)算復(fù)雜度為O(md2).

針對(duì)ISTA求解量子狀態(tài)估計(jì)優(yōu)化問題的不足,本文將改進(jìn)的快速迭代收縮閾值算法(fast iterative shrinkage thresholding algorithm,FISTA)應(yīng)用到量子態(tài)密度矩陣的估計(jì)中.FISTA是一種基于加速算子梯度估計(jì)方法優(yōu)化ISTA的算法,該算法在ISTA的基礎(chǔ)上,加入一個(gè)加速算子,該加速算子由當(dāng)前值和前一次值的線性組合構(gòu)成,以加快算法的迭代速度.本文將FISTA算法應(yīng)用于5個(gè)量子位的仿真實(shí)驗(yàn)中,并且將FISTA算法分別與迭代收縮閾值算法(ISTA)、交替方向乘子法(ADMM)、不動(dòng)點(diǎn)方程的ADMM 算法(FP–ADMM)、非精確的ADMM算法(I–ADMM)4種算法進(jìn)行性能對(duì)比實(shí)驗(yàn),并對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行分析.

2 含有稀疏干擾的量子狀態(tài)估計(jì)問題描述

考慮一個(gè)具有n比特量子系統(tǒng)的密度矩陣ρ,其本身含有稀疏干擾.此時(shí),量子密度矩陣估計(jì)問題可描述為:從測(cè)量結(jié)果b中選取的m個(gè),重構(gòu)出d×d的低秩的、含有稀疏干擾S ∈Rd×d的密度矩陣ρ.構(gòu)造觀測(cè)矩陣為A:Cd×d→Cm,則密度矩陣估計(jì)過程中的測(cè)量結(jié)果可寫為bA(ρ+S).

將密度矩陣估計(jì)問題轉(zhuǎn)化為一個(gè)帶有約束條件的目標(biāo)凸優(yōu)化問題:

其中:γ為權(quán)重因子;ρ為待估計(jì)的密度矩陣;∥·∥?為核范數(shù),定義為∥ρ∥?∑si,si為矩陣奇異值;∥·∥1為l1范數(shù).最小化∥ρ∥?和∥S∥1使密度矩陣低秩,且狀態(tài)本身干擾稀疏.為了簡(jiǎn)化,定義凸集為0,tr(ρ)1,ρHρ},ρH為ρ的共軛轉(zhuǎn)置,ρ ∈C表示滿足量子態(tài)約束,并引入示性函數(shù)

則式(1)可改寫為

對(duì)于帶有可分離的目標(biāo)函數(shù)和線性約束的凸優(yōu)化問題式(2),通過引入增廣拉格朗日函數(shù)

其中:α為懲罰參數(shù),懲罰項(xiàng)可以松弛收斂條件;∥·∥2為l2范數(shù);y ∈Rm為拉格朗日乘子.

此時(shí)問題(3)變?yōu)橐粋€(gè)無約束條件的對(duì)增廣拉格朗日函數(shù)的目標(biāo)優(yōu)化問題:

根據(jù)ADMM迭代框架,可將待優(yōu)化問題(4)分解為2個(gè)子問題:求解量子密度矩陣、稀疏干擾的優(yōu)化問題,以及使約束條件為零的拉格朗日乘子的迭代計(jì)算公式

這樣,就把一個(gè)帶有約束條件的優(yōu)化問題,轉(zhuǎn)變?yōu)闊o約束條件的2個(gè)子問題的凸優(yōu)化問題.通過求解優(yōu)化問題(5),可以求解出含有稀疏干擾的量子密度矩陣估計(jì).這個(gè)優(yōu)化問題求解性能的優(yōu)劣,取決于所采用的優(yōu)化算法.下面將先介紹迭代收縮閾值算法,然后在其基礎(chǔ)上,提出迭代收縮閾值的快速改進(jìn)算法.

3 改進(jìn)的快速迭代收縮閾值算法

3.1 迭代收縮閾值算法

在式(5)所描述的優(yōu)化問題中,l2范數(shù)∥·∥2連續(xù)可微,但是,核范數(shù)∥·∥?和l1范數(shù)∥·∥1不連續(xù)可微,對(duì)其求解比較困難.為了解決∥·∥?不連續(xù)可微問題,將式(5)中的ρk+1式乘以然后令f(ρ)∥A(ρ+Sk)這正是軟閾值函數(shù)(soft thresholding)要解決的優(yōu)化問題的形式,對(duì)其中的平滑項(xiàng)f(ρ)求梯度c1,得到梯度的迭代公式為其中wk是第k次迭代的步長(zhǎng).由于這個(gè)算法的整個(gè)過程相當(dāng)于迭代執(zhí)行軟閾值(soft thresholding)函數(shù),所以把它稱為迭代收縮閾值算法(iterative soft thresholding algorithm,ISTA).

由于量子密度矩陣有ρHρ的約束,即密度矩陣是厄米的,所以梯度也必須滿足厄米矩陣的要求,將梯度ck+1的迭代公式重新定義為

將式(6)代入密度矩陣優(yōu)化的式(5)中,可以得到密度矩陣估計(jì)的迭代公式為

其中:

UMVT為矩陣ck+1的奇異值分解;為軟閾值算子,定義為

為了書寫方便,定義D2wk/α(ck+1)為奇異值收縮算子:

對(duì)于式(5)中的Sk+1的求解,先將其中的Sk+1乘以再令f(S)∥A(ρk+1+S)?b ?這也是軟閾值函數(shù)要解決的優(yōu)化問題的形式,對(duì)其中的平滑項(xiàng)f(S)求梯度d,得到梯度的迭代公式為

由此可得ISTA算法的迭代公式為

從式(10)中,可以看出:ISTA是梯度下降法的一種延伸,每次迭代只是利用當(dāng)前點(diǎn)的信息進(jìn)行梯度估計(jì),然后分別對(duì)密度矩陣以及稀疏干擾進(jìn)行迭代估計(jì),所以算法迭代速度一般.

3.2 快速迭代收縮閾值算法

在ISTA的基礎(chǔ)上,通過分別在密度矩陣ρ和稀疏干擾S的計(jì)算公式中引入加速算子,來加速收斂速度,進(jìn)一步降低密度矩陣的估計(jì)誤差.

首先,引入加速算子z,它是由當(dāng)前值和前一個(gè)值的線性組合構(gòu)成:

其中:xk,xk?1分別代表當(dāng)前值和前一次值;xk?xk?1表示搜索方向;hk表示由當(dāng)前值xk開始,沿著xk?xk?1所構(gòu)成的搜索方向進(jìn)行迭代所需要的步長(zhǎng)因子,j是可調(diào)參數(shù);zk+1表示由當(dāng)前值xk開始,沿著xk?xk?1所構(gòu)成的搜索方向進(jìn)行步長(zhǎng)為hk所得到的下一次值.

然后,利用加速算子,分別代入式(10)中含有ρk和Sk中,重新對(duì)含有稀疏干擾的量子狀態(tài)進(jìn)行估計(jì),由此可以分別得到改進(jìn)后的狀態(tài)進(jìn)行估計(jì)迭代公式為

以及稀疏干擾的估計(jì)迭代公式為

其中:hk?1j×權(quán)重γ可以設(shè)置為wk和j需要根據(jù)實(shí)際具體情況來調(diào)節(jié)參數(shù).

由于FISTA算法是在ISTA的基礎(chǔ)上,通過加入一個(gè)由當(dāng)前值和前一次值的線性組合構(gòu)成的加速算子(11),該加速算子是隨著迭代次數(shù)k的增加,而每次增大hk(xk?xk?1),起到對(duì)當(dāng)前值與前一次值之差的進(jìn)一步的補(bǔ)償作用,因而能夠加快算法的迭代速度.理論上已經(jīng)證明[13]:ISTA和FISTA的計(jì)算復(fù)雜度都為O(md2);ISTA算法的收斂速度為O(1/k),而FISTA的收斂速度為O(1/k2).FISTA的收斂速度比ISTA快1/k倍.

4 數(shù)值仿真實(shí)驗(yàn)及其結(jié)果分析

本節(jié)將對(duì)一個(gè)5量子位系統(tǒng)的密度矩陣,采用FISTA算法進(jìn)行狀態(tài)估計(jì)的數(shù)值仿真實(shí)驗(yàn).做兩個(gè)實(shí)驗(yàn):第1個(gè)實(shí)驗(yàn)為:不同采樣率η下,對(duì)FISTA和ISTA兩種算法進(jìn)行仿真實(shí)驗(yàn),比較仿真實(shí)驗(yàn)的性能結(jié)果;第2 個(gè)實(shí)驗(yàn)為:固定采樣率50%下,對(duì)ADMM,FP–ADMM,ISTA,FISTA和I–ADMM 5種算法進(jìn)行仿真實(shí)驗(yàn),并對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行了性能的比較.

實(shí)驗(yàn)中算法性能的評(píng)估指標(biāo)為:估計(jì)出的密度矩陣ρ與真實(shí)密度矩陣¨ρ之間的歸一化誤差Error:Error

4.1 FISTA和ISTA算法的估計(jì)誤差對(duì)比

本實(shí)驗(yàn)將在采樣率η分別為25%,50%,100%,固定迭代次數(shù)為100次,分別采用FISTA以及ISTA兩種算法,對(duì)5量子位密度矩陣進(jìn)行狀態(tài)估計(jì).

兩種算法涉及的可調(diào)參數(shù)有權(quán)重γ、懲罰參數(shù)α、梯度下降步長(zhǎng)wk,FISTA多一個(gè)加速算子里的可調(diào)參數(shù)j.根據(jù)算法的收斂要求,通過實(shí)驗(yàn)結(jié)果對(duì)參數(shù)的設(shè)置進(jìn)行優(yōu)化調(diào)整,最終的參數(shù)設(shè)置如表1所示.在不同采樣率η下,兩種算法對(duì)密度矩陣的估計(jì)誤差Error隨迭代次數(shù)的變化結(jié)果如圖1所示.

表1 兩種算法對(duì)比實(shí)驗(yàn)中的最優(yōu)參數(shù)設(shè)計(jì)Table 1 Optimal algorithm setting in two experiments

圖1 兩種算法的估計(jì)誤差對(duì)比Fig.1 Comparisons of estimation errors of two algorithms

從圖1的量子狀態(tài)估計(jì)誤差的實(shí)驗(yàn)結(jié)果中可以看出:

1)在狀態(tài)估計(jì)達(dá)到穩(wěn)態(tài)之前的暫態(tài)過程中,FISTA明顯優(yōu)于ISTA.在相同采樣率和相同迭代次數(shù)下,FISTA的估計(jì)誤差一直小于ISTA的估計(jì)誤差.

2)隨著采樣率η的增加,ISTA和FISTA兩種算法的穩(wěn)態(tài)估計(jì)誤差都在降低.相同采樣率下,FISTA達(dá)到的穩(wěn)態(tài)估計(jì)誤差比ISTA的穩(wěn)態(tài)估計(jì)誤差低;采樣率為25%和50%時(shí),FISTA 比ISTA 穩(wěn)態(tài)估計(jì)誤差略低,采樣率為100%時(shí),FISTA的穩(wěn)態(tài)估計(jì)誤差明顯低于ISTA.相同采樣率下,FISTA比ISTA以較低迭代次數(shù)達(dá)到較高量子狀態(tài)重構(gòu)精度.采樣率為25%時(shí),ISTA需迭代46 次(0.1799 s)達(dá)到最低估計(jì)誤差0.0025;FISTA 需迭代20 次(0.1153 s)達(dá)到最低估計(jì)誤差0.0018,迭代次數(shù)明顯更少,估計(jì)誤差明顯更低;采樣率為50%時(shí),ISTA需迭代44次(0.1939 s)達(dá)到不高于0.0010 的估計(jì)誤差0.0010;FISTA 需迭代12 次(0.0867 s)達(dá)到不高于0.0010 的估計(jì)誤差0.0002,FISTA明顯更優(yōu);采樣率為100%時(shí),由圖上曲線直接可以看出,FISTA的估計(jì)誤差隨著迭代次數(shù)的增加在降低,一直低于ISTA的估計(jì)誤差.

3)暫態(tài)達(dá)到的最低估計(jì)誤差低于穩(wěn)態(tài)時(shí)的估計(jì)誤差.可見,并不是迭代次數(shù)越多,狀態(tài)估計(jì)達(dá)到的估計(jì)誤差就越低越好,應(yīng)在估計(jì)誤差達(dá)到最小值時(shí)就停止實(shí)驗(yàn),取當(dāng)時(shí)的最小值為狀態(tài)估計(jì)結(jié)果.

4.2 5種算法的估計(jì)誤差對(duì)比

本實(shí)驗(yàn)中將分別采用5 種算法:ADMM,FP–ADMM,ISTA,FISTA和I–ADMM,對(duì)5個(gè)量子位的密度矩陣進(jìn)行狀態(tài)估計(jì).實(shí)驗(yàn)中,采樣率η固定為50%,迭代次數(shù)選為1000次.ADMM算法中的2個(gè)參數(shù)分別取:梯度步長(zhǎng)τ1τ20.1;FP–ADMM算法中的兩個(gè)參數(shù)分別取:權(quán)重γ1、懲罰參數(shù)α0.04;I–ADMM算法中的3個(gè)參數(shù),分別取:梯度步長(zhǎng)τ10.99,τ20.6,懲罰參數(shù)α1.399.FISTA和ISTA兩個(gè)算法的參數(shù)與第4.1節(jié)中對(duì)應(yīng)的實(shí)驗(yàn)參數(shù)選擇一致.

實(shí)驗(yàn)所獲得的估計(jì)誤差Error隨迭代次數(shù)增加的變化結(jié)果如圖2所示.

圖2 5種算法的估計(jì)誤差對(duì)比Fig.2 Comparisons of estimation errors of five algorithms

從5種算法的量子狀態(tài)估計(jì)誤差的結(jié)果圖2中可以看出:

1)在5種算法中,ADMM算法是最差的,FISTA算法是最好的.同一采樣率η下,本文的FISTA算法達(dá)到最小估計(jì)誤差0.0019的量子狀態(tài)重構(gòu)所需要的最少迭代次數(shù)為13次(0.0882 s),目前最優(yōu)的求解存在稀疏干擾的量子態(tài)估計(jì)的I–ADMM優(yōu)化算法的達(dá)到最小估計(jì)誤差0.0017 所需要的最少迭代次數(shù)也為13 次(0.1077 s),最少迭代次數(shù)相同,但FISTA算法所需迭代時(shí)間0.0882 s明顯低于I–ADMM的所需迭代時(shí)間0.1077 s.

在相同的采樣率η的情況下,分別采用ADMM,FP–ADMM,ISTA,FISTA和I–ADMM 5種算法完成1000次迭代所需時(shí)間分別為:137.4342 s,11.6561 s,2.7510 s,2.7208 s,3.4330 s,FISTA 算法所需時(shí)間最短、速度最快.

由圖2可知,本文的FISTA算法估計(jì)誤差明顯低于ADMM,FP–ADMM,ISTA 3 種算法,和I–ADMM 算法估計(jì)誤差接近,由于FISTA算法迭代時(shí)間比I–ADMM算法快,故FISTA算法最優(yōu).

2)同一采樣率η下,FP–ADMM,ISTA,FISTA和I–ADMM 4種算法達(dá)到的穩(wěn)態(tài)誤差相接近,且都遠(yuǎn)低于ADMM的穩(wěn)態(tài)誤差.

ADMM只是一種求解優(yōu)化問題的計(jì)算框架,將大的全局問題分解為多個(gè)較小、較容易求解的子問題,并通過協(xié)調(diào)子問題的解而得到全局問題的解.每一個(gè)子問題和如何有效求解,需要根據(jù)f(x)和g(z)具體形式來確定.ADMM算法中求解Frobenius范數(shù)項(xiàng)最小時(shí),需要計(jì)算一個(gè)矩陣的偽逆,這樣的求解方法會(huì)導(dǎo)致巨大的運(yùn)算量.ADMM算法的重構(gòu)精度相對(duì)于完備測(cè)量重構(gòu)的精度還有待于進(jìn)一步提高,對(duì)于較高比特的量子系統(tǒng)(例如7比特),ADMM算法重構(gòu)時(shí)間要以天為單位計(jì)算,算法速度仍有待于進(jìn)一步提高.

FP–ADMM通過近鄰算子求解基于壓縮傳感的量子態(tài)估計(jì)優(yōu)化問題的最優(yōu)解滿足的不動(dòng)點(diǎn)方程,再采用迭代的方式求解最優(yōu)解,避免了基于壓縮傳感的量子態(tài)估計(jì)的ADMM算法中的大規(guī)模矩陣偽逆運(yùn)算,從而大幅度地減少在密度矩陣重構(gòu)應(yīng)用中的計(jì)算時(shí)間.

I–ADMM的主要思想是通過近鄰梯度法近似求解子問題,非精確求解密度矩陣和稀疏干擾相關(guān)的子問題從而獲得閉式解,將計(jì)算復(fù)雜度從Li的ADMM的O(d6),Zheng 的FP–ADMM 的O(md4),降低到O(md2).此外,該算法采用可調(diào)步長(zhǎng)更新拉格朗日乘子以加速收斂.

ISTA通過梯度下降和近鄰算子操作得到ADMM子問題的解,再代入ADMM迭代框架得到全局問題的解,運(yùn)算過程涉及最大的計(jì)算量為m×d2的矩陣A和d2×1的向量相乘,計(jì)算復(fù)雜度為O(md2),因此ISTA算法可以大幅度減少計(jì)算時(shí)間和存儲(chǔ)空間.

FISTA是在ISTA基礎(chǔ)上的優(yōu)化,加快了收斂速度,收斂速度由O(1/k)變?yōu)镺(1/k2).

3)ISTA和I–ADMM兩種算法的暫態(tài)達(dá)到的最低估計(jì)誤差低于穩(wěn)態(tài)時(shí)的估計(jì)誤差.可見,并不是迭代次數(shù)越多,狀態(tài)估計(jì)達(dá)到的估計(jì)誤差就越低越好,人們應(yīng)在估計(jì)誤差達(dá)到最小值時(shí)就停止實(shí)驗(yàn),取當(dāng)時(shí)的最小值為狀態(tài)估計(jì)結(jié)果.

5 結(jié)論

本文所提出的FISTA算法,可以高效精確地求解出量子狀態(tài)本身含有稀疏干擾的情況下的估計(jì)值.通過5量子位密度矩陣仿真實(shí)驗(yàn)表明,FISTA算法收斂速度更快,量子狀態(tài)估計(jì)誤差更小,更加有效地提高了求解量子狀態(tài)估計(jì)優(yōu)化問題的效率和狀態(tài)估計(jì)精度.實(shí)驗(yàn)結(jié)果還表明,在量子態(tài)估計(jì)迭代過程中,當(dāng)判斷出獲得估計(jì)的最小值時(shí),就應(yīng)當(dāng)停止迭代.