基于量子判別分析法的量子連續投資組合優化算法

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

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

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

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

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

1 馬科維茨均值方差模型

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

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

對于不設置收益目標，只追求最大利潤風險比的投資者，該模型可以轉化為：

也可以轉化為：

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

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

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

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

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

2 量子線性判別分析

2.1 線性判別分析模型

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

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

在此將特征矩陣V轉化為特征向量w分析這個問題。上述問題的對偶問題可以變為：

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

將上式對w求偏導后，得到KKT 條件：

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

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

問題轉變為求解厄米特矩陣最大特征值所對應的特征向量：

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

2.2 厄米特鏈積

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

圖1 HCP 線路圖

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

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

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

2.3 密度矩陣指數化算法

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

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

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

圖2 DME 算法線路圖

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

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

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

相位估計算法使用控制下的哈密頓模擬，DME 算法也證明出，將交換門S變為控制交換門就可以實現控制哈密頓模擬的操作[39]。

3 基于QLDA 的量子投資組合優化

3.1 模型規約

首先將式(7)按照式(2)中的均值方差構造，其模型變為：

因為該模型將解w編碼為量子態的概率幅，受量子計算限制，其模平方和為1。為方便書寫，記R=E(r)E(r)T。 這里E(r)表 示凈利率。然后將R，Σ和R1/2Σ-1R1/2譜分解得到：

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

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

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

進而有如下引理。

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

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

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

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

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

由于目標函數的縮放性質，解乘以一個正常數時，取值不變。因此：

之后，記：

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

3.2 基于QLDA 的投資組合優化算法

結合第2 部分內容，使用QLDA 計算投資組合問題。

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

2）使用DME 算法將密度算子轉化為關于R與Σ的控制哈密頓算子，以此構造HCP 線路。

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

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

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

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

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

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

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

3.3 性能分析

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

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

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

表1 不同投資組合優化算法的對比

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

4 結 束 語

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

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