Numerical Power in Regression: Sweep Operator & Cholesky in Action
Analisis regresi linier merupakan salah satu teknik fundamental dalam statistika yang digunakan untuk memodelkan hubungan antara variabel respon dan satu atau lebih variabel penjelas. Pendekatan ini tidak hanya penting secara teoretis, tetapi juga menjadi dasar bagi berbagai metode analisis data modern seperti machine learning regression, regularization, dan generalized linear models. Secara matematis, model regresi dapat diuraikan dengan pendekatan aljabar matriks, yang menjadikannya sangat bergantung pada operasi komputasi matriks seperti transpose, perkalian, dan invers.
Namun, pada praktiknya, perhitungan invers matriks terutama pada ukuran besar sering kali menghadapi dua kendala utama:
Untuk mengatasi permasalahan tersebut, dikembangkan sejumlah algoritma komputasi yang efisien dan stabil secara numerik, di antaranya:
Pendekatan Sweep Operator dan Cholesky Decomposition merupakan dua teknik inti dalam komputasi regresi modern, terutama pada perangkat lunak statistik seperti R, di mana keduanya menjadi bagian dari prosedur internal untuk menyelesaikan persamaan normal pada regresi linier.
Regresi linier merupakan metode paling mendasar dalam analisis hubungan antar variabel. Tujuannya adalah untuk memodelkan bagaimana variabel respon (y) dipengaruhi oleh satu atau lebih variabel penjelas (\(x_1, x_2, \dots, x_p\)) melalui persamaan linier.
Model regresi linier sederhana menggambarkan hubungan antara satu variabel respon \(y\) dan satu variabel penjelas \(x\) dalam bentuk:
\[ y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \quad i = 1, 2, \dots, n \]
dengan:
Pendugaan parameter \(\beta_0\) dan \(\beta_1\) dilakukan menggunakan metode kuadrat terkecil, yaitu dengan meminimalkan jumlah kuadrat galat:
\[ SSE = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 = \sum_{i=1}^{n} (y_i - \beta_0 - \beta_1 x_i)^2 \]
Turunan pertama dari fungsi tersebut terhadap \(\beta_0\) dan \(\beta_1\) menghasilkan sistem persamaan normal:
\[ n\beta_0 + \beta_1\sum x_i = \sum y_i \\ \beta_0\sum x_i + \beta_1\sum x_i^2 = \sum x_i y_i \]
yang solusinya memberikan: \[ \hat{\beta}_1 = \frac{\sum (x_i - \bar{x})(y_i - \bar{y})}{\sum (x_i - \bar{x})^2}, \quad \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x} \]
Model regresi linier berganda memperluas konsep sebelumnya dengan melibatkan lebih dari satu variabel penjelas. Secara umum dapat dituliskan sebagai:
\[ y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \dots + \beta_{p-1} x_{i,p-1} + \varepsilon_i, \quad i = 1, 2, \dots, n \]
dengan: - \(y_i\): nilai variabel respon untuk pengamatan ke-\(i\) - \(x_{ij}\): nilai variabel penjelas ke-\(j\) untuk pengamatan ke-\(i\) - \(\beta_j\): koefisien regresi yang menunjukkan perubahan rata-rata \(y\) untuk setiap kenaikan satu satuan pada \(x_j\) ketika variabel lain tetap - \(\varepsilon_i\): komponen galat acak yang mengandung semua faktor selain variabel penjelas
Secara prinsip, metode pendugaan parameter yang digunakan tetap metode kuadrat terkecil, yaitu dengan meminimalkan:
\[ SSE = \sum_{i=1}^{n} (y_i - \beta_0 - \beta_1 x_{i1} - \beta_2 x_{i2} - \dots - \beta_{p-1} x_{i,p-1})^2 \]
Turunan pertama terhadap setiap \(\beta_j\) (\(j = 0, 1, 2, \dots, p-1\)) akan menghasilkan \(p\) persamaan normal yang saling berkaitan.
Pada kasus berganda, sistem ini menjadi lebih kompleks, sehingga penulisan dan penyelesaiannya jauh lebih efisien dilakukan dalam bentuk matriks, yaitu:
\[ \mathbf{X}'\mathbf{X}\boldsymbol{\beta} = \mathbf{X}'\mathbf{y} \]
Model regresi linier berganda dapat ditulis secara ringkas dan elegan menggunakan notasi matriks. Pendekatan ini mempermudah analisis matematis dan implementasi komputasi, terutama saat jumlah pengamatan (\(n\)) dan variabel penjelas (\(p-1\)) cukup besar.
Secara umum, model regresi linier dalam bentuk matriks dinyatakan sebagai:
\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \]
dengan: - \(\mathbf{y}_{n\times1}\) : vektor respon - \(\mathbf{X}_{n\times p}\) : matriks desain (design matrix) - \(\boldsymbol{\beta}_{p\times1}\) : vektor parameter (koefisien regresi) - \(\boldsymbol{\varepsilon}_{n\times1}\) : vektor galat acak
Matriks \(\mathbf{X}\) memiliki bentuk umum:
\[ \mathbf{X} = \begin{bmatrix} 1 & X_{1,1} & X_{1,2} & \dots & X_{1,p-1} \\ 1 & X_{2,1} & X_{2,2} & \dots & X_{2,p-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & X_{n,1} & X_{n,2} & \dots & X_{n,p-1} \end{bmatrix} \]
Setiap baris \(\mathbf{x}_i'\) dari matriks \(\mathbf{X}\) merepresentasikan satu pengamatan, sedangkan setiap kolom (selain kolom pertama yang berisi 1) merepresentasikan satu variabel penjelas.
Parameter regresi \(\boldsymbol{\beta}\) ditaksir dengan metode kuadrat terkecil dengan meminimalkan:
\[ SSE = (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})'(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \]
Menurunkan fungsi ini terhadap \(\boldsymbol{\beta}\) dan menyamakannya dengan nol menghasilkan persamaan normal:
\[ \mathbf{X}'\mathbf{X}\boldsymbol{\beta} = \mathbf{X}'\mathbf{y} \]
Solusi dari persamaan normal diperoleh dengan:
\[ \hat{\boldsymbol{\beta}} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y} \]
Pendugaan ini hanya dapat dilakukan apabila \(\mathbf{X}'\mathbf{X}\) berpangkat penuh (full rank), yaitu kolom-kolom \(\mathbf{X}\) tidak saling bergantung linier.
Persamaan normal menjadi inti dari seluruh metode regresi linier. Sebagian besar algoritma komputasi—termasuk QR decomposition, Cholesky decomposition, dan Sweep Operator—pada dasarnya merupakan cara-cara berbeda untuk menyelesaikan sistem persamaan:
\[ \mathbf{X}'\mathbf{X}\boldsymbol{\beta} = \mathbf{X}'\mathbf{y} \]
# Contoh data regresi berganda
x1 <- c(1, 2, 3, 4, 5)
x2 <- c(2, 1, 4, 3, 5)
y <- c(2.3, 2.7, 3.8, 3.5, 5.1)
# Membentuk matriks X dan vektor y
X <- cbind(1, x1, x2)
y <- matrix(y, ncol = 1)
# Menghitung (X'X)^(-1) X'y secara manual
beta_hat <- solve(t(X) %*% X) %*% t(X) %*% y
beta_hat [,1]
1.3633333
x1 0.3777778
x2 0.3277778Hasil di atas akan menghasilkan vektor \(\hat{\boldsymbol{\beta}} = [\hat{\beta}_0, \hat{\beta}_1, \hat{\beta}_2]'\).
Sebelum membahas metode komputasi regresi seperti Sweep Operator dan Dekomposisi Cholesky, penting untuk memahami sifat-sifat matematis dari matriks yang muncul dalam persamaan normal, yaitu \(\mathbf{X}'\mathbf{X}\).
Sebuah matriks persegi \(A\) berukuran \(n \times n\) dikatakan simetrik apabila elemen-elemen di atas dan di bawah diagonal utamanya identik, yaitu:
\[ A = A' \quad \text{atau secara elemen: } a_{ij} = a_{ji}, \ \forall i, j \]
Contoh:
\[ A = \begin{bmatrix} 4 & 2 & 1 \\ 2 & 5 & 3 \\ 1 & 3 & 6 \end{bmatrix} \quad \Rightarrow \quad A' = A \]
Dalam konteks regresi, matriks \(\mathbf{X}'\mathbf{X}\) selalu simetrik, karena:
\[ (\mathbf{X}'\mathbf{X})' = \mathbf{X}'(\mathbf{X}' )' = \mathbf{X}'\mathbf{X} \]
Selain simetrik, \(\mathbf{X}'\mathbf{X}\) juga bersifat definit positif, artinya untuk setiap vektor tak nol \(\mathbf{z}\) berlaku:
\[ \mathbf{z}'(\mathbf{X}'\mathbf{X})\mathbf{z} > 0 \]
Sifat ini berarti bentuk kuadrat \(\mathbf{z}'A\mathbf{z}\) menghasilkan nilai positif untuk semua \(\mathbf{z} \neq 0\). Dengan kata lain, \(A\) (dalam konteks ini \(\mathbf{X}'\mathbf{X}\)) memiliki semua eigenvalue positif, sehingga matriks tersebut dapat diinvers dan stabil digunakan dalam perhitungan OLS.
# Membentuk matriks X dan menghitung X'X
x1 <- c(1, 2, 3, 4, 5)
x2 <- c(2, 1, 4, 3, 5)
X <- cbind(1, x1, x2)
XtX <- t(X) %*% X
XtX x1 x2
5 15 15
x1 15 55 53
x2 15 53 55Perhatikan bahwa elemen di atas dan di bawah diagonal utama identik — menunjukkan simetri.
# 1. Cek apakah simetrik
isSymmetric(XtX)[1] TRUE# 2. Periksa eigenvalue
eigen(XtX)$values[1] 112.1978456 2.0000000 0.8021544# 3. Uji positif definit
all(eigen(XtX)$values > 0)[1] TRUEJika semua nilai eigen bernilai positif (TRUE), maka matriks \(\mathbf{X}'\mathbf{X}\) dikatakan definit positif, sehingga aman digunakan dalam penyelesaian persamaan normal maupun untuk dekomposisi Cholesky.
Sweep Operator merupakan prosedur komputasi berbasis eliminasi Gauss yang digunakan untuk menghitung penduga parameter regresi tanpa harus melakukan invers matriks secara eksplisit.
Metode ini banyak digunakan dalam perangkat lunak statistik seperti SAS dan R karena memiliki keunggulan dari sisi efisiensi dan stabilitas numerik.
Dalam regresi linier, kita perlu menyelesaikan persamaan normal:
\[ \mathbf{X}'\mathbf{X}\boldsymbol{\beta} = \mathbf{X}'\mathbf{y} \]
Untuk menyelesaikannya tanpa menghitung \((\mathbf{X}'\mathbf{X})^{-1}\), kita dapat menggunakan prosedur sweep terhadap matriks gabungan:
\[ B = \begin{bmatrix} \mathbf{X}'\mathbf{X} & \mathbf{X}'\mathbf{y} \\ \mathbf{y}'\mathbf{X} & \mathbf{y}'\mathbf{y} \end{bmatrix} \]
Operasi sweep dilakukan terhadap elemen diagonal ke-\(k\) dari matriks simetrik \(A = [a_{ij}]\) untuk mengeliminasi kolom dan baris ke-\(k\), sehingga menghasilkan bentuk matriks baru yang berisi informasi estimasi koefisien regresi dan varians residual.
Notasi Singkat
Untuk setiap elemen diagonal ke-\(k\) yang akan di-sweep:
Proses ini dapat dilakukan secara bertahap dan reversible, artinya hasil sweep terhadap indeks tertentu dapat dikembalikan ke bentuk semula dengan melakukan sweep ulang pada indeks yang sama.
Pada contoh ini, kita tidak memulai dari data mentah, melainkan dari ringkasan matriks hasil penghitungan:
\[ \mathbf{X}'\mathbf{X}= \begin{bmatrix} 6 & 0 & 12\\ 0 & 6 & 0\\ 12 & 0 & 28 \end{bmatrix},\qquad \mathbf{X}'\mathbf{y}= \begin{bmatrix} 12\\ 2\\ 25 \end{bmatrix},\qquad \mathbf{y}'\mathbf{y}=28. \]
Asumsi urutan kolom \(\mathbf{X}\): \([\,\text{intersep},\,x_1,\,x_2\,]\).
Sebelum melakukan SWEEP, bentuk matriks gabungan:
\[ B= \begin{bmatrix} \mathbf{X}'\mathbf{X} & \mathbf{X}'\mathbf{y}\\ \mathbf{y}'\mathbf{X} & \mathbf{y}'\mathbf{y} \end{bmatrix} = \begin{bmatrix} 6 & 0 & 12 & 12\\ 0 & 6 & 0 & 2 \\ 12& 0 & 28 & 25\\ 12& 2 & 25 & 28 \end{bmatrix}. \]
Tujuan SWEEP: melakukan operasi eliminasi Gauss khusus untuk matriks simetrik agar kita memperoleh koefisien () (di kolom terakhir, kecuali baris terakhir) dan JKG/SSE (elemen kanan-bawah).
Di bawah ini kita tunjukkan praktik SWEEP langkah demi langkah. Notasi:
Pivot \(D=a_{11}=6\).
Bagi baris 1 dengan 6: \[ [6,\ 0,\ 12,\ 12] / 6 \;\Rightarrow\; [1,\ 0,\ 2,\ 2]. \]
Eliminasi baris lain:
Baris 2: \(B=a_{21}=0\)
\(\text{baris}_2 \leftarrow \text{baris}_2 - 0 \times \text{baris}_1\) (tidak berubah).
Set \(a_{21}=-0/6=0\).
Baris 3: \(B=a_{31}=12\)
\([12,0,28,25] - 12\times[1,0,2,2] = [0,0,4,1]\).
Set \(a_{31}=-12/6=-2\).
Baris 4: \(B=a_{41}=12\)
\([12,2,25,28] - 12\times[1,0,2,2] = [0,2,1,4]\).
Set \(a_{41}=-12/6=-2\).
Hasil setelah SWEEP(A, k=1): \[ A^{(1)}= \begin{bmatrix} 0.1667 & 0 & 2 & 2\\ 0 & 6 & 0 & 2\\ -2 & 0 & 4 & 1\\ -2 & 2 & 1 & 4 \end{bmatrix}. \]
Interpretasi (model intersep saja \(y=b_0\$):\)b_0 = A^{(1)}_{1,4} = 2$, \(\text{SSE} = A^{(1)}_{4,4} = 4\).
Gunakan \(A^{(1)}\) sebagai input.
Pivot \(D=a_{22}=6\).
Bagi baris 2 dengan 6: \[ [0,6,0,2]/6 \Rightarrow [0,1,0,0.3333]. \]
Eliminasi baris lain:
Hasil setelah SWEEP(A, k=1:2): \[ A^{(2)}= \begin{bmatrix} 0.1667 & 0 & 2 & 2\\ 0 & 0.1667 & 0 & 0.3333\\ -2 & 0 & 4 & 1\\ -2 & -0.3333& 1 & 3.3333 \end{bmatrix}. \]
Interpretasi (model \(y=b_0+b_1x_1\)):
\(\hat b_0 = 2.0,\ \hat b_1 = 0.3333,\ \text{SSE} = 3.3333\).
Lakukan prosedur yang sama pada pivot \(a_{33}\) dari \(A^{(2)}\).
Hasil akhirnya $ b_0 = 1.5,b_1 = 0.3333,b_2 = 0.25, . $
sweep.operator() (R)library(fastmatrix)
XX <- matrix(c(6,0,12,
0,6,0,
12,0,28), nrow=3, byrow=TRUE)
Xy <- matrix(c(12,2,25), ncol=1)
yy <- 28
A <- rbind(cbind(XX, Xy), cbind(t(Xy), yy))
# (1) k = 1 -> intersep saja
res_k1 <- sweep.operator(A, k = 1)
b0 <- res_k1[1,4]
sse0 <- res_k1[4,4]
# (2) k = 1:2 -> intersep + x1
res_k12 <- sweep.operator(A, k = 1:2)
b0_b1 <- res_k12[1:2,4]
sse1 <- res_k12[4,4]
# (3) k = 1:3 -> intersep + x1 + x2
res_k123 <- sweep.operator(A, k = 1:3)
coef_all <- res_k123[1:3,4]
sse2 <- res_k123[4,4]
list(
k1_intercept_only = c(b0 = b0, SSE = sse0),
k12_b0_b1 = c(b0 = b0_b1[1], b1 = b0_b1[2], SSE = sse1),
k123_full = c(b0 = coef_all[1], b1 = coef_all[2], b2 = coef_all[3], SSE = sse2)
)$k1_intercept_only
b0.yy SSE.yy
2 4
$k12_b0_b1
b0 b1 SSE.yy
2.0000000 0.3333333 3.3333333
$k123_full
b0 b1 b2 SSE.yy
1.5000000 0.3333333 0.2500000 3.0833333 Output utama:
| Tahap | Persamaan Model | Koefisien (urut: \(b_0, b_1, b_2\)) | JKG/SSE |
|---|---|---|---|
| 1 | \(y = b_0\) | 2.00 | 4.000 |
| 2 | \(y = b_0 + b_1x_1\) | 2.00, 0.3333 | 3.333 |
| 3 | \(y = b_0 + b_1x_1 + b_2x_2\) | 1.50, 0.3333, 0.2500 | 3.083 |
Dekomposisi Cholesky merupakan salah satu metode numerik yang sangat efisien untuk menyelesaikan sistem persamaan linear dari regresi linier, yaitu:
\[ \mathbf{X}'\mathbf{X}\boldsymbol{\beta} = \mathbf{X}'\mathbf{y} \]
Metode ini tidak menghitung invers matriks secara langsung, tetapi memfaktorkan matriks simetrik positif definit \(\mathbf{X}'\mathbf{X}\) menjadi hasil kali dua matriks segitiga:
\[ \mathbf{X}'\mathbf{X} = \mathbf{L}\mathbf{L}' \]
dengan: - \(\mathbf{L}\) : matriks segitiga bawah (lower triangular matrix), - \(\mathbf{L}'\) : transpose dari \(\mathbf{L}\).
Faktorisasi ini mempermudah penyelesaian karena sistem dapat diselesaikan dalam dua tahap substitusi: 1. \(\mathbf{L}\mathbf{z} = \mathbf{X}'\mathbf{y}\) (forward substitution) 2. \(\mathbf{L}'\boldsymbol{\beta} = \mathbf{z}\) (backward substitution)
Sehingga diperoleh vektor parameter \(\boldsymbol{\beta}\) tanpa menghitung invers matriks.
Misalkan matriks \(\mathbf{X}'\mathbf{X}\) berukuran \(3 \times 3\):
\[ \mathbf{X}'\mathbf{X} = \begin{bmatrix} 6 & 0 & 12 \\ 0 & 6 & 0 \\ 12 & 0 & 28 \end{bmatrix} \]
Karena \(\mathbf{X}'\mathbf{X}\) bersifat simetrik dan positif definit, maka dapat difaktorkan menjadi:
\[ \mathbf{L} = \begin{bmatrix} l_{11} & 0 & 0 \\ l_{21} & l_{22} & 0 \\ l_{31} & l_{32} & l_{33} \end{bmatrix} \quad \text{dengan} \quad \mathbf{X}'\mathbf{X} = \mathbf{L}\mathbf{L}'. \]
Langkah-langkahnya adalah:
\[ \begin{aligned} l_{11} &= \sqrt{a_{11}} = \sqrt{6} = 2.4495,\\ l_{21} &= \frac{a_{21}}{l_{11}} = 0,\\ l_{31} &= \frac{a_{31}}{l_{11}} = \frac{12}{2.4495} = 4.8990,\\[4pt] l_{22} &= \sqrt{a_{22} - l_{21}^2} = \sqrt{6 - 0^2} = 2.4495,\\ l_{32} &= \frac{a_{32} - l_{31}l_{21}}{l_{22}} = \frac{0 - 4.8990(0)}{2.4495} = 0,\\[4pt] l_{33} &= \sqrt{a_{33} - l_{31}^2 - l_{32}^2} = \sqrt{28 - 4.8990^2 - 0^2} = \sqrt{4.01} = 2.0025. \end{aligned} \]
Maka diperoleh: \[ \mathbf{L} = \begin{bmatrix} 2.4495 & 0 & 0\\ 0 & 2.4495 & 0\\ 4.8990 & 0 & 2.0025 \end{bmatrix}. \]
Kita tahu bahwa: \[ \mathbf{L}\mathbf{L}'\boldsymbol{\beta} = \mathbf{X}'\mathbf{y} \] Dengan data: \[ \mathbf{X}'\mathbf{y} = \begin{bmatrix} 12\\ 2\\ 25 \end{bmatrix}. \]
Langkah 1 (forward substitution):
\(\mathbf{L}\mathbf{z} = \mathbf{X}'\mathbf{y}\)
Hitung \(\mathbf{z}\) secara bertahap.
\[ \begin{aligned} 2.4495z_1 &= 12 \quad \Rightarrow \quad z_1 = 4.899,\\ 2.4495z_2 &= 2 \quad \Rightarrow \quad z_2 = 0.816,\\ 4.8990z_1 + 2.0025z_3 &= 25 \quad \Rightarrow \quad 2.0025z_3 = 25 - (4.8990)(4.899) = 25 - 24.01 = 0.99 \\ &\Rightarrow z_3 = 0.495. \end{aligned} \]
Maka diperoleh: \[ \mathbf{z} = \begin{bmatrix} 4.899\\ 0.816\\ 0.495 \end{bmatrix}. \]
Langkah 2 (backward substitution):
\(\mathbf{L}'\boldsymbol{\beta} = \mathbf{z}\)
Tuliskan sistemnya: \[ \begin{bmatrix} 2.4495 & 0 & 4.8990\\ 0 & 2.4495 & 0\\ 0 & 0 & 2.0025 \end{bmatrix} \begin{bmatrix} b_0\\b_1\\b_2 \end{bmatrix} = \begin{bmatrix} 4.899\\0.816\\0.495 \end{bmatrix}. \]
Dari bawah ke atas: \[ \begin{aligned} 2.0025b_2 &= 0.495 \Rightarrow b_2 = 0.247,\\ 2.4495b_1 &= 0.816 \Rightarrow b_1 = 0.333,\\ 2.4495b_0 + 4.899(0.247) &= 4.899 \Rightarrow b_0 = 1.500. \end{aligned} \]
Hasil akhir: \[ \hat{\boldsymbol{\beta}} = \begin{bmatrix} 1.500\\ 0.333\\ 0.250 \end{bmatrix}. \]
Koefisien tersebut identik dengan hasil dari Sweep Operator dan solusi OLS biasa.
# Matriks simetrik positif definit (X'X)
XtX <- matrix(c(6,0,12,
0,6,0,
12,0,28),
nrow = 3, byrow = TRUE)
# Vektor X'y
Xy <- matrix(c(12,2,25), ncol = 1)
# Dekomposisi Cholesky
L <- chol(XtX) # default menghasilkan upper triangular
L [,1] [,2] [,3]
[1,] 2.44949 0.00000 4.898979
[2,] 0.00000 2.44949 0.000000
[3,] 0.00000 0.00000 2.000000# Karena R menghasilkan L', ambil transpose
Lt <- t(L)
# Forward substitution: Lt %*% z = X'y
z <- forwardsolve(Lt, Xy)
# Backward substitution: L %*% beta = z
beta_hat <- backsolve(L, z)
beta_hat [,1]
[1,] 1.5000000
[2,] 0.3333333
[3,] 0.2500000