Package

library(quadprog)

Quadratic Programming

Quadratic programming (QP) adalah bidang optimasi matematis yang bertujuan untuk menemukan nilai minimum dari fungsi kuadrat tertentu, yang terbatas oleh sejumlah batasan linear atau kuadrat.Secara umum, masalah quadratic programming dapat dirumuskan sebagai berikut:

Meminimumkan \(-d^{T}b+\frac{1}{2}b^{T}Db\) dengan kendala \(A^{T}b\ge b_{0}\) (bisa ada/tidak)

Ket:

  1. b: Vektor keputusan dari inner product
  2. D: Matriks simetris dari quadratic form
  3. d: Vektor koefisien

Tanpa Kendala

Kasus: \(-8x_1-16x_2+x_1^2+4x_2^2\)

Ingat kembali inner product dan quadratic form dari video kuliah (https://youtu.be/VTk2krvzfG0) untuk mengonversi persamaan menjadi bentuk formulasi QP.

Hint:

Matriks A dan vektor b0 beranggotakan 0 semua (tidak ada kendala)

D<-matrix(c(2,0,0,8),2,2,byrow=T);D
##      [,1] [,2]
## [1,]    2    0
## [2,]    0    8
d<-c(8,16);d
## [1]  8 16
A<-matrix(c(0,0,0,0),2,2,byrow=T);A #tanpa kendala
##      [,1] [,2]
## [1,]    0    0
## [2,]    0    0
b0<-c(0,0) #tanpa kendala
solve.QP(D,d,A,b0)
## $solution
## [1] 4 2
## 
## $value
## [1] -32
## 
## $unconstrained.solution
## [1] 4 2
## 
## $iterations
## [1] 1 0
## 
## $Lagrangian
## [1] 0 0
## 
## $iact
## [1] 0

Dengan Kendala \((\ge)\)

Kasus: \(-8x_1-16x_2+x_1^2+4x_2^2\) dengan kendala \(x1+x2\ge5, x_1\ge4.5\)

Hint:

  • Matriks A berukuran (j.kendala x j. variabel), isikan nilai kendala ke-1 di baris ke-1 dst.

  • Jika salah satu kendala hanya berisi 1 variabel, maka wakili variabel lainnya dengan nilai 0. 

A<-matrix(c(1,1,1,0),2,2,byrow=T);A #kendala >=
##      [,1] [,2]
## [1,]    1    1
## [2,]    1    0
b0<-c(5,4.5);b0 #kendala >=
## [1] 5.0 4.5
solve.QP(D,d,A,b0)
## $solution
## [1] 4.5 2.0
## 
## $value
## [1] -31.75
## 
## $unconstrained.solution
## [1] 4 2
## 
## $iterations
## [1] 2 0
## 
## $Lagrangian
## [1] 0 1
## 
## $iact
## [1] 2

Dengan Kendala \((\le)\)

Kasus: \(-8x_1-16x_2+x_1^2+4x_2^2\) dengan kendala \(x1+x2\le5, x_1\le3\)

Hint:

  • \(a\le b\) sama saja dengan \(-a\ge -b\). Maka, kendalanya dapat dituliskan: \(-x1-x2\ge-5, -x_1\ge-3\) 

A<-matrix(c(-1,-1,-1,0),2,2,byrow=T);A #kendala <=
##      [,1] [,2]
## [1,]   -1   -1
## [2,]   -1    0
b0<-c(-5,-3);b0 #kendala <=
## [1] -5 -3
solve.QP(D,d,A,b0)
## $solution
## [1] 3 2
## 
## $value
## [1] -31
## 
## $unconstrained.solution
## [1] 4 2
## 
## $iterations
## [1] 2 0
## 
## $Lagrangian
## [1] 0 2
## 
## $iact
## [1] 2

OLS

Mencari š›½ yang meminimumkan \(RSS = šœ–^Tšœ– =(y-X\beta)^T(y-X\beta)=y^Ty-2y^TX\beta+\beta^T X^TX\beta\)

Yang diminumkan adalah \(-2y^TX\beta+\beta^T X^TX\beta\)

Jika diforumulasikan berdasarkan kriteria QP sebelumnya, maka:

\(min(-d^{T}b+\frac{1}{2}b^{T}Db)=min(-(2y^TX)b+\frac{1}{2}b^{T}(2X^TX)b)\), sehingga \(D=2(X^TX)\) dan \(d^T=2y^TX\)

QP Manual

Siapkan matriks X dan Y (studi kasus: mtcars)

y<-as.matrix(mtcars$mpg);y
##       [,1]
##  [1,] 21.0
##  [2,] 21.0
##  [3,] 22.8
##  [4,] 21.4
##  [5,] 18.7
##  [6,] 18.1
##  [7,] 14.3
##  [8,] 24.4
##  [9,] 22.8
## [10,] 19.2
## [11,] 17.8
## [12,] 16.4
## [13,] 17.3
## [14,] 15.2
## [15,] 10.4
## [16,] 10.4
## [17,] 14.7
## [18,] 32.4
## [19,] 30.4
## [20,] 33.9
## [21,] 21.5
## [22,] 15.5
## [23,] 15.2
## [24,] 13.3
## [25,] 19.2
## [26,] 27.3
## [27,] 26.0
## [28,] 30.4
## [29,] 15.8
## [30,] 19.7
## [31,] 15.0
## [32,] 21.4
x<-as.matrix(mtcars[,2:7]);x
##                     cyl  disp  hp drat    wt  qsec
## Mazda RX4             6 160.0 110 3.90 2.620 16.46
## Mazda RX4 Wag         6 160.0 110 3.90 2.875 17.02
## Datsun 710            4 108.0  93 3.85 2.320 18.61
## Hornet 4 Drive        6 258.0 110 3.08 3.215 19.44
## Hornet Sportabout     8 360.0 175 3.15 3.440 17.02
## Valiant               6 225.0 105 2.76 3.460 20.22
## Duster 360            8 360.0 245 3.21 3.570 15.84
## Merc 240D             4 146.7  62 3.69 3.190 20.00
## Merc 230              4 140.8  95 3.92 3.150 22.90
## Merc 280              6 167.6 123 3.92 3.440 18.30
## Merc 280C             6 167.6 123 3.92 3.440 18.90
## Merc 450SE            8 275.8 180 3.07 4.070 17.40
## Merc 450SL            8 275.8 180 3.07 3.730 17.60
## Merc 450SLC           8 275.8 180 3.07 3.780 18.00
## Cadillac Fleetwood    8 472.0 205 2.93 5.250 17.98
## Lincoln Continental   8 460.0 215 3.00 5.424 17.82
## Chrysler Imperial     8 440.0 230 3.23 5.345 17.42
## Fiat 128              4  78.7  66 4.08 2.200 19.47
## Honda Civic           4  75.7  52 4.93 1.615 18.52
## Toyota Corolla        4  71.1  65 4.22 1.835 19.90
## Toyota Corona         4 120.1  97 3.70 2.465 20.01
## Dodge Challenger      8 318.0 150 2.76 3.520 16.87
## AMC Javelin           8 304.0 150 3.15 3.435 17.30
## Camaro Z28            8 350.0 245 3.73 3.840 15.41
## Pontiac Firebird      8 400.0 175 3.08 3.845 17.05
## Fiat X1-9             4  79.0  66 4.08 1.935 18.90
## Porsche 914-2         4 120.3  91 4.43 2.140 16.70
## Lotus Europa          4  95.1 113 3.77 1.513 16.90
## Ford Pantera L        8 351.0 264 4.22 3.170 14.50
## Ferrari Dino          6 145.0 175 3.62 2.770 15.50
## Maserati Bora         8 301.0 335 3.54 3.570 14.60
## Volvo 142E            4 121.0 109 4.11 2.780 18.60
X<-cbind(rep(1,nrow(x)),x);X
##                       cyl  disp  hp drat    wt  qsec
## Mazda RX4           1   6 160.0 110 3.90 2.620 16.46
## Mazda RX4 Wag       1   6 160.0 110 3.90 2.875 17.02
## Datsun 710          1   4 108.0  93 3.85 2.320 18.61
## Hornet 4 Drive      1   6 258.0 110 3.08 3.215 19.44
## Hornet Sportabout   1   8 360.0 175 3.15 3.440 17.02
## Valiant             1   6 225.0 105 2.76 3.460 20.22
## Duster 360          1   8 360.0 245 3.21 3.570 15.84
## Merc 240D           1   4 146.7  62 3.69 3.190 20.00
## Merc 230            1   4 140.8  95 3.92 3.150 22.90
## Merc 280            1   6 167.6 123 3.92 3.440 18.30
## Merc 280C           1   6 167.6 123 3.92 3.440 18.90
## Merc 450SE          1   8 275.8 180 3.07 4.070 17.40
## Merc 450SL          1   8 275.8 180 3.07 3.730 17.60
## Merc 450SLC         1   8 275.8 180 3.07 3.780 18.00
## Cadillac Fleetwood  1   8 472.0 205 2.93 5.250 17.98
## Lincoln Continental 1   8 460.0 215 3.00 5.424 17.82
## Chrysler Imperial   1   8 440.0 230 3.23 5.345 17.42
## Fiat 128            1   4  78.7  66 4.08 2.200 19.47
## Honda Civic         1   4  75.7  52 4.93 1.615 18.52
## Toyota Corolla      1   4  71.1  65 4.22 1.835 19.90
## Toyota Corona       1   4 120.1  97 3.70 2.465 20.01
## Dodge Challenger    1   8 318.0 150 2.76 3.520 16.87
## AMC Javelin         1   8 304.0 150 3.15 3.435 17.30
## Camaro Z28          1   8 350.0 245 3.73 3.840 15.41
## Pontiac Firebird    1   8 400.0 175 3.08 3.845 17.05
## Fiat X1-9           1   4  79.0  66 4.08 1.935 18.90
## Porsche 914-2       1   4 120.3  91 4.43 2.140 16.70
## Lotus Europa        1   4  95.1 113 3.77 1.513 16.90
## Ford Pantera L      1   8 351.0 264 4.22 3.170 14.50
## Ferrari Dino        1   6 145.0 175 3.62 2.770 15.50
## Maserati Bora       1   8 301.0 335 3.54 3.570 14.60
## Volvo 142E          1   4 121.0 109 4.11 2.780 18.60

Buat matriks D, d, A, dan b

D<-2*t(X)%*%X;D
##                       cyl       disp         hp       drat         wt
##         64.000    396.000   14766.20    9388.00   230.1800   205.9040
## cyl    396.000   2648.000  103744.80   64408.00  1382.8000  1358.8080
## disp 14766.200 103744.800 4359254.94 2582728.80 50189.5920 54182.9776
## hp    9388.000  64408.000 2582728.80 1668556.00 32744.5600 32943.4880
## drat   230.180   1382.800   50189.59   32744.56   845.5814   717.4379
## wt     205.904   1358.808   54182.98   32943.49   717.4379   721.8021
## qsec  1142.320   6951.120  257603.01  162184.32  4113.8280  3656.1892
##            qsec
##        1142.320
## cyl    6951.120
## disp 257603.008
## hp   162184.320
## drat   4113.828
## wt     3656.189
## qsec  20586.960
d<-2*t(y)%*%X;d
##                cyl     disp       hp     drat       wt     qsec
## [1,] 1285.8 7387.2 257410.2 168725.4 4760.554 3819.506 23229.49
A<-t(matrix(0,dim(X)[2],dim(X)[2]));A
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]    0    0    0    0    0    0    0
## [2,]    0    0    0    0    0    0    0
## [3,]    0    0    0    0    0    0    0
## [4,]    0    0    0    0    0    0    0
## [5,]    0    0    0    0    0    0    0
## [6,]    0    0    0    0    0    0    0
## [7,]    0    0    0    0    0    0    0
b<-rep(0,dim(X)[2]);b
## [1] 0 0 0 0 0 0 0

Solusi

solve.QP(D,d,A,b)
## $solution
## [1] 26.30735899 -0.81856023  0.01320490 -0.01792993  1.32040573 -4.19083238
## [7]  0.40146117
## 
## $value
## [1] -13878.83
## 
## $unconstrained.solution
## [1] 26.30735899 -0.81856023  0.01320490 -0.01792993  1.32040573 -4.19083238
## [7]  0.40146117
## 
## $iterations
## [1] 1 0
## 
## $Lagrangian
## [1] 0 0 0 0 0 0 0
## 
## $iact
## [1] 0

Hasil lm()

ols<-lm(y~x);ols$coefficients 
## (Intercept)        xcyl       xdisp         xhp       xdrat         xwt 
## 26.30735899 -0.81856023  0.01320490 -0.01792993  1.32040573 -4.19083238 
##       xqsec 
##  0.40146117

QP manual dengan fungsi lm() menunjukkan hasil yang sama. Hal ini menjelaskan bagaimana QP digunakan dalam algoritma lm() bekerja untuk menghasilkan koefisien regresi.

LASSO

Merupakan OLS dengan kendala \(\sum_{k = 1}^{p} |\beta_p|<s\). Artinya, akan terdapat \(2^p\) kendala. Contoh untuk \(p=2\), dan \(s=1\) kombinasi yang memenuhi \(|\beta_1|+|\beta_2|\le 1\):

  • \(\beta_1+\beta_2\le 1\)

  • \(-\beta_1+\beta_2\le 1\)

  • \(\beta_1-\beta_2\le 1\)

  • \(-\beta_1-\beta_2\le 1\)

Maka, untuk kasus dengan \(p=2\), digunakan:

A

matrix(c(-1,-1,1,-1,-1,1,1,1),4,2,byrow=T)
##      [,1] [,2]
## [1,]   -1   -1
## [2,]    1   -1
## [3,]   -1    1
## [4,]    1    1

b0

rep(-1,4)
## [1] -1 -1 -1 -1

Penerapan Lasso untuk kasus mtcars sebelumnya

Model yang digunakan adalah model tanpa intersep. Berikut adalah dua pendekatan yang umum digunakan dalam penggunaan Lasso terkait dengan intersep:

  1. Menyertakan Intersep: Dalam pendekatan ini, kita menyertakan intersep sebagai variabel tambahan yang memiliki koefisien sendiri dalam model regresi. Ini berarti kita memiliki satu parameter tambahan yang harus diestimasi oleh Lasso.

  2. Tanpa Intersep: Dalam pendekatan ini, kita tidak menyertakan intersep dalam model regresi. Ini berarti bahwa kita mengasumsikan bahwa garis regresi harus melewati titik (0,0). Dalam hal ini, faktor yang berada di luar model dianggap tidak penting, artinya tidak ada parameter tambahan yang mewakili intersep yang perlu diestimasi.

Kali ini yang akan digunakan adalah skenario 2.

Pada kasus sebelumnya, terdapat 6 peubah X yang digunakan, \(p=6\), maka akan terdapat \(2^6=64\) kendala. Berikut penerapan elemen A dan b0 yang digunakan

ols.tb<-lm(y~x-1)
A<-as.matrix(expand.grid(rep(list(c(-1,1)),6)));A
##       Var1 Var2 Var3 Var4 Var5 Var6
##  [1,]   -1   -1   -1   -1   -1   -1
##  [2,]    1   -1   -1   -1   -1   -1
##  [3,]   -1    1   -1   -1   -1   -1
##  [4,]    1    1   -1   -1   -1   -1
##  [5,]   -1   -1    1   -1   -1   -1
##  [6,]    1   -1    1   -1   -1   -1
##  [7,]   -1    1    1   -1   -1   -1
##  [8,]    1    1    1   -1   -1   -1
##  [9,]   -1   -1   -1    1   -1   -1
## [10,]    1   -1   -1    1   -1   -1
## [11,]   -1    1   -1    1   -1   -1
## [12,]    1    1   -1    1   -1   -1
## [13,]   -1   -1    1    1   -1   -1
## [14,]    1   -1    1    1   -1   -1
## [15,]   -1    1    1    1   -1   -1
## [16,]    1    1    1    1   -1   -1
## [17,]   -1   -1   -1   -1    1   -1
## [18,]    1   -1   -1   -1    1   -1
## [19,]   -1    1   -1   -1    1   -1
## [20,]    1    1   -1   -1    1   -1
## [21,]   -1   -1    1   -1    1   -1
## [22,]    1   -1    1   -1    1   -1
## [23,]   -1    1    1   -1    1   -1
## [24,]    1    1    1   -1    1   -1
## [25,]   -1   -1   -1    1    1   -1
## [26,]    1   -1   -1    1    1   -1
## [27,]   -1    1   -1    1    1   -1
## [28,]    1    1   -1    1    1   -1
## [29,]   -1   -1    1    1    1   -1
## [30,]    1   -1    1    1    1   -1
## [31,]   -1    1    1    1    1   -1
## [32,]    1    1    1    1    1   -1
## [33,]   -1   -1   -1   -1   -1    1
## [34,]    1   -1   -1   -1   -1    1
## [35,]   -1    1   -1   -1   -1    1
## [36,]    1    1   -1   -1   -1    1
## [37,]   -1   -1    1   -1   -1    1
## [38,]    1   -1    1   -1   -1    1
## [39,]   -1    1    1   -1   -1    1
## [40,]    1    1    1   -1   -1    1
## [41,]   -1   -1   -1    1   -1    1
## [42,]    1   -1   -1    1   -1    1
## [43,]   -1    1   -1    1   -1    1
## [44,]    1    1   -1    1   -1    1
## [45,]   -1   -1    1    1   -1    1
## [46,]    1   -1    1    1   -1    1
## [47,]   -1    1    1    1   -1    1
## [48,]    1    1    1    1   -1    1
## [49,]   -1   -1   -1   -1    1    1
## [50,]    1   -1   -1   -1    1    1
## [51,]   -1    1   -1   -1    1    1
## [52,]    1    1   -1   -1    1    1
## [53,]   -1   -1    1   -1    1    1
## [54,]    1   -1    1   -1    1    1
## [55,]   -1    1    1   -1    1    1
## [56,]    1    1    1   -1    1    1
## [57,]   -1   -1   -1    1    1    1
## [58,]    1   -1   -1    1    1    1
## [59,]   -1    1   -1    1    1    1
## [60,]    1    1   -1    1    1    1
## [61,]   -1   -1    1    1    1    1
## [62,]    1   -1    1    1    1    1
## [63,]   -1    1    1    1    1    1
## [64,]    1    1    1    1    1    1
sigma.beta<-sum(abs(ols.tb$coefficients))
s<-0.5*sigma.beta #misalnya saja
b0<-rep(-s,2^6)

X dan Y yang digunakan merupakan adalah X dan Y hasil standarisasi

sdata<-as.data.frame(scale(data.frame(y,x)))
yscale<-as.matrix(sdata$y);yscale
##              [,1]
##  [1,]  0.15088482
##  [2,]  0.15088482
##  [3,]  0.44954345
##  [4,]  0.21725341
##  [5,] -0.23073453
##  [6,] -0.33028740
##  [7,] -0.96078893
##  [8,]  0.71501778
##  [9,]  0.44954345
## [10,] -0.14777380
## [11,] -0.38006384
## [12,] -0.61235388
## [13,] -0.46302456
## [14,] -0.81145962
## [15,] -1.60788262
## [16,] -1.60788262
## [17,] -0.89442035
## [18,]  2.04238943
## [19,]  1.71054652
## [20,]  2.29127162
## [21,]  0.23384555
## [22,] -0.76168319
## [23,] -0.81145962
## [24,] -1.12671039
## [25,] -0.14777380
## [26,]  1.19619000
## [27,]  0.98049211
## [28,]  1.71054652
## [29,] -0.71190675
## [30,] -0.06481307
## [31,] -0.84464392
## [32,]  0.21725341
xscale<-as.matrix(sdata[,2:7])
D<-2*t(xscale)%*%xscale
d<-2*t(yscale)%*%xscale

Solusi

solve.QP(D,d,t(A),b0)$solution
## [1] -0.2425580  0.2715466 -0.2039718  0.1171394 -0.6803694  0.1190301

Latihan

  1. Coba QP menggunakan model lasso dengan intersep.
  2. Coba bandingkan dengan fungsi glmnet() apakah hasilnya berbeda? Jika ya, apa yang menyebabkan hal itu? Parameter apa yang belum diakomodasi dari QP manual?