高斯過程回歸:高斯過程回歸（Gaussian Process Regre -百科知識中文網

歷史

與GPR有關的早期研究大致可分為3部分，在時間序列分析中，安德雷·柯爾莫哥洛夫（Андрей Николаевич Колмогоров）和羅伯特·維納（Robert Wiener）在二十世紀40年代提出了估計0均值平穩高斯過程信號的濾波技術。隨後出現的卡爾曼濾波（Kalman filter）提供了估計高斯隱變數（latent variable）的有效方法。在地統計學領域，1963年提出的克里金法（Kriging）首次實現了平穩隨機場的非參數回歸，其後在二十世紀90年代出現的貝葉斯克里金提供了與GPR接近或等價的高斯隨機場估計。在機器學習方面，包含無限節點的貝葉斯神經網路（Bayesian neural network）和各類核學習（kernel learning）方法例如核嶺回歸（kernel ridge regression）為GPR的核函式超參數求解帶來了啟發。

GPR作為機器學習一般方法的提出來自劍橋大學（University of Cambridge）學者Carl E. Rasmussen和愛丁堡大學（University of Edinburgh）學者Christopher K. I. Williams，二者在1996年發表的研究按函式空間觀點給出了GPR的求解系統並討論了GPR超參數的極大似然估計和蒙特卡羅方法求解。

理論

模型推導

這裡給出GPR在權重空間（weight-space）觀點和函式空間（function-space）觀點下的模型推導，前者由貝葉斯線性回歸（Bayesian linear regression）出發，通過特徵空間的映射得到GPR的預測形式，後者由高斯過程出發，以更簡潔的方式得到與前者等價的結果。

權重空間觀點

在權重空間（weight-space）觀點下，GPR可以由正態假設的貝葉斯線性回歸導出。具體地，給定相互獨立的N組學習樣本：，貝葉斯線性回歸是如下形式的多元線性回歸模型：

式中為權重係數，為殘差或噪聲。貝葉斯線性回歸假設殘差服從獨立同分布（independent and identically distributed, iid）的0均值常態分配：，由此可得貝葉斯線性回歸的似然：

在此基礎上，為模型權重賦予0均值正態先驗：。由貝葉斯定理（Bayes' theorem）可知，模型權重的後驗正比於似然和先驗的乘積：。由常態分配的共軛性（conjugacy）可知，對常態分配的似然和方差已知的常態分配先驗，其後驗也為常態分配，因此帶入似然和先驗的解析形式並按常態分配整理可得：

給定測試樣本，貝葉斯線性回歸通過對可通過邊緣化模型權重，即按其後驗積分得到測試結果的機率分布：

貝葉斯線性回歸是一個線性參數模型，為使其表示樣本間的非線性關係，可以使用給定的函式將映射至高維空間：

由於映射函式是固定的，即與模型權重無關，因此可直接帶入貝葉斯線性回歸的結果得到：

使用核方法（kernel method），即定義核函式可改寫上式：

上述改寫結果即是GPR對的均值和協方差進行預測的形式，式中表示由學習樣本得到的格拉姆矩陣（Gram matrix），即核矩陣，表示帶入測試樣本後得到的行向量和列向量：

函式空間觀點

對回歸模型，若函式的形式不是固定的，則其為潛函式（latent function）。潛函式的每個取值都是函式空間（function-space）的一個測度。GPR取該函式空間的先驗為高斯過程（Gaussian Process, GP），不失一般性，這裡表示為0均值高斯過程：

式中為學習樣本，其在高斯過程中的測度是高斯過程的有限維分布（finite-dimensional distribution），由定義可知，該有限維分布是聯合常態分配：

式中為核函式，0均值高斯過程由其核函式完全決定：

給定N組學習樣本，作為對算法的推導，假設回歸殘差服從iid常態分配：，則GPR在高斯過程先驗和常態分配似然下求解回歸模型的後驗：，並對測試樣本的測試結果進行估計（非正態似然的情形參見算法部分）。具體地，由回歸模型和高斯過程的定義，和的機率分布為：

因此二者的聯合機率分布是如下形式的聯合常態分配：

對上述聯合分布取的邊緣分布，由聯合常態分配的邊緣分布性質（marginalization property）可得：

上式是GPR的預測形式，也是函式空間後驗對測試樣本的有限維分布。式中協方差的展開使用了分塊矩陣求逆公式。注意到，上式與權重空間觀點下GPR的預測形式相同，即假設先驗和iid殘差為0均值常態分配的貝葉斯線性回歸在引入核函式後，等價於使用相同核函式和0均值高斯過程先驗的GPR 。

由上述預測形式的推導可知，GPR在假設0均值高斯過程先驗時，解得的後驗往往不是0均值的，因此0均值先驗具有泛用性。但作為說明，GPR也可使用均值不為0的高斯過程先驗，此時常見的做法是將潛函式表示為一組顯式基函式（explicit basis function）：。在給定基函式的常態分配先驗時，是具有如下形式的高斯過程：

參照先前的推導，此時GPR的預測形式為：

式中，具體步驟和討論可參見Rasmussen and Williams (2006), pp.27-30 。

核函式

GPR使用高斯過程作為先驗，即假設了學習樣本是高斯過程的採樣，因此其估計結果與核函式有密切聯繫。GPR中核函式的實際意義為協方差函式（covariance function），描述了學習樣本間的相關性，因此其不是通過核方法簡化計算的手段，而是模型假設的一部分。這裡對GPR常見的核函式進行簡單介紹，並引入核函式計算的加速方法。

核函式的選擇

核函式的選擇要求滿足Mercer定理（Mercer's theorem），即核函式在樣本空間內的任意格拉姆矩陣（核矩陣）為半正定矩陣（semi-positive definite）。

若GPR的先驗為平移不變（transformation invariant）的平穩高斯過程時，可用的核函式包括徑向基函式核（RBF kernel）、馬頓核（Matérn kernel）、指數函式核（exponential kernel）、二次有理函式核（rational quadratic kernel, RQ kernel）等，以RBF核為例，其解析形式如下：

式中為RBF核的超參數，表示其頻寬（width）和特徵長度尺度（characteristic length-scale），是表征各向異性的矩陣函式，式中的3種形式分別對應各向同性和兩類各向異性。常見地，對各向同性的情形，RBF核表示為：

GPR也可使用非平穩高斯過程，此時常見的核函式選擇為周期核（periodic kernel）與多項式函式核（polynominal kernel），二者分別賦予高斯過程周期性和旋轉不變性（rotation invariant）。

核矩陣求逆的計算簡化

由模型推導部分可知，GPR要求計算核矩陣的逆矩陣：，若使用Cholesky分解（Cholesky decomposition）求逆，則對N組學習樣本，其計算複雜度為。因此隨著學習樣本的增加，GPR的計算開銷會顯著增大。這裡簡單介紹核矩陣求逆的簡化方法以應對上述問題，一般地，這些方法除GPR外也適用於其它核學習（kernel learning）方法，例如核嶺回歸、支持向量機（Support Vector Machine, SVM）等。

1. 選取數據子集（subset of datapoints, SD）

由於GPR可以在少量學習樣本的情形下給出回歸問題的可靠估計，因此從學習樣本中合理選擇一個子集可以顯著降低核矩陣的大小，但不對GPR的表現造成顯著影響。選取SD的常見方法是通過尋找微分熵（differential entropy）最大的學習樣本作為支持向量定義SD，選取過程的計算複雜度為，在選取較小子集的情形下簡化了計算，使用SD最佳化的高斯過程建模被稱為稀疏高斯過程（sparse Gaussian process），或信息向量機（Informativevector machine, IVM）。

2. 低秩近似（low-rank approximation）

低秩近似將核矩陣近似為一系列低秩矩陣的乘積以簡化求逆計算，具體地，N×N的核矩陣可以近似表示為：，其中的大小分別為N×m、m×m、m×N，此時由Sherman-Morrison-Woodbury定理可得：

式中求逆部分的矩陣大小為m×m，因此計算複雜度減少至，在m取小值時簡化了計算。

確定低秩矩陣的常見方法包括Nyström法（Nyström method）、子集回歸法（subset ofregressors, SR）、映射過程近似（projected process approximation, PP）、BCM（Bayesian committee machine）等，這裡對前兩者進行介紹。Nyström法將核矩陣分解為4個分塊矩陣並在構建低秩近似後直接帶入GPR的預測形式： SR按與Nyström法相同的方式構建分塊矩陣，但其首先將GPR的預測形式改寫為核矩陣子集的回歸併帶入低秩矩陣，例如對GPR的均值部分有：。

算法

正態似然

GPR的求解也被稱為超參學習（hyper-parameter learning），是按貝葉斯方法通過學習樣本確定核函式中超參數的過程。由貝葉斯定理（Bayes' theorem）可知，GPR的超參數後驗有如下表示：

式中表示GPR的超參數，包括核函式的超參數和殘差的方差。為似然（liklihood），作為非參數模型，這裡的似然是對GPR的輸出邊緣化得到的邊緣似然（marginal liklihood）：由模型推導可知，在殘差服從iid常態分配時，GPR的似然服從常態分配：，此時其常見的求解方法是包含非線性最佳化的極大似然估計。

極大似然估計（Maximum Likelihood Estimation, MLE）

MLE通過令GPR的似然取極大值實現對超參數的單點估計。不失一般性，對0均值先驗的GPR，帶入聯合常態分配的解析形式並求其負自然對數可得：

式中第一項包含學習樣本，是數據擬合項，表示模型的經驗風險（empirical risk）；式中第二項僅與回歸模型有關，回歸模型的核矩陣越複雜其取值越高，表示模型的結構風險（structural risk），因此按統計學習理論，GPR是一個包含正則化的求解系統，但二者的比例，即殘差的方差由MLE給出。求解GPR似然極大值等價於求解負自然對數的極小值，因此這裡給出其對的偏導數：

常見的非線性最佳化算法均可以求解GPR的MLE問題，常見選擇為共軛梯度法（conjugate gradient method）和擬牛頓法（quasi-Newton method）。有研究使用隨機最佳化方法，例如隨機梯度下降（stochastic gradient descent）、遺傳算法（genetic algorithm）和粒子群算法（particle swarm optimization）求解GPR 。

GPR的對數似然不是凸函式，且其最佳化複雜度隨學習樣本的增加而增大，在大量學習樣本的情形下，可能會發現多個局部最優，且局部最優解的差異很大的情形，考慮解出的核函式超參數通常表示高斯過程的特徵長度尺度，多個局部最優意味著學習樣本可以按多個尺度進行回歸，其中小尺度的模型複雜度高但考慮更多局部特徵，大尺度的模型反之。

除MLE外，GPR也可通過極大後驗估計（Maximum A Posteriori estimation, MAP）求解。MAP是與MLE相近的估計方法，其最佳化目標是似然和先驗的乘積。因為MAP與MLE在正態後驗時得到的結果等價，在其他情形下的結果也相近，卻引入了先驗和額外的計算，因此在GPR研究中被較少提及。

極大偽似然估計（maximum pseudo-likelihood estimator, MPLE）

使用MPLE求解GPR的過程也被稱為交叉驗證（cross-validation），相比於MLE，其改變是對學習樣本使用留一法（Leave One Out, LOO）計算其偽似然（pseudo-likelihood）的負自然對數：

式中表示除以外的所有樣本，表示選取該矩陣的第i行第j列元素，由的形式可知其為留一數據由GPR估計的均值和方差，式中第2行為留一數據的似然，GPR的偽似然為所有留一數據似然的和。MPLE求解GPR在計算開銷上與MLE相當，因為兩者共同的計算效率瓶頸是核矩陣求逆步驟（參見理論部分）。在求解方法上，MPLE獨立第考察了每個學習樣本，因此其本質上是頻率估計（frequentist inference），有研究表明MPLE在GPR模型選擇不恰當（miss-specified）時得到的結果更好，MLE在模型選擇恰當（well-specified）時得到的結果更好。

編程實現

這裡提供一個使用GPy封裝的GPR模型的例子：

非正態似然

在似然不服從常態分配，例如處理異方差噪聲數據時，GPR對測試數據的後驗沒有解析形式，此時可以使用的解析近似（analytical approximation）方法，例如拉普拉斯近似（Laplace Approximation, LA）和期望傳播（Expectation Propagation, EP）將非正態後驗近似表示為正態後驗。LA和EP的推導可參見Rasmussen and Williams (2006), pp.41-60 。

數值方法可用在非正態似然，或核函式超參數過多，不利於MLE問題最佳化的情形下求解GPR。選擇包括蒙特卡羅方法（Monte Carlo）和變分貝葉斯估計（Variational Bayesian Inference, VBI）。VBI使用Jensen不等式得到GPR對數似然的凸函式下界作為代理損失並使用梯度算法求解，其優勢是計算複雜度接近，顯著簡化了計算。蒙特卡羅方法被認為是在GPR數值解法中擁有高準確度和泛用性的方法，其原理是使用大量計算開銷疊代生成隨機數逼近GPR後驗。蒙特卡羅求解GPR的常見方法為混合蒙特卡羅（Hybrid Monte Carlo, HMC），其步驟接近於Metropolis-Hastings算法，但通過“動量”項限制了隨機遊走的取值。HMC求解GPR的例子可參見Rasmussen (1996) 。

擴展算法

這裡給出與高斯過程回歸有關的擴展算法，但更一般地，這些方法是高斯過程模型的擴展，除回歸問題外也可被套用於其它機器學習主題。

封裝高斯過程（Warped Gaussian Process, WGP）

WGP通過對GPR進行非線性變換將其拓展至非高斯過程的情形，具體地，WGP首先在特徵空間選擇一個單調封裝函式（monotonic warping function），將學習樣本變換至封裝函式指定的潛空間後再求解包含封裝函式超參數的極大似然。WGP的中封裝函式的選擇通常為雙曲正切函式（hyperbolic tangent）線性組合，這裡給出WGP對應於GPR的模型和算法：

式中為預先給定的WGP封裝次數，在求解超參數後，類比GPR可得到WGP的預測形式。除上述形式的封裝函式外，使用更複雜的封裝函式，例如可進一步提升WPG的靈活性。在一些測試數據中，WGP的表現優於GPR 。

半參數高斯過程（Semi-parametric Gaussian Processes, SGP）

SGP是在回歸問題中將高斯過程模型與參數模型線性組合以獲得兩者優點——GPR的靈活性和參數模型在高維問題中利用數據的有效性——的算法。具體地，SGPR構建了如下求解系統：

式中為參數回歸部分，常見形式為線性回歸，第二項為0均值高斯過程先驗的GPR，最後一項為殘差。SGPR的學習分為兩步，先對學習樣本使用參數模型進行回歸，隨後使用GPR擬合參數模型的殘差。SGP被使用的常見情形是在已知學習數據和學習目標的確定性關係時考慮噪聲帶來的隨機效應，例如自動控制系統的預測問題。

深度高斯過程（deep Gaussian Process, DGP）

由性質可知，高斯過程可以視為一個擁有單隱含層和無限個隱含層節點的多層感知器（Multi-Layer Perceptron, MLP），DGP是MLP由單隱含層推廣至多隱含層時得到的高斯過程，即擁有多個隱含層和無限寬度的MLP。按定義，DGP是一組從高斯過程先驗中選擇的函式的分布：

式中表示第l層的輸出。DGP在結構上是由非參數的高斯過程基函式表示隱含層的神經網路，並通過樣本求解核函式超參數。DGP的另一種等價觀點/結構是控制MLP中的多層無限節點保持不變並在每個隱含層後加入一組加權求和，由樣本求解求和部分，即的權重。

可加高斯過程（Additive Gaussian Process, AGP）

GPR的常見核函式，例如RBF核、馬頓核等，被認為在高維特徵空間的學習問題中泛化能力較差，而AGP即是為解決GPR高維學習問題而提出的擴展方法。AGP利用核函式性質，在構建高斯過程先驗時對核函式進行了組合和結構最佳化。介紹性地，AGP在不同維度將不同結構的核函式相加作為先驗，並通過貝葉斯階層核學習（Bayesian Hierarchical Kernel Learning, HKL）更新超參數以去除對學習樣本解釋不利的結構。以RBF核為例，其1階可加高斯過程先驗可表示為：

式中表示數據維度。作為核函式結構的最佳化方法，AGP的特點在於其構建的核函式具有可解析的參數化形式，容易對模型進行解釋。在回歸問題中，AGP被認為具有與其它現代機器學習算法相當的預測能力。

有關概念

高斯過程分類（Gaussian Process Classification, GPC）

GPC是使用高斯過程作為先驗的分類器，在二元分類中，GPC在對潛函式進行估計後將其作為Sigmoid函式的輸入得到分類機率；在多元分類中，則將Sigmoid函式替換為歸一化指數函式。以二元分類為例，因為標籤數據是取值的二項分布隨機變數，所以對應的GPC似然是潛函式對學習樣本的因子乘積：

GPC通過最大化似然求解超參數並得到後驗，帶入Sigmoid函式的表達式可知，上式不是常態分配，因此不同於GPR，GPC的後驗沒有解析形式，要求使用非正態似然的求解方法，例如解析近似和蒙特卡羅方法。

克里金法（Kriging）

克里金法是與GPR相近的非參數模型，常見於隨機場的插值問題。若協方差函式的形式等價，簡單克里金（simple Kriging）和普通克里金（ordinary Kriging）的輸出與GPR在正態似然下輸出的數學期望相同，若克里金法使用高斯隨機場假設，則給出的置信區間也與GPR相同。克里金法與GPR的不同點在於，前者假設隨機場為固有平穩過程（intrinsically stationary process）並給出其對測試樣本的最優無偏估計（Best Linear Unbiased Prediction, BLUP）；後者假設隨機場為高斯過程並給出其對測試樣本的完整後驗。此外，考慮在地統計學中套用的常見情形，克里金法通常不加入噪聲，其估計結果在學習樣本處與學習目標完全相同。

有一些克里金法的擴展方法與GPR等價，例如Handcock and Stein (1993) 和Handcock and Wallis (1994) 提出的包含貝葉斯推斷的克里金法使用了各向同性的馬頓核高斯過程先驗，並在正態似然假設下按MLE求解超參數。

性質

在使用特定核函式時，GPR是一個通用近似（universal approximator），具體地，若實數域內有緊緻空間和賦范向量空間（normed vector space），則可以證明，對任意函式，等價核（Equivalent Kernel, EK）構成的再生核希爾伯特空間（Reproducing Kernel Hilbert Space, RKHS）內總是存在能夠按任意精度逼近原函式，即RKHS內總是存線上性組合：。等價核的例子包括RBF核、二次有理函式核、指數函式核等。

作為非參數高斯過程模型的性質：GPR的複雜程度取決於學習樣本，因此天然避免了過擬合（overfitting）問題。當測試樣本與學習樣本在特徵空間的距離趨於無窮時，GPR的估計結果趨於其高斯過程先驗，例如對0均值的各向同性高斯過程先驗有：。由於高斯過程可以通過協方差函式模擬隨機變數間的相關性，因此GPR可套用於平移不變、旋轉不變或周期性數據的回歸問題，但同時，由於核函式/非參數方法的局限，GPR在高維數據的學習中表現不佳，該現象有時被稱為“維數詛咒（curse of dimensionality）” 。此外由理論部分可知，GPR的計算複雜度較高。一些研究給出了計算量更低的改進算法並取得了效果，但總體而言，GPR是一個為小樣本設計的機器學習算法。

作為貝葉斯方法的性質：GPR是一個包含全貝葉斯特性（full Bayesian）的回歸模型，可以給出測試數據的完整後驗。此外，在似然為常態分配時，GPR給出的後驗是具有閉合形式的聯合常態分配，即GPR具有可解析性，該性質在非參數模型中並不多見。

與核學習方法的比較：GPR與適用於回歸問題的核學習（kernel learning）方法，例如核嶺回歸、支持向量回歸（support vector regression）等都使用了核函式與核方法，且可以通過相同的方式簡化核矩陣求逆計算，但對後者，其核函式是輸入空間向高維特徵空間的映射，而GPR中的核函式是高斯過程協方差函式的模型，表示隨機變數間的相關性。上述不同在求解中的體現是，GPR的“學習”是通過數據確定核函式超參數的過程，而核方法的“學習”，是在給定核函式超參數後求解模型權重或尋找支持向量的過程。

與人工神經網路（Artificial Neural Network, ANN）的關係：ANN是參數模型，在理論和結構上與GPR有較大差異，但二者在套用中有重疊，例子包括時間序列分析和計算機視覺問題。若學習樣本充足，ANN由於使用隨機梯度下降方法進行學習，在求解效率上優勢明顯；但在要求給出機率輸出時，GPR是更合適的方法。ANN中的多層感知器（Multi-Layer Perceptron, MLP）與高斯過程存在聯繫，高斯過程可以由擁有單隱含層和無限個隱含層節點的MLP導出。具體地，對如下預測形式的MLP：若其權重為iid隨機變數且均值為0，方差為有限值，則由中心極限定理（central limit theorem）可推出，當隱含層節點數n趨於無窮時，模型輸出的聯合分布趨於0均值的聯合常態分配：該結論不要求iid權重服從常態分配，但當權重服從iid常態分配時，有限個隱含層節點的MLP也可表示為高斯過程。另一方面，MLP也可以由高斯過程導出，由Mercer定理可知，核函式可以表示為映射函式的內積：，因此可視為以為隱含層特徵（輸出）的MLP 。

套用

GPR在各類低維的回歸問題中均可套用，尤其是時間序列數據的預測，例子包括太陽輻射的有關變數、風速、土壤溫度、冒納羅亞觀測站（Mauna Loa Observatory，MLO）的全球對流層平均二氧化碳濃度等。在圖像處理（image processing）方面，GPR被用於圖像去噪和生成超解析度圖像。在自動控制方面，GPR被用於機械臂（robotic arm）數據的實時學習，也有研究開發了GPR的機器人學習（robot learning）系統。

包含GPR的編程模組

基於Python開發的機器學習模組scikit-learn提供封裝的GPR工具GaussianProcessRegressor，其求解方法為BFGS算法（Broyden-Fletcher-Goldfarb-Shanno）的MLE，擁有自定義的求解器函式接口。在核函式方面，scikit-learn包含RBF核、馬頓核、二次有理核、指數-周期核、多項式核以及自定義核函式的API工具。

由謝菲爾德機器學習團隊（Sheffield machine learning group）開發的Python高斯過程建模工具GPy提供了GPR的完整求解系統和核函式，包括異方差噪聲GPR、隱變數模型，解析近似、VBI和蒙特卡羅求解、稀疏高斯過程和低秩近似。此外有基於GPy的模型超參數最佳化工具GPyOpt可用。

pyGPs是基於Python開發的面向對象高斯過程建模套用，包含可視化學習界面，在功能上包括了常見核函式、稀疏高斯過程和解析近似。

高斯過程回歸

歷史

理論

模型推導

核函式

算法

正態似然

非正態似然

擴展算法

有關概念

性質

套用

相關詞條

高斯過程

高斯數學

克里金法

平成系奧特曼劇場版

隨機信號分析基礎（第4版）

閃電俠[2014年美國CW新超級英雄系列電視劇]

機率論與數理統計

卡爾曼濾波

通信信號分析與處理

熱門詞條