\(Ax = b\) denklem sistemini çözmek için
\[ x = (A^TA)^{-1}A^Tb \qquad (1) \]
formülü en az kareler çözümü olarak biliniyor, bkz [2]. Bu çözüm iyi işler, fakat bazı durumlarda negatif bir tarafı var, çözüm toptan (batch) olarak yapılıyor. \(A\) içinde 100 tane satır olabilir, ona göre bir çözüm bulunur, ardından 1 tane ek veri satırı gelirse olsa 101 tane satır için tüm işlemlerin tekrar baştan yapılması gerekir. Acaba sadece o yeni verilen satır için önceki \(x\) tahminini bir şekilde güncellemek mümkün mü?
Özyineli en az kareler ile bunu yapabiliriz. Diyelim ki
\[ c_1 t + c_2 = b \]
lineer sistemini çözmek istiyoruz, yani bu bir çizgi uydurma (line fitting) olacak, kesi \(c_2\), eğim \(c_1\). Notasyon altta, altsimge \(k\) kaç tane veri satırı olduğunu gösterecek,
\[ A_kx_k \approx b_k, \quad A_k = \left[\begin{array}{cc} t_1 & 1 \\ t_2 & 1 \\ \vdots & \vdots \\ t_k & 1 \end{array}\right], \quad x_k = \left[\begin{array}{r} c_{1,k} \\ c_{2,k} \end{array}\right], \quad b_k = \left[\begin{array}{r} B_1 \\ B_2 \\ \vdots \\ B_k \end{array}\right] \]
Eğer tek istediğimiz tek boyutlu bir zaman serisi için çizgi uydurma yapmak ise \(t_1,..,t_k\) 1 ve \(k\) arası tam sayılar olurdu, bu durumda \(A_k\) iyice basitleşir. Devam edelim, eğer (1)'i üstteki format için uyarlarsak,
\[ x_k = (A_k^TA_k)^{-1}A_k^T b_k \qquad (5) \]
Yani elde \(k\) tane veri var, üstteki formülü uyguladık ve bir çözüm bulduk. Şimdi diyelim ki yeni ölçümler \((t_{k+1}, B_{k+1})\) aldık, ve
\[ x_{k+1} = (A_{k+1}^TA_{k+1})^{-1}A_{k+1}^T b_{k+1} \qquad (3) \]
hesabını yapmamız lazım. Ek notasyon;
\[ A_{k+1} = \left[\begin{array}{c} A_k \\ a_{k+1}^T \end{array}\right], \quad a_{k+1}^T = \left[\begin{array}{c} t_{k+1} \\ 1 \end{array}\right], \quad b_{k+1} = \left[\begin{array}{c} b_k \\ B_{k+1} \end{array}\right], \quad P_k = (A_k^TA_k)^{-1} \qquad (4) \]
Matris tersi \(P_k\)'nin yeni veri noktası gelince nasıl güncellendiğini görelim,
\[ P_{k+1} = (A_{k+1}^TA_{k+1})^{-1} = \bigg[ \left[\begin{array}{cc}A_k & a_{k+1} \end{array}\right] \left[\begin{array}{c}A_k \\ a_{k+1}^T \end{array}\right] \bigg]^{-1} \]
Eşitliğin sağındaki matris çarpımını yaparsak, ve \(P_k\)'yi yerine koyarsak,
\[ = [ A_k^TA_k + a_{k+1}a_{k+1}^T ]^{-1} = [ P_k + a_{k+1}a_{k+1}^T ]^{-1} \qquad (2) \]
Üstte yine sağdaki formül \((A+BCD)^{-1}\) formunda bir ters alma işlemi gibi gözüküyor; Matris Tersi Yardımcı Teorisi (Matrix Inversion Lemma) diyor ki [1, sf. 469], herhangi bir \(A,B,C,D\) için,
\[ [A + BCD]^{-1} = A^{-1} - A^{-1}B[C^{-1} + DA^{-1} B]^{-1} DA^{-1} \]
(2)'deki ifadenin üstteki forma göre paylaştırmasını şöyle yapalım, \(A = P_k\), \(B = a_{k+1}\), \(C=I\), \(D=a_{k+1}^T\). Buna göre (2) üstteki açılım üzerinden ve paylaştırılan sembollere göre şu hale gelir,
\[ P_{k+1} = P_k - P_k a_{k+1}(I + a_{k+1}^T P_k a_{k+1})^{-1} a_{k+1}^TP_k \]
Parantez içindeki büyük çarpım bir tek sayı olduğu için \(I\) değeri 1 yapılabilir,
\[ P_{k+1} = P_k - P_k a_{k+1}(1 + a_{k+1}^T P_k a_{k+1})^{-1} a_{k+1}^TP_k \qquad (6) \]
Bu durumda tersi alınan parantez içindeki tüm ifade de tek sayı demektir, ve bu tek sayının tersini almak çok basittir (\(x\) için \(1/x\)).
Nihai güncelleme formülü için devam edelim; (3) formülüne (4)'teki eşitlikleri koyalım,
\[ x_{t+1} = P_{k+1} \left[\begin{array}{cc} A_k^T & a_{k+1} \end{array}\right] \left[\begin{array}{c} b_k \\ B_{k+1} \end{array}\right] \]
\[ = P_{k+1} [A_k^Tb_k + a_{k+1}B_{k+1} ] \]
\[ (A_k^TA_k) x_k = A_k^T b_k \]
Bu sonucu iki üstteki formüle sokarsak,
\[ = P_{k+1} [A_k^TA_kx_k + a_{k+1}B_{k+1} ] \]
(4)'teki formlar üzerinden
\[ A_{k+1}^TA_{k+1} = A_k^TA_k + a_{k+1}a_{k+1}^T \]
diyebileceğimizi görmüştük, o zaman
\[ A_{k+1}^TA_{k+1}x_k = (A_k^TA_k + a_{k+1}a_{k+1}^T)x_k \]
Üç üstteki formülde yerine koyalım,
\[ = P_{k+1} [(A_k^TA_k + a_{k+1}a_{k+1}^T)x_k + a_{k+1}B_{k+1} ] \]
\[ = P_{k+1} [P_{k+1}^{-1}x_k + a_{k+1}a_{k+1}^Tx_k + a_{k+1}B_{k+1} ] \]
\[ x_{k+1} = x_k + P_{k+1}a_{k+1}a_{k+1}^Tx_k + P_{k+1}a_{k+1}B_{k+1} \]
\[ x_{k+1} = x_k + P_{k+1}a_{k+1}(a_{k+1}^Tx_k + B_{k+1}) \qquad (7) \]
Şimdi \(P_{k+1}\)'yi özyineli olarak temsil etmek şunları yapalım. \(K_{k+1} = P_{k+1}a_{k+1}\) sistemin kazanç matrisi (gain matrix) olsun, ve (6)'daki \(P_{k+1}\) eşitliği kullanarak formülü genişletelim,
\[ K_{k+1} = P_{k+1}a_{k+1} = [ P_k - P_k a_{k+1} [ 1 + a_{k+1}^T P_k a_{k+1} ]^{-1} a_{k+1}^TP_k ] a_{k+1} \]
\[ = P_ka_{k+1} - P_k a_{k+1}[a_{k+1}^T P_k a_{k+1} + 1]^{-1} a_{k+1}^TP_ka_{k+1} \]
\[ = P_ka_{k+1} \big[ I - [ a_{k+1}^T P_k a_{k+1} + 1 ]^{-1} a_{k+1}^TP_ka_{k+1} \big] \]
Eğer bu formülü aynı anda hem \((a_{k+1}^TP_ka_{k+1})\) hem de \((a_{k+1}^TP_ka_{k+1})^{-1}\) ile çarparsak (hiçbir etkisi olmayan bir işlem, birbirini iptal ediyor çünkü) bazı temizleme işlemlerini yapmak mümkün olur,
\[ = P_ka_{k+1} \big[ (a_{k+1}^T P_k a_{k+1} + 1) - a_{k+1}^TP_ka_{k+1} \big] (a_{k+1}^T P_k a_{k+1} + 1)^{-1} \]
Büyük parantez içinde sadece +1 sağ kalır, geri kalanlar iptal olur,
\[ K_{k+1} = P_ka_{k+1} (a_{k+1}^T P_k a_{k+1} + 1)^{-1} \]
Bu formülü (7) içine geri \(K_{k+1}\) olarak koyarsak,
\[ x_{k+1} = x_k + K_{k+1}(a_{k+1}^Tx_k + B_{k+1}) \]
Aynı şekilde (6) içine koyarsak,
\[ P_{k+1} = P_k - \underbrace{P_k a_{k+1}(1 + a_{k+1}^T P_k a_{k+1})^{-1}}_{K_{k+1}} a_{k+1}^TP_k \]
\[ P_{k+1} = P_k - K_{k+1}a_{k+1}^TP_k \]
Böylece \(K_{k+1},P_{k+1},x_{k+1}\) özyineli güncelleme formüllerini elde etmiş oluyoruz.
Kodlar
Güncelleme kodları alttadır,
import numpy as np
def rlse_online(aT_k1,b_k1,x,P):
K = np.dot(P,aT_k1.T)/(np.dot(np.dot(aT_k1,P),aT_k1.T)+1)
x = x +K*(b_k1-np.dot(aT_k1,x))
P = P-np.dot(K,np.dot(aT_k1,P))
return x,K,PÖrnek olarak alttaki veriyi kullanalım.
import numpy.linalg as lin
b = np.array([[3.0,4.0,6.0,3.0,8.0,7.0,5.0]]).T
A= np.ones((len(b),2)); A[:,1] = range(len(b))Özyineli olarak problemi çözelim; her veri noktasını teker teker güncelleme rutinine geçelim.
import rls
n = 2
P = np.eye(n,n)*100.
x = np.zeros((n,1))
for k in range(len(b)):
x,K,P = rls.rlse_online(np.array([[k,1]]),b[k,:],x,P)
print x[[ 0.5037057 ]
[ 3.62655923]]Üstteki sonuç bulundu. Şimdi aynı verileri en az kareler ile toptan şekilde çözelim,
import statsmodels.api as sm
y = b; x = A
f = sm.OLS(y,x).fit()
print f.params[ 3.64285714 0.5 ]Önce Toptan, Sonra Özyineli
Eğer verinin bir kısmı için toptan başlayıp sonra özyineli gitmek istersek ne yaparız? O zaman elde bir \((A_k^TA_k)^{-1}\), yani \(P_{k}\) olurdu, toptan şekilde hesaplanmış olacaktı, ve bu değerin sonraki hali için güncelleme formülünü biliyoruz, böyle devam ederdik. Tabii bu durumda \((A_k^TA_k)^{-1}\)'yi toptan hızlı hesaplamak için bir teknikten bahsetmek lazım, en az kareler rutinleri genelde bu değeri geri döndürmezler, {}'dan hatırlarsak bu hesabı direk yapmak oldukça pahalı, o yüzden QR bazlı bir yaklaşım lazım (aynen \(x\)'in kendisinin QR bazlı hesaplandığı gibi). Her \(A_k\) matrisinin bir \(A_k = QR\) açılımı olacağından hareketle,
\[ A_k^TA_k = (QR)^TQR = R^TQ^TQR = R^TR \]
O zaman
\[ (A_k^TA_k)^{-1} = (R^TR)^{-1} = R^{-1}R^{-T} \]
Şimdi verinin en son satırı hariç ilk kısmı üzerinde bu değeri hesaplayalım,
A_k = A[:-1,:]
b_k = b[:-1,:]
print A.shape, A_k.shape
q,r = lin.qr(A_k)
Pk_r = np.dot(lin.inv(r), lin.inv(r.T))
print Pk_r
Pk = lin.inv(np.dot(A_k.T,A_k))
print Pk(7, 2) (6, 2)
[[ 0.52380952 -0.14285714]
[-0.14285714 0.05714286]]
[[ 0.52380952 -0.14285714]
[-0.14285714 0.05714286]]Direk usül ve QR bazlı ters işleminin aynı sonuçlara erişildiğini görüyoruz. Toptan \(x_k\)
x_batch = np.dot(np.dot(lin.inv(r), q.T), b_k)
print x_batch.T[0][ 3.0952381 0.82857143]Şimdi yeni veri noktası ile güncelleyelim,
A_new = A[-1,:]
b_new = b[-1,:]
x_new,K_new,P_new = rls.rlse_online(A_new,b_new,x_batch.T[0],Pk_r)
print x_new[ 3.64285714 0.5 ]Aynı sonuca eriştik.
Kaynaklar
[1] Yang, Applied Numerical Methods using Matlab
[2] Bayramlı, Lineer Cebir, Ders 16