基于SOP-SeedBS的高維數據稀疏變點檢驗

銀炳皓，施三支，齊德全

（長春理工大學 數學與統計學院，長春 130022）

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

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

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

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

1.1 稀疏變點模型

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

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

假設存在υ個變點，變點檢驗的問題能被表示為經典的假設檢驗問題

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

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

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

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

1.2 稀疏最優投影（SOP）

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

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

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

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

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

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

則有：

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

懲罰函數可以將有約束的目標函數轉化為無約束的目標函數，即當得到的M不滿足約束條件時，估計的?可能會遠離最優解。通過添加一個懲罰項 -λ‖M‖1，優化問題變為以下形式：

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

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

將高維數據X沿著最優方向投影，得到一維數據，記為Y。

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

1.3 種子二元分割（SeedBS）

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

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

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

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

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

用于種子二元分割的檢驗統計量是定義一段內分割點在s的CUSUM統計量：

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

其中：

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

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

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

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

1.4 SOP-SeedBS方法

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

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

（1）對數據矩陣X進行分解并推導出最優投影方向。

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

（4）根據衰減參數a構造包含m個種子區間的區間集。

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

2 數值模擬

構造由模型1生成的具有均值變點的高維數據時，設置參數σ2=1，數據量n分別取200，500，1 000，2 000，維數p分別為200，500，1 200，將高維數據噪聲表示為ε，分別取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控制著發生變化維度數量的稀疏程度，取值分別為10，40。根據文獻[10]的建議，。根據文獻[12]的建議，設置衰減參數a= 2。另外，下面三種情形分別考慮了均值結構中發生變化的維數坐標(p)上的重疊問題，高維數據中均值結構發生變化的三種不同重疊情形說明如下：

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

情形3（S3）：在完全重疊的情況下，則對于相應的變點，變化的維數發生在相同的k個坐標上。

本文先給出了SOP-SeedBS方法在情形3下，參數設置為σ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}時，不同噪聲下的變點估計結果可視化，用ε表示高維數據的噪聲，不同ε取值對于投影后一維數據的變點檢驗結果如圖2所示，圖2的結果顯示了當ε的取值越小時，噪聲越小，分段特征越明顯。另外，當均值跳躍越大，分段特征也越明顯。后續的仿真實驗中，噪聲ε皆設置為1。

圖2 不同 ε下的投影后一維數據的變點檢驗結果

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

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

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

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

圖3給出了SOP-SeedBS方法與Inspect方法在其他參數不變的情況下，不同的n和p設置下算法平均運算時間的比較，可以看出SOP-SeedBS方法在算法運算速度上具有較大的優勢。圖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中的結果顯示，當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%，由此能體現SOP-SeedBS方法在運算時間上優于Inspect方法。

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

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

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

3 實證分析

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

圖4 前10個個體在前200個位點上的ACGH數據散點圖

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

4 結論與展望

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

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