陳柄任，袁淏木，吳涵卿，吳 磊，李 鑫，李曉瑜

(1.建信金融科技有限責(zé)任公司 上海 浦東區(qū) 200120；2.四川元匠科技有限公司 成都 611730；3.電子科技大學(xué)信息與軟件工程學(xué)院 成都 610054)

隨著量子計算的發(fā)展，其在解決特定問題時的量子優(yōu)越性正在體現(xiàn)出來。在理論方面，Shor 算法[1]在解決大素數(shù)分解問題時可以實(shí)現(xiàn)指數(shù)加速，這給現(xiàn)代密碼行業(yè)造成了沖擊。Grover 算法[2]可以二次加速數(shù)據(jù)庫中的搜索問題。在實(shí)驗(yàn)方面，谷歌懸鈴木超導(dǎo)量子計算機(jī)首次實(shí)現(xiàn)量子計算優(yōu)越性。而中國科學(xué)技術(shù)大學(xué)的“祖沖之號”[3]和“九章號”[4]分別在超導(dǎo)和光量子計算機(jī)上進(jìn)一步提升了性能。“祖沖之號”進(jìn)行1.2 小時的取樣任務(wù)需要超算花費(fèi)8 年[3]，而“九章號”的取樣速度比經(jīng)典計算機(jī)快1 024倍[4]。

現(xiàn)代金融學(xué)在解決一些特定任務(wù)時，需要極大的算力。這些算力主要被用于對歷史數(shù)據(jù)的分析、高頻交易、金融衍生品定價、投資組合優(yōu)化以及風(fēng)險管理等領(lǐng)域[5-7]。由于數(shù)據(jù)集的龐大，算法時間復(fù)雜度的略微提升，都會對運(yùn)算時間造成嚴(yán)重影響。除此之外，金融學(xué)中的算法對時間異常敏感，耗費(fèi)幾小時或幾天生成的算法可能不再實(shí)用。因此，這種特性也激發(fā)出金融學(xué)與量子計算的結(jié)合[8-10]。

本文提出了一種基于量子判別分析法（quantum linear discriminant analysis, QLDA）的投資組合優(yōu)化方案，求得無目標(biāo)利潤下馬科維茨均值方差模型[11-12]的最優(yōu)解。該解在所有滿足條件的投資方案中達(dá)到夏普率[13]最大。對比經(jīng)典方法，可實(shí)現(xiàn)準(zhǔn)指數(shù)加速。

1 馬科維茨均值方差模型

投資組合優(yōu)化可以描述為如下問題：投資者將一筆資金進(jìn)行投資，在初期用來購買一些證券，然后在末期賣出。投資者需要在眾多證券中選擇哪些證券進(jìn)行購買，并決定分配在這些證券上資金的份額。投資者有兩個決策目標(biāo)：收益率最高；風(fēng)險最低。馬科維茨均值方差模型[11-12]針對投資組合優(yōu)化問題的建模表示為：

式中，w=(w1,w2,···wN)T，N表示可投資資產(chǎn)數(shù)量；wi表 示投資第i個資產(chǎn)的資金占總資金的份額，其解空間是連續(xù)的，因此這種模型被稱為連續(xù)馬科維茨均值方差模型； Σ代表收益率的協(xié)方差矩陣，大小為N×N維；E(r)為每個資產(chǎn)收益率的期望； μ為預(yù)先決定的收益目標(biāo)。馬科維茨均值方差模型的目的在于在給定收益目標(biāo)后，求得一種資產(chǎn)配置方案使得投資風(fēng)險最小。

對于不設(shè)置收益目標(biāo)，只追求最大利潤風(fēng)險比的投資者，該模型可以轉(zhuǎn)化為：

也可以轉(zhuǎn)化為：

如果利率r代表資產(chǎn)利率減去無風(fēng)險利率的凈值，式(2)的目標(biāo)函數(shù)被稱為夏普率[13]，其代表投資人每多承擔(dān)一分風(fēng)險，就可以多拿到幾分超額報酬。夏普率是衡量一個投資方案優(yōu)劣的重要指標(biāo)之一。式(3)的線性形式更加普遍運(yùn)用于各個優(yōu)化模型中[14-16]，其中 λ ∈[0, 1] 表示風(fēng)險偏好，其值越大，投資策略更加偏向于考慮風(fēng)險。

馬科維茨均值方差模型在不同的市場或假設(shè)下也會有細(xì)微的差別。如在一個只允許做多，不允許做空的市場下，wi被約束在[0, 1]間。而在允許做空的市場下就沒有這個假設(shè)。在小額投資中，wi不認(rèn)為是第i個資產(chǎn)的資金占總資金的份額，而是購買第i個資產(chǎn)的手?jǐn)?shù)。在這種情況下，wi是離散的，這種模型稱為離散馬科維茨均值方差模型。其約束條件將會調(diào)整為。在允許做空市場中wi∈{-1,0,1}， 而在不允許做空的市場中wi∈{0,1}。D為所有資產(chǎn)的總手?jǐn)?shù)。連續(xù)型的馬科維茨均值方差模型被廣泛應(yīng)用于實(shí)際場景，投資機(jī)構(gòu)根據(jù)其結(jié)果精度，配置與精度相對應(yīng)的資金規(guī)模。離散型的馬科維茨均值方差模型更貼近真實(shí)情況，因?yàn)樵诠善笔袌鲋?，投資者必須以“手”為最小單位，進(jìn)行購買。離散模型可以通過增加可能離散解的個數(shù)來逼近連續(xù)模型，但算法的時間復(fù)雜度會大大提升。

求解連續(xù)馬科維茨均值方差模型的經(jīng)典算法包含拉格朗日乘子法[14]、混合灰色關(guān)聯(lián)度法[15]、遺傳算法[16]。這些算法的前提(需要)計算收益率的協(xié)方差矩陣，時間復(fù)雜度為ON2M，其中M為歷史時間節(jié)點(diǎn)數(shù)。在求解拉格朗日乘子法時，需要計算矩陣乘積和求逆，其時間復(fù)雜度為O(N3)。因此這些算法的時間復(fù)雜至少是多項(xiàng)式級別的。而求解離散馬科維茨均值方差模型的確定性經(jīng)典算法是Brute-Force 算法，其時間復(fù)雜度為O(2N)。除此之外，GW 優(yōu)化[17]算法被認(rèn)為有解決離散馬科維茨均值方差模型的潛力，GW 算法可以以多項(xiàng)式時間解決相似的最大割問題，其近似比可以達(dá)到最優(yōu)解的80%[18]。即使對于馬科維茲模型的延伸模型[19]，求解離散問題也是一個NP 完全問題[20-21]。

解決馬科維茨均值方差模型的最著名的兩個量子算法是量子退火（quantum annealing, QA）[22-23]與量子近似優(yōu)化算法（quantum approximation optimization algorithm, QAOA）[24-25]。這兩個算法都通過解決組合優(yōu)化問題來解決離散馬科維茨模型。QA 利用量子的隧穿效應(yīng)，進(jìn)而加速優(yōu)化時間[26-28]。DWave 量子退火機(jī)可以實(shí)現(xiàn)伊辛模型的優(yōu)化。伊辛模型哈密頓函數(shù)表示為：

而將組合優(yōu)化問題式(3)寫為二次型形式后，可以發(fā)現(xiàn)該模型恰好就是伊辛模型。QAOA 將可行解編碼為量子態(tài)后放入量子參數(shù)線路（quantum parameterized circuit, QPC）[29]，然后利用變分量子特征值求解算法，運(yùn)用經(jīng)典的非線性優(yōu)化器求得最優(yōu)參數(shù)，最后對量子態(tài)進(jìn)行測量。除這兩個算法外，文獻(xiàn)[30]基于HHL 算法，設(shè)計出復(fù)雜度為poly(log(NM))的組合優(yōu)化算法[31]，而文獻(xiàn)[32]將馬科維茲模型規(guī)約為二階錐優(yōu)化，然后用量子優(yōu)化算法求解，其時間復(fù)雜度是O(N1.5)。

2 量子線性判別分析

2.1 線性判別分析模型

線性判別分析（linear discriminant analysis,LDA）是一種監(jiān)督學(xué)習(xí)下的降維方法[33]。假設(shè)數(shù)據(jù)集是F:M×N，數(shù)據(jù)表示為 {xi∈RN: 1 ≤i≤M}。這M個數(shù)據(jù)被分在k個類中。記第c個類數(shù)據(jù)的均值為μc， 記所有數(shù)據(jù)的均值為xˉ。則有類內(nèi)散度矩陣SW和 類間散度矩陣SB分別為：

LDA 的目的是尋找一個N×K特征矩陣V，使得降維后的M×K數(shù)據(jù)FV類內(nèi)散度矩陣足夠小，類間散度矩陣足夠大。因此LDA 問題變成求解優(yōu)化問題：

在此將特征矩陣V轉(zhuǎn)化為特征向量w分析這個問題。上述問題的對偶問題可以變?yōu)椋?/p>

這是因?yàn)槭?7)中w可以隨意縮放，對于任意常數(shù)k滿 足J(w)=J(kw)。 因此增加條件wTSWw=1后，依舊與原問題等價。該問題使用拉格朗日乘子法求解：

將上式對w求偏導(dǎo)后，得到KKT 條件：

對比式(7)，上式等價于：

根據(jù)式(7)、式(10)和式(11)得到w就 是最大特征值所對應(yīng)的特征向量。然后進(jìn)行進(jìn)一步轉(zhuǎn)化：

問題轉(zhuǎn)變?yōu)榍蠼舛蛎滋鼐仃囎畲筇卣髦邓鶎?yīng)的特征向量：

因此，如果想求解特征矩陣V，就需要求解式(13)最大的K個特征值對應(yīng)的特征向量v，然后將其轉(zhuǎn)化為向量w后，拼接為V。

2.2 厄米特鏈積

文獻(xiàn)[34]提出了一種方案來求解厄米特矩陣的特征向量，這種方案稱為厄米特鏈積（hermitian chain product, HCP）。HCP 的線路如圖1 所示，與HHL 算法十分相似，其目的是將f1(A1)I(f1(A1))?編碼于量子密度算子內(nèi)。

圖1 HCP 線路圖

量子線路的初態(tài)被儲存在3 個寄存器中，第1 個儲存器儲存單量子比特 |0〉；第2 個寄存器可以設(shè)計任意態(tài)為初態(tài)，在理論計算上沒有區(qū)別；第3個寄存器儲存可變的密度算子 ρ0， ρ0在第一輪中為：

量子線路被分為3 部分，第一部分是相位估計門[35-36]，其哈密頓模擬的哈密頓量就是厄米特矩陣A1。 相位估計門可以將矩陣A1的特征向量與其特征值進(jìn)行糾纏。第二部分是一個控制旋轉(zhuǎn)門，其旋轉(zhuǎn)角度與函數(shù)f1有關(guān)，第三部分為第一部分的逆操作。最后測量第一個量子比特，如果結(jié)果為1，那么第三寄存器的密度算子則為ρ1=f1(A1)I(f1(A1))?。

如果以 ρ1作 為第三寄存器的輸入，以A2和f2作為參數(shù)構(gòu)造HCP 線路，那么第三寄存器的量子態(tài)表示為ρ2=f2(A2)f1(A1)I(f2(A2)f1(A1))?。這與式(13)的形式一致。

2.3 密度矩陣指數(shù)化算法

在量子判別分析[34]中，散度矩陣SW和SB通過量子隨機(jī)存儲器（QRAM[37-38]）編碼為密度算子：

式中，A和B是配平常數(shù)。而在厄米特鏈積中，SW和SB需 要 作 為 哈 密 頓 算 子 e-iSWt和 e-iSBt而 非 密度算子輸入給相位估計門。因此，需要使用量子主成分分析（quantum principal component analysis,QPCA）[39]算法中的密度矩陣指數(shù)化算法（density matrix exponentration, DME）將密度算子轉(zhuǎn)化為哈密頓量。

DME 算法的線路如圖2 所示，其輸入態(tài)包含n個 輔助態(tài)ρ（ 被指數(shù)化的態(tài)）和一個目標(biāo)態(tài) σ（被作用于哈密頓模擬e-iρt的態(tài)）。n越大模擬精度越高。其時間復(fù)雜度為， ?代表誤差。

圖2 DME 算法線路圖

整個電路由n個以交換門為哈密頓量的模擬組成，模擬時間為 Δt=t/n。S表示量子交換門。該算法基于公式：

接著利用Hadamard 引理可以得到：

因此利用交換門的哈密頓模擬可以實(shí)現(xiàn)密度算子的指數(shù)化。最后根據(jù)Trotter 公式[40]將上述算法遞歸n次后，對第一個量子比特取偏跡后為e-iρtσeiρt。因此通過復(fù)制多個密度算子，可以將其指數(shù)化。

相位估計算法使用控制下的哈密頓模擬，DME 算法也證明出，將交換門S變?yōu)榭刂平粨Q門就可以實(shí)現(xiàn)控制哈密頓模擬的操作[39]。

3 基于QLDA 的量子投資組合優(yōu)化

3.1 模型規(guī)約

首先將式(7)按照式(2)中的均值方差構(gòu)造，其模型變?yōu)椋?/p>

因?yàn)樵撃Ｐ蛯⒔鈝編碼為量子態(tài)的概率幅，受量子計算限制，其模平方和為1。為方便書寫，記R=E(r)E(r)T。 這里E(r)表 示凈利率。然后將R，Σ和R1/2Σ-1R1/2譜分解得到：

式中，xi是 時間節(jié)點(diǎn)i時N個股票的利潤向量。根據(jù)譜分解，不難看出這3 個矩陣是半正定的厄米特矩陣。使用QLDA 方法可以求解式(19)。接著將式(19)規(guī)約為式(2)，并證明以下定理。

定理 如果已知式(19)的解，那么在多項(xiàng)式時間內(nèi)可以求得式(2)的解。

在證明之前首先引入模型：

進(jìn)而有如下引理。

引理1 如果有解w滿足式(19)，那么w或 -w其中之一滿足式(21)。

證明：使用反證法進(jìn)行證明，如果w或 -w都不滿足式(21)，那么存在w′滿 足，使得下列兩式之一成立：

將不等式左右取平方，有：

則w不滿足式(19)條件的最大值，與原假設(shè)矛盾。

引理2 如果有解w滿足式(21)，那么有解u=滿足式(2)。

由于目標(biāo)函數(shù)的縮放性質(zhì)，解乘以一個正常數(shù)時，取值不變。因此：

之后，記：

根據(jù)式(25)和式(26)，式(28)滿足：

3.2 基于QLDA 的投資組合優(yōu)化算法

結(jié)合第2 部分內(nèi)容，使用QLDA 計算投資組合問題。

1）使用QRAM 按照式(20)將矩陣R與 Σ編碼為密度算子。其構(gòu)造方法由文獻(xiàn)[30]給出。

2）使用DME 算法將密度算子轉(zhuǎn)化為關(guān)于R與Σ的控制哈密頓算子，以此構(gòu)造HCP 線路。

3）使用HCP 算法，其鏈長為2，參數(shù)為f1(x)=x-1/2，A1=Σ，f2(x)=x1/2，A2=R。結(jié) 果 得到形式為R1/2Σ-1R1/2的密度算子。

4）使用DME 算法將R1/2Σ-1R1/2密度算子轉(zhuǎn)化為控制哈密頓算子，并以此構(gòu)造相位估計算法，并且以R1/2Σ-1R1/2作為密度算子輸入。測量占比最多的量子態(tài)則對應(yīng)最大特征值的特征向量v。

5）使用DME 算法將R轉(zhuǎn)化為控制哈密頓算子，并構(gòu)造HHL 線路，輸入的函數(shù)為f1(x)=x-1/2，輸入的量子態(tài)為v。 輸出的量子態(tài)w滿足R1/2w=v。

6）使用Swap-Test 方法[41]或經(jīng)典方法計算vTE(r)≥0， 然后用經(jīng)典方法計算最優(yōu)配置。

資產(chǎn)配置的信息在這個過程中被提取在概率幅中，如何將其概率幅提取為經(jīng)典信息是量子計算的解決難題。一種直觀的方法是對量子態(tài)進(jìn)行M次測量，然后統(tǒng)計每個基態(tài)的次數(shù)Mi。 那么有|vi|=。文獻(xiàn)[30]給出一種判斷vi的符號的方法，利用E(ri)的 符號來判斷vi的符號，它的意義是對于正利潤股票進(jìn)行做多，對負(fù)利潤股票進(jìn)行做空。這在大多數(shù)情況是正確的，但是對于相關(guān)性比較強(qiáng)的股票組合，會存在誤差。配置態(tài) |v〉的其他作用包括計算風(fēng)險 〈v|Σ|v〉或 者計算 〈v|v′〉來衡量其他投資策略的優(yōu)劣。這兩個算法無須測量都可以使用Swap-Test 方法[41]來計算。

在投資組合模型中，R只有一個非0 特征值，而 Σ也不一定滿秩，也可能存在多個0 特征值。為了避免這些0 特征值帶來的干擾，本算法在控制旋轉(zhuǎn)門中規(guī)定一個下限κ ，小于κ 的特征值將不會引起目標(biāo)比特的翻轉(zhuǎn)，進(jìn)而在最終測量結(jié)果為1 時，不會出現(xiàn)小于 κ特征值的基矢。特別地，在這種情況下，一個帶有0 特征值矩陣A的逆表示為：

式中， λi為A的 特征值；ei為 λi的特征向量。

3.3 性能分析

基于QLDA 算法可以實(shí)現(xiàn)準(zhǔn)指數(shù)加速?？紤]誤差項(xiàng)，其總時間復(fù)雜度與QLDA 一致，為O(poly(log(MN))/?3)。其中將數(shù)據(jù)存儲到QRAM 的時間復(fù)雜度是O(log(NM))；HCP 算法不計入誤差項(xiàng)其時間復(fù)雜度也是O(log(NM))；DME 算法的時間復(fù)雜度是O(n)， 因?yàn)榕cN無關(guān)，被認(rèn)為是常數(shù)項(xiàng)；相位估計算法的時間復(fù)雜度是O(poly(log(MN)))；HHL 算法的時間復(fù)雜度是O(log(NM))；Swap-Test 算法的時間復(fù)雜度是O(log(N))。因此所有的量子算法可以實(shí)現(xiàn)指數(shù)加速。但在最后一步，還需要將vi求 和，這一步的復(fù)雜度是O(N)?？紤]到目前應(yīng)用投資組合資產(chǎn)種類小于100，即使將這個值擴(kuò)大到 106，普通的經(jīng)典計算機(jī)也可以瞬間完成，因此，這部分時間可以忽略不計。

在QRAM 中，算法需要提前將數(shù)據(jù)存儲在經(jīng)典RAM 中，然后使用Oracle 查詢獲取感興趣的數(shù)值，這一部分的準(zhǔn)備時間是O(NM)且是一次性的，對比其他算法也可以實(shí)現(xiàn)二次加速。

表1 描述了不同投資組合優(yōu)化問題算法的時間復(fù)雜度。由于部分算法分析時間復(fù)雜度沒有考慮誤差項(xiàng)，因此在比較時認(rèn)為誤差項(xiàng)為常數(shù)項(xiàng)。數(shù)，即迭代線路的深度。從表中可見QLDA 算法與HHL 算法維持在同一水平。HHL 方法雖然沒在文中特別說明，但也需要對向量上的每個元素進(jìn)行求和，以滿足約束條件。此表沒有對消耗量子比特數(shù)目進(jìn)行對比，這是因?yàn)镼LDA 和HHL 需要使用Oracle 利用輔助比特構(gòu)造QRAM，而二階錐優(yōu)化、QAOA、QA 需要重復(fù)算法，其消耗的量子比特難以估算。

表1 不同投資組合優(yōu)化算法的對比

在制備QRAM 過程中，QLDA 與HHL 算法消耗的量子比特是一致的。在算法過程中，HHL 需要求解 3N個線性方程組，其HHL 線路需要的量子比特大致是HCP 的3 倍。因此QLDA 算法相比HHL 更加經(jīng)濟(jì)。

4 結(jié) 束 語

本文基于QLDA 的HCP 算法與QPCA 的DME 算法設(shè)計了無預(yù)先獲利目標(biāo)下的量子投資組合優(yōu)化算法，算法可以求得滿足最大夏普率的投資組合方案解，并且其算法的時間復(fù)雜度約為poly(log(NM))，相比于經(jīng)典解決方案可以達(dá)到準(zhǔn)指數(shù)加速。這在投資組合優(yōu)化場景中，可以使投資組合優(yōu)化算法所考慮的股票數(shù)量大大提高。

然而，本算法和諸多基于DME 的算法在實(shí)施時有著局限性。DME 算法需要n>100時，才有逼近指數(shù)化算子的效果。最后需要重復(fù)n次HCP 實(shí)驗(yàn)來獲得n個R1/2Σ-1R1/2算子才可以實(shí)施相位估計算法。那么一開始在QRAM 中準(zhǔn)備R， Σ的數(shù)量時將達(dá)到O(n2)。 盡管在理論分析中因?yàn)閚與資產(chǎn)數(shù)量N無關(guān)，所以被認(rèn)為是常數(shù)項(xiàng)，但在試驗(yàn)中這個常數(shù)值將達(dá)到1 04，是不可忽視的。如何在真實(shí)量子計算機(jī)上完成DME 算法也是挑戰(zhàn)之一。