一類保序最優化問題的迭代算法

鄧文麗，朱瑩瑩

（江西師范大學 數學與信息科學學院，南昌 330027）

0 引言

約束條件下的統計推斷已經成為統計分析的一個重要領域，保序回歸是約束條件下的統計推斷的一類典型問題。保序回歸的研究從20世紀中葉就已經開始，到20世紀60年代，保序回歸得到了廣泛的重視。Bartholoremew,Barlow等對此進行了研究，并且于1972年聯合出版了保序回歸的第一本專著 《The theory and application of Isotonic Regression》，使得保序回歸問題得到了進一步的重視和發展，成為熱門話題。

在回歸分析問題中，假定自變量為 z=（z1，z2，…，zn），因變量為 y=（y1，y2，…，yn）。 在很多情況下，研究者并沒有關于回歸函數的很多信息，而只是一些關于變量的約束條件，比較典型的一類就是，假設z的元素是經過排序的，z1≤z2≤…≤zn，yi會隨著zi遞增，這樣一種情形被稱為保序回歸；yi會隨著zi遞減，被稱為逆序回歸。兩種情況合在一起，被稱為單調回歸。詳細介紹可參見文獻 de Leeuw（2005）。

令 x=（x1，x2，…，xn）表示 y 所對應的回歸擬合值，保序回歸可以表述為：

尋找 x*=（x1*，z2*，…，zn*）∈G，G={x∈Rn|x1≤x2≤…≤xn}使得

其中 ω=（ω1，ω2，…，ωn），ωi≥0 （i=1，2，…，k）是給定的權數，ω1+ω2+…+ωn=1。這時x*被稱為y的保序回歸。

保序回歸問題的求解過程中，有一個最基本的定理是：如果yi＞yi+1，那么它們對應的擬合值滿足x^i=x^i+1。因此，保序回歸的求解歸結為找到指標集{1，2，…，n}的一個合理劃分，B1，B2，…Bk，k≤n，使得 B1+B2+…+Bk={1，2，…，n}，對劃分中的任意一塊Bj而言，任給i∈Bj，xi都是相等的，在滿足這樣條件的劃分中，尋找（1）的解。 Ayer et al.（1995）提出的 PAVA（Pool Adjecent Violators Algorithm）算法使問題（1）從理論上得到了很好的解決。PAVA算法求得的“最優”劃分B1，B2，…Bk，k≤n，對劃分中的任意一塊 Bj而言，任給 i∈Jj，。

本文將經典的保序回歸問題拓展到絕對損失下，去解決下面的最優化問題：

尋找 x*∈G，G={x∈Rn|x1≤x2≤…≤xn}使得

其中 ω=（ω1，ω2，…，ωn），ωi≥0 （i=1，2，…，k）是給定的權數，ω1+ω2+…+ωn=1。

本文將根據PAVA思想過程給出一種迭代算法，通過迭代過程分析得出最優化問題（2）的解，并給予證明。

1 利用迭代算法求解最優化問題

J=（B1，B2，…，Bk），B1，B2，…，Bk，k≤n 是指標集{1，2，…，n}的一個劃分。 若 B={p，p+1，…，q}，1≤p≤q≤n，J中包含 p-1的元素，用B-表示；包含q+1的元素用B+表示。若p=1時，B-=φ，若 q=n 時，則 B+=φ。 用［AM（B）L，AM（B）U］記指標集 B所對應y的元素的加權中位數所在集合，有

類似PAVA算法的思想，下面針對本文討論的最優化問題（2）給出一種迭代算法。

G 表示可行域{x1，x2，…，xn|x1≤x2≤…≤xn}，如果（y1，y2，…，yn）∈G，則=yi，i=1，2，…，n 否則，令 J={{1}，{2}，…，{n}}，B={1}，進入迭代算法。

（1）若 B+=φ，則迭代結束，xi*∈［AM（B）L，AM（B）U］，i∈B，B∈J；否則將 AM（B）L，AM（B）U和 AM（B+）L，AM（B+）U進行比較。

若 AM（B）U≤AM（B+）L，B=B，回到初始；

若 AM（B）U＞AM（B+）L，但 AM（B）L≤AM（B+）U，即

若 AM（B）L＞AM（B+）U進入第（2）步。

（2）J=J{B，B+}∪（B∪B+），令 B=B∪B+，計算AM（B）L=AM（B）U。

（3）若 B-=φ，進入第（1）步；否則進行下面的判斷：

若 AM（B-）U≤AM（B）L，回到第（1）步；

若 AM（B-）U＞AM（B）L，但 AM（B-）L≤AM（B）U，即

進入第（1）步。

AM（B-）L＞AM（B）U，進入第（4）步；

（4）J=J{B，B-}∪（B∪B-），B=B∪B+，進入第（3）步。

迭代的結果，指標集{1，2，…，n}被劃分成了 B1，B2，…，Bk，k≤n 是指標集的一個劃分。 B1+B2+…+Bk={1，2，…，n}，而且最優解屬于集合

下面證明（5）所表示的就是最優化問題（2）的解集。

2 最優解的相關證明

定理 1 假定指標集 B=（p，p+1，…，q），?i∈B，yi∈R，yi所對應的權數記為 wi，yp≥yp+1≥…≥yq。 數列{yi，i∈B}的加權中位數所在的集合表示為［AM（B）L，AM（B）U］，［AM（B）L，AM（B）U］如（3）、（4）所示，那么

①使得

定理 2 若 B=（p，p+1， …，q），1≤p≤q≤n 任給 p≤i≤+1≥yi+2≥…≥yq，記

①若 T≠φ，選取

可使得

證明：①式顯然成立，下面證②。

根據定理1中的結論①，又因為yp≥yp+1≥…≥yi，yi+1≥yi+2≥…≥yq， 那么使得成立的一定滿足使得成立的一定滿足…=。 所以使得（7）式成立的滿足

因為 T=φ，即 AM（Ui（B））U＜AM（Li（B））L，顯然使得（7）式成立的{，i∈B}還應滿足

定理 3 若 B={p，p+1，…q}，1≤p≤q≤n，任給 p≤i≤q，Li（B）{p，p+1，…，i}。 記記

①若T≠φ，選取

可使得

例：已知 y1=4，y2=5，y3=1，y4=6，y5=8，y6=7 并且 wi=1/6，i=1，…，6。 求 x* 使在滿足 x1≤x2≤…≤x6的限制下達到最小。

解：y1=4＜y2=5，y2=5＞y3=1，所以 B1={1}，B2={2，3}

取 AM（B1）L=AM（B1）U=4，AM（B2）L=AM（B2）U=5，由此重新記 AM（B1）L=AM（B1）U=4，AM（B2）L=4，AM（B2）U=5；AM（B2）U＜y4＜y5；而 y5＞y6，所以 B3={4}，AM（B3）L=AM（B4）U=6；B4={5，6}，AM（B4）L=AM（B4）U=8，顯然 AM（B4）L＞AM（B3）U=6。 所以選取，J={{1}，{2，3}，{4}，{5，6}}［4，5］，∈［7，8］，

用上述迭代計算該模型的解的次數為O（n），當n不是很大時，其運算并不復雜，當n很大時，利用計算機輔助仍然可以比較快得出結果。因此，加權最小絕對偏差在保序限制下迭代計算仍是比較實用的，有奇異值出現時，其所得解也比加權最小平方損失下的解更穩健。

[1]Michael J BEST，Nilotpal Chakravarti.Active Set Algorithms for Isotonic Regression;A Unifying Framework[J].Mathematical Programming，1990，（47）.

[2]De Leeuw，Jan，Hornik，Kurt，Mair，Patrick.Isotone Optimization in R:Pool-Adjacent-Violators Algorithm （PAVA） and Active Set Methods[M].UC Los Angeles:Department of Statistics，2009.

[3]De Leeuw.Monotonic Regression，Encyclopedia of Statistic in Behavioral Science[M].New York:Wiely，2005.

[4]史寧中.保序回歸與最大似然估計[J].應用概率統計，1993，（19）.

[5]刑務強.保序回歸的研究及應用[D].西北工業大學，2002.

[6]趙選民，刑務強.凸規劃下的保序回歸[J].數理統計與管理，2003，（22）.