求解整數線性規劃問題的量子近似優化算法

戚 晗，何婉瑩，邱 濤，Abdullah Gani

（1. 沈陽航空航天大學 計算機學院，沈陽 110136；2. 馬來西亞沙巴大學 計算機與信息學院，亞庇 87000）

近年來，物聯網、移動邊緣計算、車載自組織網絡（vehicular ad hoc network，VANET）發展迅速，隨之而來的一些基站和服務器選址、資源分配調度等問題受到廣泛關注。這些問題可以理解為在有限的可供選擇的方案中，尋找滿足一定約束的最好方案，因此類似此類問題可以歸結為整數線性規劃問題［1］。整數線性規劃問題屬于線性規劃中的一種，其中又分為二元整數線性規劃、混合整數線性規劃等。隨著問題規模的擴大，在約束條件的限制下尋找到整數線性規劃問題最優解的時間復雜度較高。

隨著“噪聲中尺度量子”時代［2］的快速發展，量子優化算法逐漸受到關注，量子近似優化算法（quantum approximate optimization algo‐rithm， QAOA）在2014 年由Farhi 等提出［3］，被認為是在近年可以實現量子霸權的算法之一［4］。Willsch 等［5］、Herrman 等［6］以最大切割問題和2-SAT 問題為例分析了QAOA 的性能，發現在D-Wave量子退火計算機上的性能優于量子計算機模擬器。Bengtsson 等［7］使用超導量子處理器實現QAOA，可以提高找到精確覆蓋問題最優解的成功概率。QAOA 變分態的對稱性會使算法本身有一定的局限性，但是基于伊辛模型的QAOA 優于原始QAOA，且優于Williamson 算法［8］。目前，QAOA 被用于解決許多組合優化問題，使用QAOA 找到問題全局最優解的概率比較依賴于量子線路的迭代次數（即步長P）。然而，當QAOA 的步長P 很低時，找到近似最優解的概率很小。Azad 等［9］在使用QAOA 解決車輛路徑規劃問題時，迭代24次找到最優解的概率約為30%，當迭代次數更高時，概率并沒有達到較高的水平。Zhang等［10］在使用QAOA 解決最小頂點覆蓋問題時，迭代2 次時，概率低于20%；迭代10 次時，概率可以達到約80%。Vikstl 等［11］在使用QAOA 解決精確覆蓋問題時，在迭代2 次時的單次測量下，找到最優解的概率為8.97%，需要重復多次測量才能提高概率。此外，QAOA 還被應用于哈密頓回路問題［12］和最大權獨立集問題［13］。盡管成功概率因不同問題而異，但QAOA 無疑是解決組合優化問題的一種較好的方法。2020 年，Choi等［13］針對無線自組織網絡中簇頭選擇的問題，采用量子經典混合方法求解簇頭選擇策略。為了讓所選擇的簇頭具有較高的能效，將此問題定義為最大加權獨立集問題，運用QAOA 求解最佳策略。2022 年，Fan 等［14］使用QAOA 解決圖計算中最短路徑問題，結果顯示，QAOA 在參數選擇和步長選擇上優于經典算法和基于Grover 的最短路徑算法，且需要的量子比特數更少。2023 年，Soloviev 等［15］為了解決許多傳統方法無法解決的貝葉斯網絡結構學習問題，采用基于3n（n-1）/2 個量子比特的QAOA 進行求解（其中n 為貝葉斯網絡中的節點數量）。結果顯示，在任何數據集上的性能都要優于經典和其他量子方法。在縮短整數線性規劃問題的求解時間方面，QAOA 是一種有意義且可行的解決方法。

本文針對整數線性規劃問題，提出一種基于改進目標哈密頓量的QAOA 的量子線路方法，提高在低迭代時找到最優解的概率并縮短程序執行時間。通過構造與整數線性規劃問題對應的數學模型，將經典伊辛模型量子化，得到本文所要解決問題的目標哈密頓量。推導哈密頓量對應的量子門線路，并減少線路中使用的量子門數量。本文使用本源量子的py‐Qpanda框架進行模擬，驗證了改進目標哈密頓量對找到基態概率和執行時間的影響。

1 資源分配問題的整數線性規劃公式

為了便于理解這一類場景下的整數線性規劃問題，本文以車載自組織網絡為例進行說明。通過車載自組織網絡實現車輛之間通信可以有效地避免出現交通擁堵和事故，因為司機可以獲得超出可視范圍內的實時路況信息。但是遠距離車輛無法直接通信，需要合作通信來獲取信息，合作通信需要使用中繼車輛節點的資源，但每輛車的剩余資源、移動速度和位置都是不固定的。在這種情況下，需要考慮中繼節點如何選擇。如果中繼節點集中在一輛車上，那么會造成這輛車過載，從而導致通信失敗，所以中繼節點的負載平衡十分重要，需要合理分配中繼節點的資源。針對資源分配下的整數線性規劃問題，本文以源節點代表需要請求計算資源的設備，以目標節點代表提供計算資源的設備。

假設目標節點具有相同的資源，源節點的資源需求量不同，那么影響延遲的主要因素即為通信距離。但是當目標節點的資源不足以處理所有請求時，會產生排隊延遲?，F在用G表示整個網絡，V 表示所有源節點的集合，R 表示所有目標節點的集合，d 表示源節點和目標節點之間的通信距離。資源分配的目的是實現負載均衡的同時盡量減少通信延遲。每個源節點的資源需求表示為vl，xvr為二元決策變量。如果源節點與目標可以相連，則xvr為1，否則為0。此問題可以表述為

可選取的目標節點不確定是否每個都可以與全部源節點連接，但是需要確?？傆幸粋€目標節點可以為其提供資源。通信條件相同的情況下，訪問距離越短，訪問延遲越小，且距離越短，可以連接的概率也越大，所以將源節點與目標節點之間的距離表示為dvr， xvr同樣為二元決策變量。如果源節點與目標連接，xvr為1，否則為0，此問題可以表述為

約束條件可以表述為

式中：n 為源節點的數目。式（1）和（2）為本文中的目標，式（3）為約束條件，即一個源節點只能分配給一個目標，不能分配給多個目標，規劃范圍內的所有源節點都要與相對應的目標連接，即所有源節點的請求都會受到相對應的目標的處理。與源節點相連的目標應該處理其所有的請求服務；一個目標節點可管理多個源節點；不考慮源節點與源節點之間、目標與目標之間的鏈路。

2 基于伊辛模型的量子近似優化算法

對于同一個問題，當設置的目標哈密頓量不同時，結果是不同的。在解決整數線性規劃問題時，設置了兩個不同的目標哈密頓量。一種是僅為問題函數（用Hp表示）求解目標哈密頓量；另一種是將問題函數和變形約束求和，并忽略公式中的常數，以獲得改進的目標哈密頓量（用Hc表示）。以上兩種方法用于設置不同的目標哈密頓量，為比較實驗做準備。

QAOA 的核心是通過從初始哈密頓量的基態，經過P 步迭代演化至目標哈密頓量的基態，在這個過程中需要使用經典計算來完成參數的更新。圖1 為QAOA 的算法結構圖，它的線路是以初始量子態為生成元的酉變換和以目標量子態為生成元的酉變換乘積的累積。Hb為初始哈密頓量，以Hb為生成元的酉變換等于e-iHbβi，Hp為目標哈密頓量，以Hp為生成元的酉變換為e-iHpγi。其中βi和γi所代表的是不同的變分參數。

圖1 QAOA的算法結構圖

對于目標哈密頓量的設定，采用伊辛模型來實現，本文所描述問題的目標哈密頓量可以表示為HA+HB［16］，其中

使用二元決策變量xr∈{0，1}代替自旋變量sr∈{-1，1}，即

2.1 設置目標哈密頓量HP

在資源分配整數線性規劃問題中，主要考慮兩個方面，即通信延遲和工作負載。本文將通信延遲轉化為通信距離，以HA為目標哈密頓量；將工作負載轉化為源節點的需求量，以HB為目標哈密頓量。將公式展開，得到問題中

讓A、B、C、D、E為

式中：C 和E 是常數，將自旋變量變為自旋泡利矩陣si→σzi，即得到最終目標哈密頓量

由于絕熱量子計算和量子門電路模型的計算能力相同［17］，本文將使用門電路構造QAOA 的線路，經過推導，最終得到基于伊辛模型的目標哈密頓量是19 個CNOT 門、CZ 門和RZ門的組合，如式（9）所示

2.2 改進目標哈密頓量Hc

根據資源分配問題的定義，本文將目標哈密頓量改進為原始目標函數與約束條件函數之和，可以表述為

其中A'、B'、C'、D'、E'為

式中：C'和E'為常數，它只與哈密頓量的相位有關，對本征態沒有影響［10］，所以可以進一步將其簡化為

將自旋變量變為自旋泡利矩陣si→σzi，即得到最終目標哈密頓量，如式（13）所示

經過推導，目標函數Hc可以表示為5 個CNOT門、RZ門組合，如式（14）所示

2.3 初始哈密頓量

設定初始哈密頓量為泡利X算符在每個量子位上的和，如式（15）所示

其基態是泡利算符對應特征能量的張量積，如式（16）所示

經過推導初始哈密頓量是H 門操作，原始QAOA和改進QAOA的初始哈密頓量均為Hb。

3 仿真實驗

本文使用基于python 的本源量子的pyQ‐panda 框架來執行QAOA，并解決整數線性規劃問題下的資源分配優化問題實例（4，2）。其中，4 表示源節點的數量，2 表示目標節點的數量。實驗所需的實際量子位N根據源節點和目標節點之間的通信鏈路確定。當源節點請求目標節點的資源時，中間距離對通信延遲和連接的可能性都有影響。如果距離超過其覆蓋區域，則無法利用目標節點的資源。

例如（4，2），使用式（17）來描述問題（它是歸一化數據），相應連接如圖2 所示，覆蓋區域分別表示E和F的可連接范圍。

圖2 問題示例（4，2）的連接圖

如圖2 所示，A 和F 之間的距離太大，超出了F 的通信覆蓋范圍，從而阻止了鏈路連接。實際距離應為∞，因此式（17）中使用1 表示無連接。為了編碼這個問題，將每個量子位分配給圖2中的每個鏈路，即如果使用了鏈路，則意味著對應的二進制決策變量為1，否則為0。為了使用QAOA，必須采取以下步驟：

（1）初始化線路，并加入初始哈密頓量基態所對應的量子門；

（2）通過初始化待優化的參數β和γ來確定量子門的初始旋轉角度，根據參數生成更新的量子線路，然后測量期望值；

（3）將當前的參數通過經典計算機再次優化會得到一組新的參數，通過新的參數計算新的損失值，直到滿足結束條件，即哈密頓量從初始基態演化到目標基態。

3.1 QAOA中參數更新方法對比

如何進行經典計算來更新參數至關重要，雖然也有使用張量網絡技術［18］來尋找參數的策略，但是最常用的方法是使用優化器計算。本文對比了AdaGrad 優化器［19］、Adam［20］、Mo‐mentum［21］、RMSProp［22］和Vanilla Gradient De‐scent［23］5種優化器分別在原始和改進QAOA 中的性能。通過對比這5 種優化器，可以看出哪一種更新參數的方法最適用于本文所提出的問題。本文使用損失函數的收斂值來評價不同優化器對算法的影響。

首先針對原始目標哈密頓量Hp進行實驗，不同優化器的性能如圖3 所示，在步長P=2 且優化器迭代次數為200時，5種優化器所形成的損失函數值均不斷減小，其中Momentum 優化器的整體損失函數值要低于其他優化器，但其前期波動較大，其余4 個優化器的損失函數收斂值相近?？梢钥闯龀薃daGrad 優化器之外，其余的損失函數收斂很快，在開始迭代到50 次左右就可以穩定收斂在固定數值，且保持在-0.900 3，但是前期同樣出現較小波動。AdaGrad 優化器雖然收斂不如其他優化器快，但是并沒有任何波動。綜合來說Momentum優化器雖然損失函數值更低，但是由于其波動較為嚴重，所以不適用于目標哈密頓量為Hp時的QAOA。

圖3 不同優化器性能

改進目標哈密頓量Hc在不同優化方法下損失函數的變化如圖4所示。在相同步長和迭代次數時，Hc的損失函數值均低于Hp，可以看到Adam、Momentum、RMSProp 和Vanilla Gra‐dient Descent 4 種優化器在迭代次數為50 時均可收斂，AdaGrad 優化器在迭代200 次時還沒有明顯收斂趨勢，但當迭代次數為500 時可以收斂。Momentum 優化器在收斂前一直存在明顯波動，這是因為引入歷史梯度信息動量。Adam 優化器是基于Momentum 優化器和RM‐SProp 優化器提出的，在迭代次數為20～30 也有小范圍波動，之后達到收斂狀態。RMSProp優化器和Vanilla Gradient Descent 優化器沒有波動，并且可以以較快的速度達到收斂值，其中Vanilla Gradient Descent 優化器的初始值要低于RMSProp 優化器。這與哈密頓量Hp存在一定差別，這是因為Hc與Hp本身存在差異，對于不同的目標哈密頓量，最合適的優化方法是不確定的。

圖4 Hc在不同優化方法下損失函數的變化

3.2 QAOA的不同目標哈密頓量性能對比

在實例（4，2）中，通信連接如圖2 所示，邊緣映射序列為｛A?E，B?E，A?F，C?E，C?F，D?E，D?F｝。一般來說，基于門電路的QAOA 使用量子門形成絕熱演化過程，使得哈密頓量從初始基態演化到目標基態。最后測量并得到哈密頓量的目標基態，在這個問題中，對應的基態是|11100101〉，映射序列是｛A?E，B?E，C?F，D?F｝。當P=2時，找到目標基態的概率可以達到54.156 3%，對于低迭代級別來說這是非常高的成功率。然而，隨著哈密頓量的改進，概率進一步增加到82.9%。概率分布如圖5 所示。

圖5 不同哈密頓量求得最優解的概率

從圖5 可以看出，目標哈密頓量為Hp時，所有解中概率較高的為“1011010”和“1100101”，雖然正確解“1100101”的概率高于“1011010”，但是極容易陷入局部最優解。而當目標哈密頓量為Hc時，正確解“1100101”的概率遠高于其他解，與目標哈密頓量為Hp時相比，更加容易跳出局部最優，以一個較高的概率得到問題的正確解。

QAOA 的執行時間與量子門的數量密切相關。與原始目標哈密頓量Hp相比，改進目標哈密頓量Hc減少了量子門的數量，這影響了找到最優解的時間。如圖6 所示，在同一設備上執行QAOA，隨著步長P 的增加，時間顯著增加，尤其是對于Hp?？梢钥闯?，盡管Hc的時間也隨著P 的增加而增加，但增加的幅度沒有Hp的大。在P=10 時，Hc的執行時間為Hp的19.16%，P=2 時Hc的執行時間為Hp的24.62%，Hc的平均執行時間為Hp的20.8%。這意味著當步長P 較大時，改進目標哈密頓量Hc的優勢更明顯。

圖6 不同哈密頓量的執行時間

3.3 時間復雜度分析

QAOA 是經典算法和量子算法的結合。經典部分通過經典優化器優化參數，其時間復雜度為O（poly（q）），其中q 是經典優化器迭代次數。本文中經典部分主要對比AdaGrad、Adam、Momentum、RMSProp 和Vanilla Gradi‐ent Descent 5 種優化器，其中Vanilla Gradient Descent每次更新需要計算整個數據集的梯度，且需要選擇合適的學習率，學習率太小會導致收斂緩慢，太大則會導致收斂效果波動。Mo‐mentum 動量法是在隨機梯度下降法（Stochas‐tic Gradient Descent， SGD）基礎上提出的，SGD 具有較快的下降速度，但因為方差頻繁更新，容易出現劇烈波動的情況，但也正是因為SGD 的波動性使其可能收斂到更好的局部最優。Momentum 動量法則可以抑制SGD 的搖擺情況，在減輕波動情況的同時，更容易收斂到局部最優。在本文實驗中Momentum優化器前期波動，最終收斂值更低。AdaGrad 在訓練過程中累計平方梯度越來越大，導致學習率過早減少，迭代時收斂緩慢。在本文實驗中AdaGrad 優化器在迭代200 次時沒有達到明顯收斂效果，速度緩慢。RMSProp 解決AdaGrad學習率急劇下降的問題，加速效果更好，收斂速度快。在本文實驗中，RMSProp 優化器收斂速度最快。Adam 是Momentum 動量法和RM‐SProp 的結合體，下降速度快，參數更新穩定，但是可能存在小波動。在本文實驗中Adam 優化器在迭代過程中出現小范圍波動，下降速度最快。

量子部分需要完成從初始狀態到目標狀態的演化，這部分的復雜度相當于量子態演化的時間，即O（poly（P）），其中P 是迭代次數（即步長）。QAOA 是從絕熱算法演變而來的，該系統可以通過使用絕熱算法來保證開始狀態是基態，并且最終狀態也是基態，但是成本相對較高（例如時間成本）。相比之下，QAOA 運行時間相對較短，但獲得基態的概率也受到限制。因此，如果增加迭代次數P，總體成功概率也會得到一定的提高。

根據上述分析，QAOA 的時間復雜度如式（18）所示

從式（18）可以看出，QAOA 在時間復雜度上的性能要優于其他經典算法，因為QAOA 的復雜性與所需要求解問題的大小無關，而經典算法的復雜性隨著問題規模的增加呈指數增長。

4 結論

本文利用QAOA 解決了與資源分配相關的整數線性規劃問題，將節點間的鏈路映射到量子位來編碼哈密頓量Hc、Hp和Hb。通過從原始的目標哈密頓量Hp轉換到改進的目標哈密頓量Hc，量子電路的運行時間和所需量子門的數量都顯著減少。此外，本文分析并比較了原始目標哈密頓量Hp和改進的目標哈密頓量Hc對找到近似最優解的影響。實驗結果證明了QAOA 在解決這些問題方面的有效性，即使有少量的迭代，該算法也有很高的成功率。例如，當目標哈密頓量是Hp時，成功率是54.156 3%，而當它是Hc時，概率增加到82.9%。這種低迭代水平轉化為實現算法所需的低電路深度，證實了其在近期量子機器上的可行性。由于QAOA 的時間復雜性不受問題大小的影響，因此它更適合大規模和多節點場景。此外，哈密頓量Hc的電路需要更少的量子門，這可以縮短執行時間，平均執行時間為Hp電路的20.8%。在未來，將會對QAOA 的內部結構進行改進，發現更有效的參數更新技術，并將不同的模型應用于不同的組合優化問題。