Fast Domain Generalization with Kernel Methods – Part 1


Sekarang saya sedang melakukan riset tentang bagaimana mempercepat algoritma machine learning dalam konteks domain adaptation/generalization (bisa lihat di bab 2 pada tulisan saya sebelumnya untuk memahami istilah tersebut). Basis dari ide project kali ini berdasarkan paper-paper berikut:

  1. Domain Generalization via Invariant Feature Representation, ICML 2013
  2. Using the Nystrom Method to Speed Up Machines, NIPS 2000
  3. Correlated random features for fast semi-supervised learning, NIPS 2013
  4. A Simple Algorithm for Semi-supervised Learning with Improved Generalization Error Bound, ICML 2012
  5. Random Features for Large-Scale Kernel Machines, NIPS 2007

Rencananya adalah bagaimana mempercepat metode [1] dengan salah satu atau kombinasi dari metode [2] – [5].

Setelah saya eksplorasi, secara umum inti dari metode-metode dibahas pada [2]-[5] adalah mengaproksimasi matriks kernel dengan matriks yang berukuran lebih kecil. Oleh karena itu, penting untuk memahami apa itu kernel berikut properti-propertinya, yang merupakan komponen utama dari berbagai algoritma machine learning arguably paling efektif saat ini, seperti Support Vector Machines, Kernel PCA, Kernel Ridge Regression, etc). Kumpulan dari algoritma-algoritma berbasis kernel belakangan ini populer dengan istilah kernel machines.

Walaupun metode kernel ini bukan hal yang benar-benar baru bagi saya, dalam 1 minggu pertama sebagian besar waktu saya habiskan untuk mengeksplorasi kernel agar benar-benar memahami sifat-sifatnya lebih lanjut.

Di tulisan ini saya ingin membahas definisi beserta beberapa contoh fungsi kernel, dan aplikasinya pada algoritma ridge regression.

1. Kernels & Feature map

Yang dimaksud dengan kernel dalam konteks ini adalah sejenis fungsi kesamaan (similarity measure) antara 2 buah objek (bukan kernel di sistem operasi ataupun di konteks lainnya) yang dapat didefinisikan sebagai berikut

\displaystyle k: \mathcal{X} \times \mathcal{X} \rightarrow \mathbb{R}: (x, x') \mapsto k(x,x')\quad \quad \quad (1),

yaitu fungsi yang memetakan 2 buah objek dari ruang data \mathcal{X} ke dalam bilangan real \mathbb{R}. Semakin besar nilai outputnya berarti kedua objek tersebut semakin sama, i.e., kebalikan dari fungsi jarak.

Misalkan \mathcal{X} merupakan ruang vektor, salah satu contoh kernel adalah dot product

\displaystyle k(x,x) = \langle x, x' \rangle.

Namun demikian, ruang input \mathcal{X} tidak selalu berbentuk ruang vektor, bisa berbentuk apa saja asal masih berupa non-empty set. Poin saya adalah bisa jadi \mathcal{X} tidak ‘dikaruniai’ dengan dot product atau fungsi sejenis yang menyatakan kesamaan.

Untuk mengatasi hal tersebut, trik yang dapat digunakan adalah memetakan input ke dalam ruang yang pasti memiliki operasi dot product, yang kita namakan sebagai feature space, dengan fungsi

\displaystyle \Phi: \mathcal{X} \rightarrow \mathcal{H}: x \mapsto \Phi(x)\quad \quad \quad (2).

Fungsi ini juga biasa dikenal sebagai feature map. Maka, sekarang akan selalu aman apabila fungsi kernel didefinisikan di ruang \mathcal{H} dalam bentuk

k(x,x') = \langle \Phi(x), \Phi(x') \rangle_{\mathcal{H}} \quad \quad \quad (3)

karena \langle \cdot, \cdot \rangle_{\mathcal{H}} akan selalu ada.

Bahkan, walaupun ruang input \mathcal{X} berbentuk ruang vektor yang ‘dikarunai’ dot product, pemetaan (2) masih, di hampir semua kasus, diperlukan terutama pada problem klasifikasi atau regresi jika ruang input \mathcal{X} terlalu kompleks / non-linear untuk dikategorisasi. Salah contoh dari feature map adalah polynomial map. Misalkan \mathcal{X} \subseteq \mathbb{R}^2, polynomial map didefinisikan sebagai berikut

\displaystyle \Phi: \mathbb{R}^{2} \rightarrow \mathbb{R}^{3}: [x_1, x_2]^{\top} \mapsto [x_1^2, x_2^2, \sqrt{2} x_1 x_2]^{\top} \quad \quad \quad (4).

Dalam banyak kasus, trik pemetaan ini (terutama dalam bentuk fungsi non-linear) terbukti efektif. Namun demikian, pemetaan ini biasanya berakhir pada dimensionality explosion, i.e., dimensi feature space menjadi sangat besar, mungkin saja tak hingga, sehingga menambah kompleksitas komputasi.

Sebagai ilustrasi, untuk ruang input \mathcal{X} berdimensi d, i.e., \mathcal{X} \subseteq \mathbb{R}^d, polynomial map orde ke-2 pada persamaan (4) akan menjadi

\displaystyle \Phi(x) = [ x^2_d, ..., x^2_1, \sqrt{2} x_d x_{d-1}, ..., \sqrt{2} x_d x_1, \sqrt{2} x_{d-1} x_{d-2}, ..., \sqrt{2} x_{d-1} x_1,..., \sqrt{2} x_2 x_1 ] \quad \quad \quad (5).

Dapat dibayangkan laju pertambahan kompleksitas komputasi terhadap penambahan dimensi input yang tentunya memperberat pekerjaan komputer. Penggunaan kernel memungkinkan kita untuk dapat memanfaatkan mapping  seperti yang telah dijelaskan tanpa harus mengkalkulasi feature map secara eksplisit.

Pertanyaan pertama adalah apakah setiap kernel berasosiasi dengan feature map? Jawabannya adalah tidak. Hanya kernel yang memenuhi Mercer’s condition dan/atau positive definiteness (ekuivalen satu sama lain) yang memiliki pasangan feature map.

Teorema 1 (Mercer’s condition) Misalkan \mathcal{X} \subset \mathbb{R}^d sebuah himpunan dan k: \mathcal{X} \times \mathcal{X} \rightarrow \mathbb{R} merupakan sebuah fungsi kontinu dan simetris. k dapat diekspresikan dengan sebuah ekspansi konvergen dalam bentuk

\displaystyle k(x, y) = \sum_{n=0}^{\infty} a_n \phi_n(x) \phi_n(y)

dimana a_n > 0, jika dan hanya jika untuk setiap fungsi c: \mathcal{X} \rightarrow \mathbb{R}, c \in L_2(\mathcal{X}) memenuhi kondisi di bawah ini

\displaystyle \int \int_{\mathcal{X} \times \mathcal{X}} c(x) c(y) k(x,y)dx dy \geq 0.

Definisi 1 (Positive Definite Symmetric Kernel) Kernel k: \mathcal{X} \times \mathcal{X} \rightarrow \mathbb{R} disebut sebagai positive definite symmetric (PDS) kernel jika untuk setiap \{x_1, ...., x_N \} \subseteq \mathcal{X}, matriks \mathbf{K} = [K(x_i, x_j)]_{ij} \in \mathbb{R}^{N \times N}, yang disebut sebagai Gram matrix atau kernel matrix, bersifat simetris dan positive semi-definite.

Matriks \mathbf{K} merupakan matriks positive semi-definite  jika salah satu kondisi di bawah ini terpenuhi:

  • eigenvalues dari \mathbf{K} tidak negatif;
  • diberikan vektor \mathbf{c} \in \mathbb{R}^N,

\displaystyle \mathbf{c}^{\top} \mathbf{K} \mathbf{c} = \sum_{i,j=1}^{N} c_i c_j K(x_i, x_j) \geq 0.

2. CONTOH KERNEL

Berikut ini merupakan beberapa contoh dari fungsi kernel k: \mathcal{X} \times \mathcal{X} \rightarrow \mathbb{R} yang memenuhi Teorema 1 dan/atau Definisi 1, \forall x, x' \in \mathcal{X}

1. Polynomial Kernel

\displaystyle k(x,x') = (\langle x, x' \rangle + c)^m

Ini merupakan kernel yang berasosiasi dengan polynomial feature map pada persamaan (4), namun dihitung dalam orde ke-m.

2. Gaussian Kernel atau Radial Basis Function (RBF)

\displaystyle k(x,x') = \exp ( - \frac{\| x - x' \|}{2 \sigma^2})

dimana \sigma>0 merupakan parameter yang disebut sebagai bandwidth. Kernel ini paling sering digunakan di berbagai aplikasi.

3. Sigmoid Kernel

\displaystyle k(x,x') = \tanh ( a \langle x,x' \rangle + b )

3. LINEAR & RIDGE REGRESSION

Linear Regression

Sebelum menggunakan kernel, ada baiknya saya membahas terlebih dahulu salah satu algoritma paling dasar di statistics dan juga machine learning, yaitu Linear Regression (LR). Bahkan mungkin di SMA atau di tingkat awal perguruan tinggi LR sudah diajarkan. Selanjutnya kita akan masuk ke bagian bagaimana memodifikasi algoritma ini dengan menggunakan kernel.

Diberikan sebuah himpunan data latih berlabel \{ \mathbf{x}^{(i)}, y^{(i)}\}_{i=1}^{n} dimana \mathbf{x} \in \mathcal{X} \subseteq \mathbb{R}^d dan y \in \mathbb{R}, tujuan dari linear regression adalah untuk mencari fungsi predictor \hat{g}(\mathbf{x}') \approx y', dimana (\mathbf{x}', y') \not\in \{ \mathbf{x}^{(i)}, y^{(i)} \}_{i=1}^{n}, i.e., dari data yang tidak terdapat pada data latih.

Metode ini disebut linear  karena diawal kita mengasumsikan model linear f: \mathbb{R}^d \rightarrow \mathbb{R} yang dapat diekspresikan sebagai berikut:

\displaystyle f(\mathbf{x}) = \sum_{j=1}^{d} x_j \beta_j + \beta_0 \quad \quad \quad (6)

dimana \boldsymbol{\beta} = [ \beta_0, ..., \beta_d ]^{\top} merupakan parameter yang biasa disebut sebagai coefficients. Mencari model yang optimal berkorespondensi dengan menghitung nilai yang tepat untuk \boldsymbol{\beta}. Perhatikan bahwa persamaan (6) hanyalah merupakan persamaan garis, namun di dimensi d, dimana \{ \beta_j\}_{j=1}^{d} merupakan slope dan \beta_0 merupakan intercept (untuk menyederhanakan notasi, mulai sekarang saya akan meniadakan \beta_0 di persamaan-persamaan berikutnya).

Bagaimanakah menghitung koefisien \boldsymbol{\beta} \in \mathbb{R}^{d} ? Pertama-tama kita perlu mendefinisikan sebuah loss function yang menyatakan “seberapa dekat” fungsi f(\mathbf{x}) dengan label y. Untuk kasus ini, kita menggunakan Residual Sum of Square (RSS) yang didefinisikan sebagai berikut:

\displaystyle RSS^{lr}(\boldsymbol{\beta}) = \sum_{i=1}^{n} \left( \sum_{j=1}^{d} x^{(i)}_j \beta_j - y^{(i)} \right)^2 \quad \quad \quad (7)

Algoritma Linear Regression merupakan algoritma untuk mencari \boldsymbol{\beta} sehingga kita mendapatkan nilai minimum untuk RSS^{lr}(\boldsymbol{\beta}), i.e.,

\displaystyle \hat{\boldsymbol{\beta}} := \arg \min_{\boldsymbol{\beta}} RSS^{lr}(\beta) \quad \quad \quad (8)

Dengan kalkulus, parameter optimal \hat{\boldsymbol{\beta}} dapat dihitung dengan turunan \frac{\partial RSS^{lr}( \boldsymbol{\beta} )}{\partial \boldsymbol{\beta}} = 0. Sebelum menghitung turunan tsb, saya akan mengubah RSS^{lr}(\cdot) ke dalam bentuk matriks-vektor

\displaystyle RSS^{lr}(\boldsymbol{\beta}) = \| \mathbf{X} \boldsymbol{\beta} - \mathbf{y} \|^2_2 \quad \quad \quad (9)

dimana notasi \| \mathbf{v} \|_2 = \sqrt{ v_1 + ... + v_d} menyatakan panjang vektor, i.e., vector norm di ruang l_2 (Euclidean norm). Sekarang kita akan menghitung turunan pertama dari persamaan (9)

\displaystyle \frac{\partial RSS^{lr}(\boldsymbol{\beta}) }{\partial\boldsymbol{\beta}} = 2 \mathbf{X}^{\top} \left( \mathbf{X} \boldsymbol{\beta} - \mathbf{y} \right)=0 \quad \quad \quad (10)

maka

\displaystyle \hat{\boldsymbol{\beta}}^{lr} = (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{y} \quad \quad \quad (11)

dimana setiap baris ke-i dari matriks \mathbf{X} \in \mathbb{R}^{n \times d} merupakan data latih \mathbf{x}^{(i)}, dan setiap elemen ke-i dari vektor \mathbf{y} \in \mathbb{R}^{n} berisi label y^{(i)}. Untuk memahami dari mana datangnya persamaan (10) & (11), sedikit pemahaman tentang linear algebra dan vector calculus akan sangat membantu. Perlu diperhatikan bahwa kita akan mendapatkan solusi unik dari persamaan (11) jika matriks \mathbf{X}^{\top} \mathbf{X} \in \mathbb{R}^{d \times d}merupakan matrix full rank dan dapat di-inverse.

Selama permasalahan yang dihadapi memiliki karakteristik linear, i.e., hubungan antara \mathbf{x} dan y cukup linear, maka LR akan memberikan solusi yang reasonable, yaitu

\displaystyle g(\mathbf{x}') = \hat{\boldsymbol{\beta}}^{\top} \mathbf{x}' \approx y'

Ridge Regression

Masalah yang sering timbul pada Linear Regression adalah overfitting. Fungsi predictor g(\cdot) akan akurat pada data latih \mathbf{x} namun tidak akurat pada data tes \mathbf{x}'. Salah satu trik yang dapat dilakukan untuk mengurangi overfitting adalah dengan mencari parameter optimal \hat{\boldsymbol{\beta}}^{lr} dengan panjang yang sekecil mungkin. Pembahasan mengenai mengapa trik ini dapat menghasilkan model regresi yang lebih baik cukup panjang dan intimidating. Perlu pemahaman lebih lanjut mengenai bias-variance tradeoff. Untuk kasus ini, pada intinya kita menginginkan model yang memiliki low variance  (berkorespondensi dengan generalization) yang dapat dicapai dengan \hat{\boldsymbol{\beta}}^{ridge} yang lebih pendek.

Secara matematis, Linear Regression akan menjadi Ridge Regression (RR) dengan sedikit mengubah persamaan (9) menjadi

\displaystyle RSS^{rr}(\boldsymbol{\beta}) = \| \mathbf{X} \boldsymbol{\beta} - \mathbf{y} \|^2_2+ \lambda \| \boldsymbol{\beta}\|^2_2 \quad \quad \quad (12)

dimana \lambda > 0 merupakan konstanta yang mengontrol ‘seberapa banyak’ panjang vektor \boldsymbol{\beta} akan ‘dikecilkan’. Dengan cara yang sama seperti LR, nilai optimal \hat{\boldsymbol{\beta}} dari persamaan (12) didapatkan dengan persamaan

\displaystyle \hat{\boldsymbol{\beta}}^{rr} = (\mathbf{X}^{\top} \mathbf{X} + \lambda \mathbf{I}_d)^{-1} \mathbf{X}^{\top} \mathbf{y} \quad \quad \quad (13)

dimana \mathbf{I}_d merupakan matriks identitas berdimensi d \times d. Perhatikan bahwa matriks (\mathbf{X}^{\top} \mathbf{X} + \lambda \mathbf{I}_d) akan selalu dapat di-inverse. Perhatikan pula bahwa jika \lambda = 0 maka persamaan (13) akan sama dengan persamaan (11), yang dapat diartikan: LR merupakan kasus spesial dari RR.

Walaupun pada umumnya RR dapat menghasilkan predictor g(\cdot) yang lebih baik dibandingkan LR, kita harus melakukan tuning terhadap nilai \lambda > 0 yang bisa jadi berbeda-beda bergantung dengan data yang diobservasi. Proses ini seringkali dikenal dengan istilah hyper-parameter tuning.

4. Kernel Ridge Regression

Sekarang kita akan mengaplikasikan kernel pada Ridge Regression sehingga menghasilkan metode yang disebut sebagai Kernel Ridge Regression (KRR). Lebih spesifiknya, ini dilakukan dengan memodifikasi persamaan (13) sedemikian rupa sehingga kita dapat memasukkan kernel tanpa mengubah arti persamaannya.

Terlebih dahulu, dengan sedikit manipulasi aljabar matriks-vektor, perhatikan bahwa

\displaystyle (\mathbf{X}^{\top} \mathbf{X} + \lambda \mathbf{I}_d) \mathbf{X}^{\top} = \mathbf{X}^{\top} \mathbf{X} \mathbf{X}^{\top} + \lambda \mathbf{X}^{\top}=\mathbf{X}^{\top} (\mathbf{X} \mathbf{X}^{\top} + \lambda \mathbf{I}_n)

sehingga persamaan di bawah ini terpenuhi:

\displaystyle (\mathbf{X}^{\top} \mathbf{X} + \lambda \mathbf{I}_d)^{-1} \mathbf{X}^{\top}=\mathbf{X}^{\top} (\mathbf{X} \mathbf{X}^{\top} + \lambda \mathbf{I}_n)^{-1}.

Dengan demikian, kita akan mudah mengamati bahwa persamaan di bawah ini ekuivalen dengan persamaan (13)

\hat{ \boldsymbol{\beta} }^{krr} = \mathbf{X}^{\top} ( \mathbf{X} \mathbf{X}^{\top}+\lambda \mathbf{I}_n )^{-1} \mathbf{y} \quad \quad \quad (14)

Sampai saat ini, belum terlalu terlihat jelas bagaimana caranya kita dapat menggunakan kernel untuk menghitung nilai optimal \hat{\boldsymbol{\beta}}.

Ambil sebuah non-linear feature map \phi : \mathbb{R}^d \rightarrow \mathcal{H}. Jika \phi(\cdot) diaplikasikan pada \{ \mathbf{x}^{(i)}\}, maka kita dapat membuat matriks baru \mathbf{\Phi} \in \mathbb{R}^{n \times d}, dimana baris ke-i berisi vektor \phi(\mathbf{x}^{(i)}). Dengan demikian, kita dapat mengubah persamaan (14) menjadi sebagai berikut:

\hat{ \boldsymbol{\beta} }_{\mathcal{H}}^{krr} = \mathbf{\Phi}^{\top} ( \mathbf{\Phi} \mathbf{\Phi}^{\top}+\lambda \mathbf{I}_n )^{-1} \mathbf{y} \quad \quad \quad (15)

Perbedaan antara (14) dan (15) adalah (15) bekerja pada ruang feature \mathcal{H}, sedangkan (14) bekerja pada ruang data \mathbb{R}^d. Dengan menggunakan solusi (15), kita dapat menghitung fungsi prediksi KRR

g(\mathbf{x}') = \hat{\boldsymbol{\beta}}^{krr \top} \mathbf{x}'

g(\mathbf{x}') = \left( \mathbf{\Phi}^{\top} ( \mathbf{\Phi} \mathbf{\Phi}^{\top}+\lambda \mathbf{I}_n )^{-1} \mathbf{y} \right)^{\top} \mathbf{x}'

g(\mathbf{x}') = \mathbf{y}^{\top} \left( \mathbf{\Phi} \mathbf{\Phi}^{\top} + \lambda \mathbf{I}_n \right)^{-1} \mathbf{\Phi}\mathbf{x}' \quad \quad \quad (16)

Namun demikian, persamaan (16) masih memiliki masalah apabila dimensi dari ruang \mathcal{H} tidak terhingga, i.e., vektor \phi(\mathbf{x}) memiliki panjang tak terhingga, sehingga \mathbf{\Phi} \in \mathbb{R}^{\infty \times \infty} yang berarti tidak mungkin bisa dihitung oleh komputer !

Sekarang kita akan mensubstitusi beberapa elemen dengan menggunakan kernel. Pertama-tama, perhatikan bahwa matriks \mathbf{\Phi} \mathbf{\Phi}^{\top} \in \mathbb{R}^{n \times n} merupakan matriks Gram \mathbf{K} = [k( \mathbf{x}_i, \mathbf{x}_i) ]_{ij}, dimana kernel k didefinisikan seperti pada persamaan (3). Dan juga vektor \mathbf{\Phi} \mathbf{x}' \in \mathbb{R}^{n} dapat dibuat menjadi vektor kernel  \kappa(\mathbf{x}') = [ k(\mathbf{x}', \mathbf{x}_1) , ...,k(\mathbf{x}', \mathbf{x}_n))]^{\top}. Dengan demikian, bentuk final fungsi predictor dari KRR dapat dinyatakan dengan

g(\mathbf{x}') = \mathbf{y}^{\top} ( \mathbf{K} + \lambda \mathbf{I}_n)^{-1} \kappa(\mathbf{x}') \quad \quad \quad (17)

dimana komputasi eksplisit dari \phi: \mathbb{R}^{d} \rightarrow \mathcal{H} tidak lagi diperlukan.

5. Eksperimen pada Wine dataset

Saya melakukan eksperimen untuk membandingkan algoritma LR, RR, dan KRR. Dataset yang digunakan untuk mengevaluasi algoritma2 tsb adalah Wine Quality Datasets. Di sana terdapat 2 jenis datasets: red wine dan white wine. Dalam kesempatan ini saya hanya mengunakan white wine.

Data white wine terdiri dari 11 atribut yang mendeskripsikan karakteristik dari white wine, dan 4,898 instances. Kemudian dataset dibagi menjadi training set dan test set, dimana training set terdiri dari 90% sampel dari keseluruhan data yang diambil secara acak, dan 10% sisanya menjadi test set.

Tugas kita adalah untuk memprediksi kualitas dari White wine tersebut dengan skor berupa bilangan real berkisar antara 0 - 10, dimana 0: terburuk, 10: terbaik. Hasil prediksi tersebut (\hat{\mathbf{y}}) kemudian dibandingkan dengan skor sebenarnya yang didapatkan dari human expert (\mathbf{y}) dengan menggunakan Mean Square Error (MSE):

\displaystyle MSE(\hat{\mathbf{y}}, \mathbf{y}) = \frac{1}{n} \| \hat{\mathbf{y}} - \mathbf{y} \|_2^2

Evaluasi dilakukan 10 kali untuk setiap algoritma. Khusus untuk KRR, fungsi kernel yang digunakan adalah Gaussian kernel. Semua hyper-parameter , e.g., ridge penalty (\lambda) dan kernel bandwidth (\sigma), di-tuning dengan menggunakan grid search. Hasil evaluasi berupa rata-rata dan standar deviasi dari MSE untuk tiap algoritma setelah 10 kali eksekusi ditampilkan pada tabel di bawah ini:

LR

RR (\lambda=10)

KRR (\lambda=10, \sigma=1.4)

Training MSE

 0.716 \pm 0.004  0.718 \pm 0.004  0.590 \pm 0.004

Test MSE

  0.735 \pm 0.037  0.719 \pm 0.039  0.663 \pm 0.032

Dari tabel tersebut terlihat jelas bahwa KRR memberikan MSE terkecil dibandingkan 2 metode lainnya, yang berarti KRR merupakan metode yang terbaik untuk kasus ini. Jika kita membandingkan LR dan RR, RR lebih baik dari LR pada data tes namun tidak begitu signifikan. Pada data latih, RR bahkan performanya sedikit menurut bila dibandingkan LR.

Namun demikian, KRR membutuhkan waktu sekitar 21 detik untuk 1 kali eksekusi dimana LR dan RR hanya membutuhkan waktu kira 0.02 detik, yang berarti KRR ~1000 kali lebih lambat dibandingkan LR atau RR. Kondisi ini akan memburuk apabila data latih bertambah. Ini dikarenakan komputasi dari kernel pada KRR yang cukup time consuming.

Secara umum, memang algoritma yang memanfaatkan kernel tidak terlalu cocok untuk large scale problem. Saat ini sudah terdapat berbagai studi bagaimana mempercepat metode kernel tanpa banyak kehilangan performa prediksi. Beberapa cara diantaranya dengan menggunakan metode Nystrom  atau random Fourier features, yang kemungkinan menjadi pembahasan utama saya pada part berikutnya.

Source code

Saya melakukan eksperimen di atas dengan menggunakan MATLAB. Berikut code hasil implementasi dari LR, RR, dan KRR yang saya tulis.

Main file (main_regress.m)


clear; clc;
% main_regress.m
% - The main file for running either LR, RR, or KRR
% - Reading the White wine dataset downloaded from http://www3.dsi.uminho.pt/pcortez/wine/

%Read the white wine dataset
allData = dlmread('Wine/winequality-white.csv',';',1,0);
allX = allData(:,1:end-1);
N = size(allX,1);
Ns = int32(N*0.9); %number of training data
Nt = N - Ns; % number of test data
allY = allData(:,end);

MSE_train_list = [];
MSE_test_list = [];

r = 10;
for iter = 1:r
 % Get training dataset
 idx_s = randperm(N, Ns);
 Xtrain = allX(idx_s,:);
 ytrain = allY(idx_s);

% Get test dataset
 idx_t = setdiff([1:N], idx_s);
 Xtest = allX(idx_t,:);
 ytest = allY(idx_t);

% Pre-process the data into representations that have zero-mean and unit-variance
 Xtrain = zscore(Xtrain);
 Xtest = zscore(Xtest);
 ytrain = zscore(ytrain);
 ytest = zscore(ytest);

start_t = cputime;
 % % == MAIN ALGORITHM (choose one of them !) ==
 % % linear regression
 % [gtrain, gtest, beta] = Regressor.lrPredict(Xtrain, ytrain, Xtest);

% ridge regression
 lambda = 10; %tuned by grid search
 [gtrain, gtest, beta] = Regressor.rrPredict(Xtrain, ytrain, Xtest, lambda);

% % ridge regression
 % lambda = 10; %tuned by grid search
 % sigma = 1.4; %tuned by grid search
 % [gtrain, gtest] = Regressor.krrPredict(Xtrain, ytrain, Xtest, lambda,sigma);

% == END MAIN ALG ==
 elapsed_t = cputime - start_t;
 fprintf('Elapsed time (in seconds) : %f\n',elapsed_t);

MSE_train = mean( (gtrain - ytrain).^2 );
 MSE_test = mean( (gtest - ytest).^2 );

% recap over r runs
 MSE_train_list = [MSE_train_list; MSE_train];
 MSE_test_list = [MSE_test_list; MSE_test];

fprintf('Mean Squared Error (MSE) : %f (train), %f (test) \n', MSE_train, MSE_test);
end

%Compute mean and std over r runs
mean(MSE_train_list)
std(MSE_train_list)

mean(MSE_test_list)
std(MSE_test_list)


Regressor.m


% Regressor.m
% - The class containing functions that run either LR, RR, or KRR
% - The Gaussian kernel function is also included
classdef Regressor 
 methods (Static)
 function [gtrain, gtest, beta] = lrPredict(Xs,ys, Xt)
 % Linear Regression
 % Input : 
 % - Xs : training datapoint matrix with dimension n x (d +1)
 % - ys : training label vector with dimension n 
 % - Xt : test datapoint matrix with dimension n x (d + 1)
 % Output:
 % - beta : the optimal LR parameter vector
 % - gtrain : the optimal predicted training label vector
 % - gtest : the optimal predicted test label vector

% get the number of rows
 ns = size(Xs,1); 
 nt = size(Xt,1);

% add vector [1, 1, ..., 1]^{\top} to the first column of X 
 % to consider the "intercept" parameter
 Xs = [ones(ns,1) Xs]; 
 Xt = [ones(nt,1) Xt];


% compute the optimal coefficients 
 beta = inv(Xs'*Xs)*Xs'*ys;
 
 % compute the predicted labels
 gtrain = Xs * beta;
 gtest = Xt * beta;

end
 function [gtrain, gtest, beta] = rrPredict(Xs,ys, Xt, lambda)
 % Ridge Regression
 % Input : 
 % - Xs : training datapoint matrix with dimension n x (d +1)
 % - ys : training label vector with dimension n 
 % - Xt : test datapoint matrix with dimension n x (d + 1)
 % - lambda : ridge penalty 
 
 % Output:
 % - beta : the optimal RR parameter vector
 % - gtrain : the optimal predicted training label vector
 % - gtest : the optimal predicted test label vector

% get the number of rows
 ns = size(Xs,1); 
 nt = size(Xt,1);

% add vector [1, 1, ..., 1]^{\top} to the first column of X 
 % to consider the "intercept" parameter
 Xs = [ones(ns,1) Xs]; 
 Xt = [ones(nt,1) Xt];


% compute the optimal coefficients
 d = size(Xs,2); % get the data dimension
 Id = eye(d,d); % create an identity matrix with dimension d x d
 beta = inv(Xs'*Xs + lambda*Id)*Xs'*ys;
 
 % compute the predicted labels
 gtrain = Xs * beta;
 gtest = Xt * beta;

end
 function [gtrain, gtest] = krrPredict(Xs,ys, Xt, lambda, sigma)
 % Kernel Ridge Regression
 % Input : 
 % - Xs : training datapoint matrix with dimension n x (d +1)
 % - ys : training label vector with dimension n 
 % - Xt : test datapoint matrix with dimension n x (d + 1)
 % - lambda: ridge penalty
 % - sigma : Gaussian kernel bandwidth
 % Output:
 % - gtrain : the optimal predicted training label vector
 % - gtest : the optimal predicted test label vector

% get the number of rows
 ns = size(Xs,1); 
 nt = size(Xt,1);

% add vector [1, 1, ..., 1]^{\top} to the first column of X 
 % to consider the "intercept" parameter
 Xs = [ones(ns,1) Xs]; 
 Xt = [ones(nt,1) Xt];

% compute the training kernel matrix and training-test kernel matrix
 Kss = Regressor.gaussKernel(Xs',Xs',sigma);
 Kst = Regressor.gaussKernel(Xs',Xt',sigma);


% compute the predicted labels
 In = eye(ns,ns); % create an identity matrix with dimension ns x ns


 gtrain = ys'*inv(Kss + lambda * In)*Kss;
 gtrain = gtrain';
 gtest = ys'*inv(Kss + lambda * In)*Kst;
 gtest = gtest';
 end
 function K = gaussKernel(X,Y, sigma)
 % X : d x NcolX
 % Y : d x NcolY
 Xsq = sum(X.^2);
 Ysq = sum(Y.^2);
 NcolX = size(X,2);
 NcolY = size(Y,2);
 Dsq = (ones(NcolY,1)*Xsq)' + ones(NcolX,1)*Ysq - 2*X'*Y;
 K = exp(-Dsq/(2*sigma.^2));
 end
 end
end

2 thoughts on “Fast Domain Generalization with Kernel Methods – Part 1

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s