GPR モデルのブロック座標降下近似

観測値の数が多い場合、新しいデータによる予測とパラメーター推定に厳密法を使用すると、計算に時間がかかる可能性があります (厳密 GPR 法を参照してください)。予測の場合にこの問題への対処に役立つ近似法の 1 つは、ブロック座標降下 (BCD) 法です。BCD 法を使用して予測を行うには、名前と値の引数 PredictMethod="bcd" で予測法を指定して fitrgp を呼び出してから、関数 predict を使用します。

BCD 法の考え方では、次を

$$\widehat{\alpha }={\left(K(X,X)+{\sigma }^2{I}_N\right)}^{-1}\left(y-H\beta \right)$$

厳密法とは異なる方法で計算します。BCD では、次の最適化問題を解くことにより $$\widehat{\alpha }$$ を推定します。

$$\widehat{\alpha }={\text{arg min}}_{\alpha }f(\alpha )$$

ここで

$$f(\alpha )=\frac{1}{2}{\alpha }^T\left[K(X,X)+{\sigma }^2{I}_N\right]\alpha -{\alpha }^T\left(y-H\beta \right).$$

fitrgp は、ブロック座標降下 (BCD) を使用してこの最適化問題を解きます。サポート ベクター マシン (SVM) に精通していると、SVM の当てはめに使用する逐次最小最適化 (SMO) と反復単一データ アルゴリズム (ISDA) が BCD の特殊なケースであることがわかります。fitrgp は、ブロック選択およびブロック更新という 2 つのステップで各 BCD を更新します。ブロック選択フェーズでは、fitrgp はランダムな方法とグリーディな方法を組み合わせて、次を最適化するための一連の係数 $$\alpha$$ を選択します。ブロック更新フェーズでは、他の係数 $$\alpha$$ を前の値に固定したまま、選択した係数 $$\alpha$$ が最適化されます。この 2 つのステップは、収束するまで繰り返されます。BCD では、$$n*n$$ の行列 $$K\left(X,X\right)$$ をメモリに格納することはなく、必要に応じて行列 $$K\left(X,X\right)$$ の一部を計算するだけです。この結果、BCD は PredictMethod="exact" よりメモリ効率が高くなります。

実際の経験によると、$$\alpha$$ を計算するための最適化問題を非常に高い精度で解く必要はありません。たとえば、$$\left|\right|\nabla f(\alpha )\left|\right|$$ が $$\eta \left|\right|\nabla f({\alpha }_0)\left|\right|$$ より小さくなったときに反復を停止すれば十分です。ここで、$${\alpha }_0$$ は $$\alpha$$ の初期値、$$\eta$$ は小さい定数 ($$\eta ={10}^{-3}$$ など) です。BCD では $$\alpha$$ を異なる方法で計算することを除き、PredictMethod="exact" と PredictMethod="bcd" の結果は同じです。したがって、PredictMethod="bcd" は高いメモリ効率で PredictMethod="exact" を実行するための方法であると考えることができます。n が大きい場合、BCD はしばしば PredictMethod="exact" より高速です。ただし、PredictMethod="bcd" オプションを使用して予測区間を計算することはできません。これは、予測区間を計算するには、すべてのテスト点について $$\alpha$$ を計算するために解くものと同じような問題を解く必要があり、計算に時間がかかるためです。

BCD 近似の使用による GPR モデルの当てはめ

この例では、ブロック座標降下 (BCD) 近似を使用して、多数の観測値があるデータにガウス過程回帰 (GPR) モデルを当てはめる方法を示します。

小さい n - PredictMethod "exact" と "bcd" は同じ結果を生成

小さい標本データ セットを生成します。

rng(0,"twister")       
n = 1000;    
X = linspace(0,1,n)';
X = [X,X.^2];    
y = 1 + X*[1;2] + sin(20*X*[1;-2])./(X(:,1)+1) + 0.2*randn(n,1);

FitMethod="exact" と PredictMethod="exact" を使用して GPR モデルを作成します。

gpr = fitrgp(X,y,KernelFunction="squaredexponential", ...
    FitMethod="exact",PredictMethod="exact");    

FitMethod="exact" と PredictMethod="bcd" を使用して別の GPR モデルを作成します。

gprbcd = fitrgp(X,y,KernelFunction="squaredexponential", ...
    FitMethod="exact",PredictMethod="bcd",BlockSize=200);

PredictMethod="exact" と PredictMethod="bcd" は等価であるはずです。これらは、同じ $$\underset{}{\overset{ˆ}{\alpha }}$$ を異なる方法で計算しただけです。fitrgp を呼び出すときに名前と値の引数 Verbose=1 を使用すると、各反復を表示することもできます。

2 つの方法を使用して当てはめた値を計算し、比較します。

ypred = resubPredict(gpr);
ypredbcd = resubPredict(gprbcd);
max(abs(ypred-ypredbcd))

ans = 5.8853e-04

当てはめた値は互いに近くなっています。

次に、PredictMethod="exact" および PredictMethod="bcd" の当てはめた値とともにデータをプロットします。

figure
plot(y,".")
hold on
plot(ypred,"-",LineWidth=2)
plot(ypredbcd,"--",LineWidth=2)
hold off
legend("Data","PredictMethod Exact","PredictMethod BCD", ...
    Location="best")
xlabel("Observation index")
ylabel("Response value")

PredictMethod="exact" と PredictMethod="bcd" はほぼ同じ近似を生成することがわかります。

大きい n - FitMethod="sd" と PredictMethod="bcd" を使用

前のものと似ている大きいデータ セットを生成します。

rng(0,"twister")        
n = 50000;    
X = linspace(0,1,n)';
X = [X,X.^2];    
y = 1 + X*[1;2] + sin(20*X*[1;-2])./(X(:,1)+1) + 0.2*randn(n,1);

n = 50000 の場合、行列 K(X, X) は 50000 行 50000 列になります。K(X, X) をメモリに格納するには、約 20 GB の RAM が必要です。このデータ セットに GPR モデルを当てはめるため、m = 2000 個の点によるランダムなサブセットを選択して FitMethod="sd" を使用します。fitrgp の呼び出しにおける名前と値の引数 ActiveSetSize では、アクティブ セットのサイズ "m" を指定します。$$\alpha$$ を計算するため、BlockSize を 5000 にして PredictMethod="bcd" を使用します。fitrgp の名前と値の引数 BlockSize では、各反復で最適化するベクトル $$\alpha$$ の要素数を指定します。5000 という BlockSize は、使用しているコンピューターで 5000 行 5000 列の行列をメモリに格納できると仮定しています。

gprbcd = fitrgp(X,y,KernelFunction="squaredexponential", ...
    FitMethod="sd",ActiveSetSize=2000,PredictMethod="bcd",BlockSize=5000);     

データと当てはめた値をプロットします。

figure
plot(y,".")
hold on    
plot(resubPredict(gprbcd),"-",LineWidth=2)
hold off
legend("Data","PredictMethod BCD",Location="best")
xlabel("Observation index")
ylabel("Response value")

プロットは、"n" が小さい場合のものと同じようになります。

参照

[1] Grippo, L. and M. Sciandrone. "On the convergence of the block nonlinear Gauss-Seidel method under convex constraints." Operations Research Letters. Vol. 26, pp. 127–136, 2000.

[2] Bo, L. and C. Sminchisescu. "Greed block coordinate descent for large scale Gaussian process regression." In Proceedings of the Twenty Fourth Conference on Uncertainty in Artificial Intelligence (UAI2008): http://arxiv.org/abs/1206.3238, 2012.

参考

fitrgp | predict

関連するトピック

• ガウス過程回帰モデル
• 厳密 GPR 法
• GPR モデルのデータ サブセット近似
• GPR モデルの回帰変数サブセット近似
• GPR モデルの完全独立条件近似