基于SOP-SeedBS的高維數(shù)據(jù)稀疏變點檢驗

銀炳皓，施三支，齊德全

（長春理工大學 數(shù)學與統(tǒng)計學院，長春 130022）

變點估計和檢驗問題是統(tǒng)計學中的一個經(jīng)典問題，用于識別時間序列中的突變，發(fā)生這種突然變化的時間點稱為“變點”。通過估計變點的位置，可以將原始數(shù)據(jù)集劃分為較短的段，然后使用平穩(wěn)時間序列的方法進行分析，解決各個領(lǐng)域相關(guān)的問題。變點問題在金融行業(yè)[1]、互聯(lián)網(wǎng)安全[2]、生物醫(yī)學[3]等領(lǐng)域都有著廣泛的應(yīng)用，例如，變點估計可用于發(fā)出有關(guān)異常金融事件的警報，檢測網(wǎng)絡(luò)上的分布式拒絕服務(wù)攻擊，精確定位人腦中腦波活動的開始等。隨著時代的發(fā)展和科技的進步，現(xiàn)實生活中的數(shù)據(jù)多以高維的形式出現(xiàn)，高維數(shù)據(jù)流導致了數(shù)據(jù)量的激增，這就要求變點分析方法能有所進步，并且在運算時間和方法準確性上達到一個良好的平衡效果。

經(jīng)典的變點檢驗關(guān)注的是單變量時間序列，例如2012年Killick[4]提出的PELT方法是通過找到變點可能的數(shù)量和位置的成本函數(shù)最小值，來推斷變點的最佳數(shù)量和位置。2014年Frick和Munk等人[5]提出的SMUCE算法針對指數(shù)族回歸中的變點問題，引入了一種新的估計量；2014年Fryzlewicz[6]提出的 WBS 算法和 2018 年 Baranowski等人[7]提出的最窄超閾值法（NOT）提高了二元分割算法的性能。然而，經(jīng)典的單變量變點檢驗方法通常不直接適用于現(xiàn)代應(yīng)用中經(jīng)常遇到的高維數(shù)據(jù)集。因此，學者們提出了幾種方法來估計和檢驗高維數(shù)據(jù)中的變點，并且為了解決維數(shù)災難的問題，現(xiàn)有的幾種高維變點檢驗方法假設(shè)了變化信號具有某種形式的稀疏性，稀疏性假設(shè)大大降低了原始高維問題的復雜性。例如，在 2015年 Jirak[8]使用 CUSUM 統(tǒng)計量對不同子部分進行L∞聚合，對于稀疏變點的效果很好；2015 年 Cho 和 Fryzlowicz[9]的稀疏二元分割算法，也考慮了稀疏性；2018年Wang和Samworth[10]基于投影的方法中假設(shè)了一個變點前后的均值差異只在小部分維數(shù)坐標中是非零的；2019 年 Enikeeva和 Harchaoui[11]提出了基于具有稀疏性的掃描統(tǒng)計量的方法。關(guān)于高維數(shù)據(jù)的變點檢驗問題的工具仍在不斷創(chuàng)新和發(fā)展，但大多不具有普適性，且運算時間較長。

本文通過研究具有稀疏變點的高維數(shù)據(jù)的特點，提出了一種基于稀疏最優(yōu)投影的多變點檢驗方法，SOP-SeedBS方法能夠檢驗出數(shù)據(jù)的維度p在小于或大于等于數(shù)據(jù)量n的情況下變點的位置及個數(shù)，并研究了變點間隔大小對該檢驗方法穩(wěn)定性的影響。通過仿真模擬，在多種情況下與Wang和Samworth的Inspect方法（2018）進行了比較，結(jié)果表明，SOP-SeedBS方法具有較高的準確度，并在計算時間方面具有明顯的優(yōu)勢。

1 SOP-SeedBS高維稀疏多變點檢驗

1.1 稀疏變點模型

設(shè)X1,…,Xn是n個相互獨立的p維隨機向量，Xt～Np(μt,σ2Ip)，1 ≤t≤n，σ2是方差，記

其中，p可能與序列的長度n相等，甚至大于序列的長度n。

假設(shè)存在υ個變點，變點檢驗的問題能被表示為經(jīng)典的假設(shè)檢驗問題

假設(shè)備擇假設(shè)成立，數(shù)據(jù)可劃分成υ+1段的分段常數(shù)均值結(jié)構(gòu)，換句話說，如果有υ個變點 1≤z1＜z2＜…＜zυ，即?0 ≤i≤υ。記z0=0,zυ+1=n， 記θ(i)=μ(i)-μ(i-1),i=1,…,υ。

本文假設(shè)均值的變化是稀疏的，即一個變點前后的均值差異只在小部分維度（p）中是非零的。記‖θ(i)‖0表示向量θ(i)中非零元素的個數(shù)，則‖θ(i)‖0≤k，k∈{1,…,p}且k?p。本文只考慮參數(shù)μ(i)的檢驗問題，其他參數(shù)皆視為討厭參數(shù)。假設(shè)變點位置z1，z2，…，zυ滿足：

其中，τ是變點的最小間隔，0＜τ＜1。均值變化的大小滿足：

其中，?是均值變化大小的閾值。

1.2 稀疏最優(yōu)投影（SOP）

Wang等人[10]（2018）通過對高維數(shù)據(jù)進行高維CUSUM變換進而推導出數(shù)據(jù)矩陣的最優(yōu)投影方向，并解決了計算第一個左奇異向量引出的凸優(yōu)化問題。因此對于高維稀疏變點數(shù)據(jù)，本文根據(jù)類似的思想，并基于奇異值分解和凸松弛優(yōu)化方法將高維數(shù)據(jù)投影成一維數(shù)據(jù)，由于本文研究的高維數(shù)據(jù)中均值變化是稀疏的，從而稱該方法為稀疏最優(yōu)投影（Sparse Optimal Projection）。

對于任意的a∈Rp，滿足 ‖a‖2=1，其中

選擇a=θ/‖θ‖2來最大限度地擴大兩個數(shù)據(jù)片段之間的均值差異的大小。然而θ在實踐中通常是未知的，本文的想法是尋找一個最優(yōu)的投影方向v=θ/‖θ‖2，文獻[10]啟發(fā)了本文對數(shù)據(jù)矩陣X進行分解，再利用奇異值分解得到一個合適的投影方向后，將數(shù)據(jù)投影到一維空間。

假設(shè)均值變化只出現(xiàn)在數(shù)據(jù)矩陣X的k個列中，且在沒有發(fā)生均值變化的p-k個列中的元素均值為μt=0，t∈{1,…,n}，k∈{1,…,p} 。對數(shù)據(jù)矩陣X進行分解，得到X=μ+W，其中μ=(μ1,…,μn)T∈Rn×p是均值矩陣，W=(W1,…,Wn)T∈Rn×p，W1,…,Wn是獨立的Np(0,σ2Ip)的隨機噪聲。當σ已知，零假設(shè)成立時，μ為恒定的常數(shù)向量，備擇假設(shè)下只存在一個變點z時，矩陣X的元素均值可以計算出：

其中，Xr,j為矩陣X中第r行，第j列的元素；t∈{1,…,n} ；j∈ {1,…,p}。記前z行向量元素相同，后n-z行向量元素相同。

這里St,j為第j列中變點z前后分段數(shù)據(jù)之和，t∈{1,…,n}，j∈{1,…,p} ，即：

則有：

文獻[10]中，計算矩陣XT的第一個左奇異向量是一個組合優(yōu)化問題。公式（10）是一個計算復雜度很高的凸優(yōu)化問題，解決這個問題的一種方法是使用優(yōu)化問題中的凸松弛得到M?來代替X，具體過程如下：

懲罰函數(shù)可以將有約束的目標函數(shù)轉(zhuǎn)化為無約束的目標函數(shù)，即當?shù)玫降腗不滿足約束條件時，估計的?可能會遠離最優(yōu)解。通過添加一個懲罰項 -λ‖M‖1，優(yōu)化問題變?yōu)橐韵滦问剑?/p>

其中，λ為軟閾值。則有：

記是的第一個左奇異向量，即：

將高維數(shù)據(jù)X沿著最優(yōu)方向投影，得到一維數(shù)據(jù)，記為Y。

如果數(shù)據(jù)存在變點，投影后的一維數(shù)據(jù)呈現(xiàn)出了分段常數(shù)的特征，適用于此種數(shù)據(jù)的幾種經(jīng)典的單變量多變點檢驗方法，包括2012年Killick[4]給 出 的 PELT 方 法 、2014 年 Fryzlewicz[6]提出的 WBS方法、2018年 Baranowski等人[7]給出的最窄超閾值方法（NOT）以及2020年Kovács等人[12]提出的SeedBS方法。PELT方法是通過找到變點可能的數(shù)量和位置的成本函數(shù)最小值來推斷變點的最佳數(shù)量和位置，計算復雜度較高。WBS方法和NOT方法都是基于二元分割思想的改進算法，NOT方法中隨機抽取間隔的想法與WBS方法類似，但NOT方法更關(guān)注數(shù)據(jù)中可能存在變點的最小間隔。WBS和NOT方法中用于搜索變點的間隔都是通過隨機抽樣構(gòu)造的，在數(shù)據(jù)量較大的情況下，需要構(gòu)造的隨機間隔數(shù)量會很多，且會出現(xiàn)許多長度與原始序列相同或者重復出現(xiàn)的間隔，這增加了算法復雜度和運算時間，特別是高維稀疏變點數(shù)據(jù)中變點較少的情況下。2020 年 Kovács等人[12]提出的種子二元分割（SeedBS）是一種解決大規(guī)模數(shù)據(jù)變點檢驗問題的方法，該方法能夠有效地降低計算復雜度，提高計算效率。因此，本文選擇SeedBS算法來對投影后的一維數(shù)據(jù)進行多變點檢驗。

1.3 種子二元分割（SeedBS）

種子二元分割方法是通過構(gòu)造一個具有確定性，且可以預先計算的搜索區(qū)間集來搜索變點。區(qū)間集的構(gòu)造獨立于潛在的變點數(shù)量，并能達到搜索間隔的總長度盡可能小的效果。

其中，μt是t時刻的信號；σ是恒定標準差；εt是正態(tài)噪聲，記有υ個變點1≤z1＜z2＜ ...＜zυ，則μt可以被分為υ+1段，其中第i段中。

記a∈(1/2,1]為一個給定的衰減參數(shù)，對于任意的q：1≤q≤log1/a(n)，將第q層定義為由nq個區(qū)間組成，初始長度為lq且進行確定性位移sq的集合Lq：

例如當a=1/2時，得到一個重疊的二元分隔區(qū)間。來自前一層的每個間隔都被分成左、右和重疊的中間間隔，每個間隔的大小都是前一層的一半。當間隔的長度在層上迅速衰減時，隨著間隔長度變小，間隔的數(shù)量迅速增加，還保證了在同一層內(nèi)的連續(xù)間隔之間有足夠大的重疊。生成的間隔如圖1所示[12]。

圖1 a=1/2時生成搜索間隔的可視化圖

用于種子二元分割的檢驗統(tǒng)計量是定義一段內(nèi)分割點在s的CUSUM統(tǒng)計量：

其中，0≤l＜r≤n和h=r-l為整數(shù)，達到CUSUM統(tǒng)計量絕對值最大值的位置s記為最好的分割點。因此，計算m個種子間隔的檢驗統(tǒng)計量，得到一個候選變點集，記為，如果候選變點集合中對應(yīng)間隔的檢驗統(tǒng)計量超過了閾值ζn，就認為得到一個最終的變點估計集合O(0,n]，并將聲明為一個變點：

其中：

閾值ζn可以通過強化施瓦茨信息準則計算：

其中，α≥ 1；是模型（2）分段參數(shù)μ(i)的最大似然估計；nυ表示需要估計參數(shù)的總個數(shù)，包括變點位置。另外，當α=1時相當于施瓦茨信息準則。

在具有固定的方差σ2和分段常數(shù)均值的單變量高斯信號模型中，也可以通過貝葉斯信息準則（BIC）來計算閾值：

其中，σ2是方差；n為樣本量；T={z1,…,zυ} ? {1,…,n} ；|T|為集合T的基數(shù)。通過BIC準則計算出的閾值篩選出的變點個數(shù)要多于強化施瓦茨信息準則篩選出的變點個數(shù)。

1.4 SOP-SeedBS方法

本文結(jié)合了類似 Wang等人[10]（2018）的投影思想和 Kovács等人[12]（2020）提出的種子二元分割（SeedBS）方法，進行高維數(shù)據(jù)的稀疏變點檢驗，先將高維稀疏變點數(shù)據(jù)沿著最優(yōu)的投影方向投影到一維，再利用SeedBS方法對一維數(shù)據(jù)進行多變點檢驗，SOP-SeedBS方法具體步驟如下：

輸入：數(shù)據(jù)矩陣X=(X1,...,Xn)T∈Rn×p，軟閾值λ，衰減參數(shù)足夠大的變點上界K，max閾值ζn。

（1）對數(shù)據(jù)矩陣X進行分解并推導出最優(yōu)投影方向。

（3）對進行SVD分解，得到第一個左奇異向量，計算投影數(shù)據(jù)

（4）根據(jù)衰減參數(shù)a構(gòu)造包含m個種子區(qū)間的區(qū)間集。

（5）計算m個種子間隔的檢驗統(tǒng)計量，并根據(jù)閾值ζn確定變點位置：

2 數(shù)值模擬

構(gòu)造由模型1生成的具有均值變點的高維數(shù)據(jù)時，設(shè)置參數(shù)σ2=1，數(shù)據(jù)量n分別取200，500，1 000，2 000，維數(shù)p分別為200，500，1 200，將高維數(shù)據(jù)噪聲表示為ε，分別取0.1，0.5，0.7，1，變點位置滿足z/n={0 .25,0.5,0.75}，定義?(i)為第i個變點的均值跳躍大小，=(?(1),?(2),?(3))=(0.4,0.8,1.2)，?(0)=0，，其中k控制著發(fā)生變化維度數(shù)量的稀疏程度，取值分別為10，40。根據(jù)文獻[10]的建議，。根據(jù)文獻[12]的建議，設(shè)置衰減參數(shù)a= 2。另外，下面三種情形分別考慮了均值結(jié)構(gòu)中發(fā)生變化的維數(shù)坐標(p)上的重疊問題，高維數(shù)據(jù)中均值結(jié)構(gòu)發(fā)生變化的三種不同重疊情形說明如下：

情形1（S1）：在沒有重疊的情況下，對于相應(yīng)的變點，發(fā)生變化的k個維數(shù)坐標是來自完全不相交的坐標集；

情形3（S3）：在完全重疊的情況下，則對于相應(yīng)的變點，變化的維數(shù)發(fā)生在相同的k個坐標上。

本文先給出了SOP-SeedBS方法在情形3下，參數(shù)設(shè)置為σ2=1，n=2 000，p=200，k=40，??=(?(1),?(2),?(3))=(0.4,0.8,1.2)，變點位置z/n={0.25,0.5,0.75}時，不同噪聲下的變點估計結(jié)果可視化，用ε表示高維數(shù)據(jù)的噪聲，不同ε取值對于投影后一維數(shù)據(jù)的變點檢驗結(jié)果如圖2所示，圖2的結(jié)果顯示了當ε的取值越小時，噪聲越小，分段特征越明顯。另外，當均值跳躍越大，分段特征也越明顯。后續(xù)的仿真實驗中，噪聲ε皆設(shè)置為1。

圖2 不同 ε下的投影后一維數(shù)據(jù)的變點檢驗結(jié)果

表1給出了σ2=1時，SOP-SeedBS方法和Inspect方法的比較結(jié)果，兩種方法在四組不同的參數(shù)設(shè)置M1、M2、M3、M4下，且每組參數(shù)設(shè)置都具有三種不同的情形，通過100次運行算法（兩種方法通過相同參數(shù)設(shè)置下生成的高維數(shù)據(jù)進行比較）檢測到的變點數(shù)量的頻數(shù)和平均Hausdorff距離dH分別來度量估計變點的數(shù)量和估計變點位置的準確度。M1的參數(shù)設(shè)置為n=2 000，p=200，k=10，=(0.8,1.6,2.4)；M2的參數(shù)設(shè)置為n=2 000，p=200，k=40，=(0.4,0.8,1.2)；M3的參數(shù)設(shè)置為n=500，p=500，k=40，=(0.4,0.8,1.2)；M4的參數(shù)設(shè)置為n=200，p=500，k=40，=(0.4,0.8,1.2)。表中Z1，Z2，Z3，Z4，Z5表示估計到的變點個數(shù)，取值分別為1，2，3，4，5。

表1 兩種方法在不同參數(shù)設(shè)置下進行100次模擬的估計變點頻數(shù)和平均Hausdorff距離的比較

從表1可以看出，在σ2=1，n=2 000，p=200，k=10的參數(shù)設(shè)置下，SOP-SeedBS方法在三種不同均值變化重疊情況下的結(jié)果都優(yōu)于Inspect方法，進行100次模擬的平均Hausdorff距離分別為 6.87，2.32，0.47，而 Inspect方法分別為12.72，26.01，13.91；第二組實驗將第一組實驗中的發(fā)生變化的維數(shù)k=10替換為k=40，第一個變點的均值跳躍大小由0.8變?yōu)?.4，SOPSeedBS方法在三種情形下進行100次模擬得到的平均Hausdorff距離分別為4.48，2.42，0.34。雖然兩組實驗在變點頻數(shù)上同樣具有良好的檢驗結(jié)果，但第二組結(jié)果與第一組相比，在均值跳躍變小的情況下，第二組結(jié)果中情形1和情形3的平均Hausdorff距離要比第一組更小，說明k的大小也會對實驗結(jié)果產(chǎn)生影響。

從SOP-SeedBS方法的第三組和第四組實驗結(jié)果可以看出，在其他參數(shù)固定的情況下，n越大，對比實驗結(jié)果顯示n較大時檢驗結(jié)果較好，說明由于樣本量的增加，適當增加種子區(qū)間的個數(shù)有利于進行變點檢驗。本文還通過模擬實驗發(fā)現(xiàn)，當其他參數(shù)固定時，維度p變大，對比實驗的結(jié)果表現(xiàn)出檢驗效果相近，因此維度的增加對于變點檢驗效果影響不大。

圖3給出了SOP-SeedBS方法與Inspect方法在其他參數(shù)不變的情況下，不同的n和p設(shè)置下算法平均運算時間的比較，可以看出SOP-SeedBS方法在算法運算速度上具有較大的優(yōu)勢。圖3（a）給出了當p=200，n分別為 500、1 000、2 000時，兩種方法100次模擬的平均運算時間；圖3（b）給出了當p=500，n分別為 200、500、1 000時，兩種方法進行100次模擬的平均運算時間。

圖3 SOP-SeedBS方法與Inspect方法計算時間比較

圖3中的結(jié)果顯示，當p=200，n=2 000時，Inspect方法進行100次模擬的平均運算時間達到了13.7 s，而SOP-SeedBS方法的平均運算時間僅為0.14 s，約是Inspect方法運算時間的1%。當p=500，n=500，Inspect方法進行 100次模擬的平均運算時間是13.6 s，而SOP-SeedBS方法的平均運算時間僅為0.15 s，約是Inspect方法運算時間的1%，由此能體現(xiàn)SOP-SeedBS方法在運算時間上優(yōu)于Inspect方法。

為了研究變點間隔對實驗結(jié)果準確性的影響，從前面的對比實驗表1中發(fā)現(xiàn)，當均值變化為情形1時，SOP-SeedBS方法的效果不如情形2和情形3的情況，這也符合不同變點在不同稀疏維數(shù)坐標中發(fā)生變化時，信號強度會減弱的預期效果。因此，表2顯示的是在情形1下的模擬結(jié)果，即發(fā)生變化的k個維數(shù)坐標是來自完全不相交的坐標集。表中Z1、Z2、Z3、Z4表示估計到的變點個數(shù)，取值分別為1，2，3，4。

表2 不同變點間隔在p=200，σ2=1時四組不同參數(shù)設(shè)置下100次模擬的結(jié)果

表2中有四組不同的參數(shù)設(shè)置B1,B2,B3,B4，其中B1的參數(shù)設(shè)置為zi+1-zi=20，表示相鄰變點之間間隔20個數(shù)據(jù)，信號的均值跳躍大小為=(0.8,1.6,2.4)；B2的參數(shù)設(shè)置為zi+1-zi=30,=(0.8,1.6,2.4)；B3參數(shù)設(shè)置為zi+1-zi=40,=(0.8,1.6,2.4)；B4的參數(shù)設(shè)置為zi+1-zi=50,=(0.8,1.6,2.4)。從表2的結(jié)果可以看出，通過模型1生成的四組不同參數(shù)設(shè)置下的變點數(shù)據(jù)，在σ2=1，p=200，n=2 000，k=40，變點間隔為20時，真實變點個數(shù)為3，三個變點的均值跳躍大小分別為 0.6、1.2、1.8，SOP-SeedBS方法在100次模擬中估計正確變點個數(shù)次數(shù)僅為 31 次，而在σ2=1、p=200、n=2 000、k=60時，100次模擬中估計正確變點個數(shù)次數(shù)僅有61次。由最后一組實驗可以看出當變點之間的間隔達到50個數(shù)據(jù)點及以上時，SOP-SeedBS方法不僅能穩(wěn)定地檢測到正確的變點個數(shù)，在準確性上也能達到很好的效果。

3 實證分析

本節(jié)研究了R中ecp包的膀胱腫瘤微陣列數(shù)據(jù)集（Bladder Tumor Micro-Array Data），目的是檢驗樣本目標區(qū)域染色體拷貝數(shù)的變化，常常用于醫(yī)學中檢驗癌癥的研究，染色體拷貝數(shù)變化反映了DNA的擴增和驟減。從醫(yī)學的角度來說，擴增的片段可能存在致癌基因，驟減的片段可能存在抑癌基因。本文對數(shù)據(jù)集中43個個體（p）的基因組上的2 215個不同基因位點（n）的對數(shù)強度比測量值進行檢驗，圖4繪制了前10個個體在前200個基因位點上的對數(shù)強度比散點圖，并采用SOP-SeedBS方法估計變點位置（其中虛線表示檢測到的變點位置），同時給出了文獻[10]和文獻[14]中Inspect方法和 Adapt-WBS方法的檢驗結(jié)果，來驗證SOP-SeedBS方法的有效性。

圖4 前10個個體在前200個位點上的ACGH數(shù)據(jù)散點圖

從圖4中可以看出，SOP-SeedBS方法檢測到的變點結(jié)果為：15，26，28，52（56），73，91，102，134，147（146），177（174），183（180），185（186），188（191），與文獻[10]和文獻[14]中檢驗到的變點相近，其他無法檢測到的變點，例如位點30到40之間的數(shù)據(jù)出現(xiàn)的分段特征不明顯，許多個體在位點110～130之間測量值分布較為分散，沒有呈現(xiàn)出明顯的分段特征，位點155附近的測量值分布非常分散。SOP-SeedBS算法的結(jié)果與Inspect的結(jié)果的相似度約為80%。

4 結(jié)論與展望

本文提出SOP-SeedBS方法結(jié)合了奇異值分解和投影的思想，通過分析高維稀疏變點數(shù)據(jù)的特點，對數(shù)據(jù)矩陣進行奇異值分解找到一個最優(yōu)的投影方向，將存在稀疏變點的高維數(shù)據(jù)沿著最優(yōu)的投影方向投影成一維數(shù)據(jù)，再利用SeedBS方法對一維數(shù)據(jù)進行多變點檢驗。通過仿真實驗，分別用變點檢驗頻數(shù)和平均Hausdorff距離評估了變點檢驗的準確度。本文提出的SOP-SeedBS方法在運算時間、算法準確度上都優(yōu)于Inspect方法，特別是在時間上具有較大的優(yōu)勢。最后將SOP-SeedBS方法應(yīng)用于ecp包中的膀胱腫瘤微陣列數(shù)據(jù)集中，檢驗結(jié)果與文獻[14]中給出的Inspect方法的檢驗結(jié)果有約80%的相似度。

本文是先將高維數(shù)據(jù)投影到一維，再進行變點檢驗。進一步還可以研究能否在保持數(shù)據(jù)維度不變的情況下，通過某種方式，例如引入一種針對高維數(shù)據(jù)的檢驗統(tǒng)計量，將SeedBS方法直接用于高維數(shù)據(jù)的變點檢驗問題中。