王芙蓉 楊帆 張亞 李世中 王鶴峰

1) (中北大學(xué)機電工程學(xué)院, 太原 030051)

2) (清華大學(xué)物理系, 北京 100084)

3) (西安交通大學(xué)理學(xué)院, 西安 710049)

在大數(shù)據(jù)時代, 高效的數(shù)據(jù)處理至關(guān)重要, 量子計算具有平行計算能力, 為方便處理數(shù)據(jù)提供了新的解決途徑.本文提出了一個基于奇異值分解的矩陣低秩近似量子算法, 復(fù)雜度為 O [log(pq)].在核磁共振量子計算系統(tǒng)完成了算法的原理演示, 選擇一個 8 ×8 維的圖像矩陣, 實現(xiàn)共振躍遷算法的哈密頓量 H 的時間演化,用量子態(tài)層析法分別讀出密度矩陣的不同成分, 對密度矩陣進行重構(gòu), 保真度為99.84%, 在誤差范圍內(nèi)驗證了本文提出的矩陣低秩近似量子算法的正確性.而通過奇異值分解計算低秩矩陣的經(jīng)典算法的復(fù)雜度是O[poly(pq)], 量子算法與經(jīng)典算法相比, 實現(xiàn)了指數(shù)加速.

1 引 言

隨著互聯(lián)網(wǎng)的技術(shù)不斷提高, 大數(shù)據(jù)時代正在到來.矩陣作為處理大數(shù)據(jù)的基本工具之一, 常被用于解決許多實際問題.大多數(shù)情況下, 矩陣表示在空間和時間復(fù)雜度上會隨著數(shù)據(jù)的規(guī)模呈二次方增長, 這導(dǎo)致數(shù)據(jù)處理困難.因此, 去除冗雜信息, 只保留有用的信息, 構(gòu)造低秩矩陣, 快速高效地近似一個目標矩陣, 可提高機器學(xué)習(xí)和數(shù)據(jù)管理的效果.目前, 低秩矩陣已被廣泛用于模式識別[1]、圖像壓縮和檢測[2,3]、圖像處理[4,5]、人臉識別[6]等應(yīng)用領(lǐng)域.

矩陣的低秩近似是一種稀疏表示形式, 可以在保持原矩陣的諸多性質(zhì)同時, 減少冗余和噪聲, 降低了數(shù)據(jù)的存儲空間和計算復(fù)雜度.常用的矩陣低秩近似方法包括: 低秩矩陣恢復(fù)[7], 主成分分析法(PCA)[8]、非負矩陣分解(NMF)[9]、奇異值分解(SVD)、核映射中的核主成分分析(KPCA)和拉普拉斯特征映射(LE)等.

2014年, Tipping和Bishop[10]采用PCA方法降秩, 他們將圖像分割成8 × 8個不重疊的塊, 進行均勻量化將比特平均分配, 用變換后的前4個主成分對整個圖像進行重構(gòu), 模型似然度最大化后,將數(shù)據(jù)分配給重構(gòu)誤差最小的分量進行圖像編碼,最終的比特率為0.5比特/像素, 誤差為 6.2×10-2.與PCA法相比, KPCA法既可以保持全局屬性又可以獲取非線性信息.2017年, 徐夢珂[11]提出基于二維核主成分分析法的拉普拉斯特征映射(2DKPCA+LE)算法進行人臉識別.在此方法中, 用2D-PCA算法訓(xùn)練樣本去相關(guān)性, 之后利用核主成分分析法提取非線性特征, 最后用拉普拉斯特征映射再一次降秩, 其算法識別率高, 計算復(fù)雜度低.2001年, Tsuge S[12]等提出將NMF[13]應(yīng)用于逐項文檔矩陣中文檔向量的降秩.NMF將非負矩陣分解為兩個非負矩陣, 將分解后的矩陣之一視為基本向量.通過將文檔向量投影到由這些基本向量形成的較低維空間上來實現(xiàn)降秩.大多矩陣的低秩近似方法都避免不了求解特征值或奇異值分解, 因此研究基于奇異值分解的低秩近似方法是提高大數(shù)據(jù)處理效率的一個核心問題.

量子計算基于量子力學(xué)原理, 利用量子疊加性、量子糾纏等特性進行計算, 在處理特定問題時相比經(jīng)典計算能指數(shù)加速.利用疊加態(tài)原理, 將某些類型的數(shù)據(jù)存入量子系統(tǒng)所需要的物理資源也遠小于經(jīng)典系統(tǒng).近年來, 許多求解線性代數(shù)問題的量子算法被提出, 例如量子線性方程組算法[14]、量子支持向量機[15]、量子主成分分析[16]、量子奇異值分解算法[17]等, 也提出了許多量子數(shù)據(jù)存儲、讀取的方案.在某些情況下, 這些量子算法和方案在時間復(fù)雜度或空間復(fù)雜度上, 相比經(jīng)典計算機有著指數(shù)級優(yōu)勢.目前只擁有中等規(guī)模、帶噪聲的量子計算機, 大多數(shù)基于通用量子計算模型的量子算法還無法實現(xiàn).利用有限的物理資源對量子算法進行驗證, 是目前量子計算的重要研究方向.

本文提出了一種量子低秩近似算法, 該方法保留了圖像數(shù)據(jù)的主要成分而去掉那些相對不重要或者由噪聲引入的成分, 利用主成分分析的原理,實現(xiàn)了對圖像數(shù)據(jù)矩陣“奇異值分解”、“奇異值過濾”、“矩陣重構(gòu)”的步驟, 相比于經(jīng)典計算, 本量子算法復(fù)雜度有指數(shù)加速.首次利用該量子算法替代了以往的相位估計算法, 大量地減少了對輔助比特的需求, 成功地實現(xiàn)了對圖像矩陣的低秩近似, 驗證了算法的正確性.

2 理 論

2.1 矩陣的低秩近似與奇異值分解

設(shè) A 是一個秩為r的任意矩陣, 若r遠小于A的行數(shù)和列數(shù), 則稱 A 是低秩矩陣.給定一個秩為r的矩陣 A , 欲求其最近似矩陣該問題可形式化為

任何實矩陣 A ∈Rm×n都可以分解為

其中, U ∈Rm×m和 V ∈Rn×n都是幺正矩陣, Σ 是對角元為從大到小依次排列的奇異值 σi的對角矩陣.對矩陣 A 進行奇異值分解后, 將矩陣 Σ 中的r-k 個最小的奇異值置零獲得矩陣 Σk, 僅保留最大的k個奇異值, 則

(3)式就是(1)式的最優(yōu)解(Eckart-Young-Mirsky定理), 其中 Uk和 Vk分別是(2)式中前k列組成的矩陣[18].

2.2 量子低秩近似算法

假設(shè)對于矩陣 A =[aij]∈Rp×q, 有1個量子操作黑箱(oracle) OA:

即給定i和j, OA能給出矩陣A對應(yīng)的元素值.考慮到一般圖像在量子態(tài)上的編碼方式, 使用振幅編碼制備如下量子態(tài):

以上是將圖像制備成量子態(tài)的過程, 也就是初始化的過程.該步可利用量子隨機存儲器(QRAM)[19,20],采用 O [log(pq)] 個步驟實現(xiàn), 當然也可以采用其他初始化方法[21].

對于矩陣 A 的奇異值分解:

若只選取由大到小前 r′個大于閾值 τ 的奇異值重構(gòu)出低秩近似矩陣 A′, 則有

其中, F作用于矩陣 A 時, 會將第 r′+1 到r個成分奇異值取反, 而其他不發(fā)生改變.初始化完成后的量子態(tài):

則F可以利用相位估計算法[22]和閾值 τ 對應(yīng)的反相算符 Oτ高效實現(xiàn).相位估計算法可以總結(jié)為以下算符:

其中寄存器B和C分別用來存儲厄米矩陣A的本征值的估計值和量子態(tài) | ψ〉 ,是量子傅里葉變換算符的逆算符, H?t是t比特Hadamard門.Oτ可以認為是一個量子過濾器算符,

即對寄存器C中所有量子態(tài)二進制表示小于 τ 的態(tài)相位取反, 而其他則不變.

(7)式需要計算兩個矩陣的差, 早期的酉算子乘積量子計算模型無法直接計算[23], 利用對偶量子計算(DQC)[24-26], 可以在量子計算機上實現(xiàn)(7)式, 對偶量子計算在包括計算空間和輔助比特的大空間中是酉的, 酉演化算符可以寫成:

之后對輔助量子比特進行測量, 當測量結(jié)果為0時, 就得到了非幺正算符 ( U0+U1)/2 作用在量子態(tài)上的結(jié)果.

其中 ρAA?是與厄米矩陣 A A?對應(yīng)的密度矩陣, 二者矩陣元相同, 僅相差一個常系數(shù).

綜上, 量子圖像低秩近似算法可以概括如下.

輸入圖像矩陣 A 對應(yīng)的量子態(tài) | ψA〉 、幺正演化算符、閾值 τ 對應(yīng)的反相算符 Oτ.

輸出A 矩陣的低秩近似矩陣 A′對應(yīng)的量子態(tài)

算法步驟

1) 在3個量子寄存器上準備初始態(tài):

2) 對初始態(tài)運行量子相位估計算法, 其中幺正演化算符為.通過求偏跡制備出密度算符 ρAA?, 再根據(jù)Lloyd等[16]在2014年提出的方法, 用 ρAA?有效地實現(xiàn)算符, 之后得到量子態(tài):

3) 利用DQC實現(xiàn)(7)式.令 U0=I , U1=Oτ2,得到:

其中前 r′個奇異值 σk大于等于設(shè)定的閾值 τ.

4) 運行步驟2)的逆運算, 并測量寄存器R.當測量結(jié)果為0時, 忽略寄存器R和C, 剩余部分就是所需要的低秩近似矩陣 A′對應(yīng)的量子態(tài) | ψA′〉 :

現(xiàn)在分析這一算法的復(fù)雜度.利用QRAM制備初始態(tài) | ψ0〉 , 需要的時間是 O [log(pq)].相位估計算法中, 一般有 t0=O(κ/ε).考慮到目標是奇異值大于 τ 的部分, 則 O (κ)=O(σmax/τ).對于低秩近似問題, 通常要求近似后的矩陣和初始矩陣差別很小, 即通常不隨矩陣維度呈多項式增長.相位估計算法的復(fù)雜度一般是而該算法中相位估計算法中 ε 的大小對于遠離閾值 τ 的奇異值幾乎不產(chǎn)生影響, 只會對接近 τ 的部分產(chǎn)生誤差, 因此常數(shù)誤差不會對大多數(shù)被過濾的成分產(chǎn)生影響, 不影響算法輸出的低秩性.綜上所述, 該量子低秩近似算法的復(fù)雜度為 O [log(pq)].作為對比, 經(jīng)典計算機解決基于奇異值分解的低秩近似問題需要的時間是 O [poly(pq)].本量子算法的復(fù)雜度有指數(shù)提高.

2.3 代替相位估計算法的共振躍遷量子算法

本文的量子低秩近似算法, 在實驗上實現(xiàn)有以下困難.首先是制備初始態(tài) | ψA〉.在有QRAM的情況下, 可以高效的制備初始態(tài), 但目前為止還沒有真正實現(xiàn)QRAM功能的量子器件, 因此要選擇其他方式制備.其次是相位估計算法需要一個額外的寄存器存放本征值的估計值, 需要大量的量子比特, 這也是所有基于相位估計的算法的共同問題,使得目前在硬件上實現(xiàn)該類算法非常困難; 除此,實現(xiàn)密度算符 ρAA?和演化算符也需要大量的輔助比特.在設(shè)計算法時, 通常假設(shè)擁有足夠多的比特資源, 只考慮時間復(fù)雜度的優(yōu)勢, 但目前量子計算機仍在研發(fā)當中, 很難提供足夠多的比特資源.

針對以上問題, 分別提出了對應(yīng)的解決方案:1)利用梯度下降優(yōu)化脈沖(GRAPE)[28,29], 將量子態(tài)從贗純態(tài)演化到初始態(tài) | ψA〉 , 算法直接從初始態(tài)|ψA〉開始, 而制備該初始態(tài)的方式不會影響到本文算法的正確性.2)為了減少量子比特數(shù), 選擇量子共振躍遷算法來替代相位估計算法以減少輔助比特的數(shù)量, 只用一個輔助比特即可實現(xiàn)這一步驟.

相位估計在該問題中的作用是先區(qū)分目標奇異向量(需要保留的)和非目標奇異向量(不需要保留的), 以便之后通過測量進行篩選, 量子共振躍遷算法也可以實現(xiàn)對特定的量子態(tài)的選擇.這里,對量子共振躍遷算法簡單介紹(圖1).

圖1 共振躍遷.一個本征值未知的系統(tǒng)與另一個二能級輔助系統(tǒng)存在相互作用.H s 是左邊系統(tǒng)的哈密頓量, 假設(shè)E0=0.當輔助系統(tǒng)的激發(fā)態(tài)能級 ε 接近任意特征值Ei時, 兩個系統(tǒng)之間會發(fā)生共振躍遷Fig.1.Resonant transition.A system with unknown eigenvalues interacts with another two-level auxiliary system.Hs is the Hamiltonian of the system on the left.Assuming E0=0 , when the excited state energy level ε of the auxiliary system is close to any eigenvalue E i , a resonance transition will occur between the two systems.

在兩個原本獨立的量子系統(tǒng)中, 產(chǎn)生了比較微弱(相比系統(tǒng)本身的能量)的相互作用, 例如兩個分子距離逐漸靠近等, 若兩個系統(tǒng)的一些能級能量接近, 則部分能量會在兩個系統(tǒng)的這些能級上躍遷.基于這一物理現(xiàn)象, Wang[30]設(shè)計了一種算法,并在核磁量子計算平臺上成功計算了2比特低能等效水分子哈密頓量的基態(tài)[31].該算法的步驟如下.

1) 將量子態(tài)初始化到 | 0〉?|Φ〉.

2) 模擬哈密頓量 H 并進行時間演化, H 具體形式由下式給出:

其中, σ 是泡利算符, c是兩個系統(tǒng)的相互作用強度, B 是相互算作算符, 需要求解的系統(tǒng)哈密頓量 為 Hs, 在 HT中 引 入 參 考 點,HT=ε0|0〉〈0|?|Φ〉〈Φ|+|1〉〈1|?Hs, 其中 | Φ〉 是初始態(tài).

3) 演化一段時間后測量第1個輔助比特, 若發(fā)生共振, 則測量結(jié)果可能為1或0.當結(jié)果為1時得到目標態(tài).結(jié)果為0時重復(fù)上一步, 再測量,直到測量結(jié)果為1.

理論上共振躍遷算法可以替代相位估計算法,利用較少的輔助比特實現(xiàn)目標.以圖2(a)的?；諡槔? 利用共振躍遷算法進行數(shù)值模擬, 結(jié)果如圖2(b)—(d)所示.

圖2 共振躍遷量子算法數(shù)值模擬 (a) 3 00×300 的原圖像; (b), (c), (d)分別是前10, 30, 50個奇異值恢復(fù)的圖像.50個奇異值已經(jīng)可以很好地恢復(fù)圖像, 可以辨認?；盏募毠?jié)、文字Fig.2.Numerical simulation of resonance transition by quantum algorithm.Panel (a) is the original image of 300×300.Panels (b), (c) and (d) are the images recovered by the first 10, 30 and 50 singular values respectively.50 singular values can be a good restoration of the image, the details and words of the school emblem are legible.

從數(shù)值模擬的結(jié)果來看, 共振躍遷量子算法可以很好地實現(xiàn)預(yù)期目標.

3 實驗演示及結(jié)果分析

3.1 實驗參數(shù)和步驟

為了減少量子比特數(shù), 首先選擇1個 8 ×8 維的圖像矩陣, 并預(yù)先進行厄米化, 表示該矩陣需要的量子比特數(shù)是3.之后利用共振躍遷算法, 只需要1個輔助比特, 共計需要4量子比特, 可以在核磁共振量子計算平臺進行演示.原始矩陣 M 是由圖3所示數(shù)字1的灰度圖像在Matlab中對應(yīng)的8×8維實矩陣:

圖3 數(shù)字1的灰度圖像, 共 8 ×8 個像素Fig.3.A grayscale image of the number 1, with a total of 8×8pixels.

考慮到輔助比特只有1個, 令閾值τ=(σ1+σ2)/2并對 M 厄米化后得到 A.

本次實驗平臺是核磁共振量子計算平臺, 實驗中使用13C 標記的巴豆酸作為4比特樣品, 溶解在d6-丙酮中, 整個過程中1H 解耦.該分子的結(jié)構(gòu)和參數(shù)如圖4所示.C1—C4表示4個量子位, 選擇C1作為輔助量子位.弱耦合近似下的內(nèi)部哈密頓量為圖中 vi為化學(xué)位移, Jij為第i和第j個核之間的耦合強度.實驗在室溫(T = 296.5 K)下的Bruker DRX 400-MHz譜儀上進行.實驗過程如下.

圖4 1 3C 標記巴豆酸的分子結(jié)構(gòu)和參數(shù).C1是輔助量子比特, C2—C4是目標量子比特.在整個實驗中, 1 H 是解耦的.對角元素和非對角元素是化學(xué)位移和J耦合(單位:Hz).最下面是相干弛豫時間 T 2 (單位: s).Fig.4.Molecular structure and parameters of crotonic acid are labeled 1 3C.C1 is the auxiliary qubit and C2-C4 is the target qubit.1 H is decoupled throughout the experiment.The diagonal and off-diagonal elements are chemical shifts and J coupling (in Hertz).At the bottom is the coherent relaxation time T 2 (in seconds).

第1步制備贗純態(tài).這里使用的是空間平均法[32-34], 即利用梯度磁場和多個幺正演化算符實現(xiàn), 當然也可以采用其他方法[35,36].完成后密度算符有如下形式:

其中單位矩陣I不產(chǎn)生有效信號, 極化度 ? ≈10-5.制備贗純態(tài)后, 利用GRAPE將量子態(tài)從 | 0000〉 演化到需要的態(tài) | 0〉?|ψA〉 , 完成初始態(tài)的制備.

第2步實現(xiàn)共振躍遷算法的哈密頓量 H 的時間演化.利用形狀脈沖能夠?qū)崿F(xiàn)時間演化算符e-iHt0, 其中 H 由(17)式給出, 考慮到這里的目標態(tài)是初始態(tài)的一個很好的近似, 令 c = 0.01 ,B=I?3, Hs=A , ε =1 并且 ε0=σ1-1.將演化算符多次作用在初始態(tài)上, 也就是讓 | 0〉? | ψA〉 在H下演化一段時間.

最后, 利用量子態(tài)層析方法對密度矩陣進行讀出[37-41], 需要用到多個讀出脈沖分別讀出密度矩陣的不同成分, 最后對密度矩陣重構(gòu).

本實驗的簡要流程如圖5所示.

圖5 簡要量子電路圖.制備好贗純態(tài)后, 將經(jīng)過GRAPE優(yōu)化好的形狀脈沖作用于第2個寄存器, 完成初始化.之后將時間演化算符 e -iHt0 作用于初始態(tài)N次, 再測量第1個比特, 結(jié)果為 | 1〉 時 得到目標態(tài)|ψA′〉Fig.5.A brief quantum circuit diagram.After the pseudomorphic state is prepared, the shape pulse optimized by GRAPE is applied to the second register to complete the initialization.Then, the time evolution operator e -iHt0 is applied to the initial state for N times, and the first bit is measured.When the result is | 1〉 , the target state | ψA′〉 is obtained.

3.2 實驗結(jié)果及分析

對實驗讀出數(shù)據(jù)進行處理并重構(gòu)密度矩陣后得到的結(jié)果如圖6所示, 另一個奇異向量可以通過完全相同的方法獲得, 只需要交換厄米化時矩陣M 與 M?的順序.通過第1個奇異值對應(yīng)的奇異向量對圖像進行重構(gòu), 從圖6能看出實驗結(jié)果和理論值符合得很好.

圖6 (a)通過實驗得到的密度矩陣; (b)理論計算的結(jié)果; (c)最大的奇異值恢復(fù)的矩陣對應(yīng)圖像; 由于輔助比特是 | 1〉 時才給出有意義的結(jié)果, 所以這里只考慮其為 | 1〉 的子空間, 忽略其余部分Fig.6.(a) Density matrix obtained by experiment; (b) the result of theoretical calculation; (c) the corresponding image of the matrix with the maximum singular value recovery.Since the auxiliary bit is | 1〉 and gives a meaningful result, only consider subspaces whose values are | 1〉 and the rest are ignored.

通過計算保真度F來定量判斷實驗結(jié)果的準確性:

其中 ρ 和 σ 分別是實驗值和理論值, 得到保真度為99.84%.

實驗結(jié)果與理論計算的誤差主要來源于兩方面.一方面是核磁共振譜儀有一定的漲落, 導(dǎo)致脈沖會有一定的誤差, 這是由于實驗設(shè)備造成的.除此之外, 在制備初始態(tài)、實現(xiàn)時間演化算符時, 均利用了梯度下降算法優(yōu)化脈沖, 每個優(yōu)化后的脈沖與目標存在很小的誤差.綜合以上兩點, 實驗結(jié)果在誤差范圍內(nèi)驗證了本文提出的量子矩陣低秩近似算法的正確性.

4 總 結(jié)

本文提出了一個基于奇異值分解的矩陣低秩近似量子算法.對于矩陣 A ∈Rp×q, 求解其低秩矩陣 A′∈Rp×q, 在經(jīng)典計算機上計算耗時一般是O[poly(pq)], 而本文的量子算法, 在量子計算機上只需要耗時 O [log(pq)] 即可解決該問題, 對比經(jīng)典計算機有指數(shù)加速.

本文第2節(jié)中提出的算法用到了相位估計, 在比特數(shù)較多的情況下, 量子主成分分析、量子推薦算法等量子算法都可以實現(xiàn)同樣的目標, 在時間復(fù)雜度上基本相同.實驗中利用了共振躍遷的算法,在這一類問題中等效替代了相位估計算法, 大大減少了對輔助比特的需求, 只通過1個輔助比特, 保留1個奇異值就較好地恢復(fù)了圖像, 這是其他基于相位估計的算法目前不能實現(xiàn)的.

以數(shù)字1的灰度圖像在Matlab中對應(yīng)的矩陣為原始矩陣, 在核磁量子計算平臺求解了該矩陣的低秩近似并展示了近似后的矩陣對應(yīng)的圖像.實驗結(jié)果和理論值相比保真度為99.84%, 充分證明了本文算法的正確性.

在降維、圖像處理等多個機器學(xué)習(xí)相關(guān)領(lǐng)域,矩陣的低秩近似問題具有重要意義.本文提出的量子算法只需要少量的輔助量子比特即可在O[log(pq)]時間內(nèi)給出一個矩陣的低秩近似, 在解決經(jīng)典算法耗時久的問題的同時大大減少了對量子比特資源的需求, 對量子計算機的實際應(yīng)用有重要的意義.

感謝清華大學(xué)物理系龍桂魯教授、南京師范大學(xué)計算機與電子信息學(xué)院段博佳老師提供的寶貴意見.