Pengantar
Dalam praktikum Optimisasi Statistika ini, kita akan mengeksplorasi
berbagai metode penyelesaian sistem persamaan linear (SPL). Kita akan
menggunakan pendekatan Eliminasi Gauss, Gauss-Jordan, Iterasi Jacobi,
dan Gauss-Seidel. Kedua metode eliminasi (Gauss dan Gauss-Jordan)
bertujuan untuk mengubah matriks SPL menjadi bentuk yang lebih sederhana
untuk menemukan solusi, sementara metode iterasi (Jacobi dan
Gauss-Seidel) mengandalkan pendekatan berulang untuk mendekati solusi.
Studi kasus pada data trees akan menunjukkan penerapan
metode-metode ini dalam konteks regresi linear. Dengan pendekatan
praktis ini, kita akan menggali lebih dalam ke dalam teori dan
implementasi praktis dari metode-metode tersebut dalam pemecahan masalah
nyata dalam statistika.
Eliminasi Gauss
Poin utama yang dilakukan pada Eliminasi Gauss adalah:
- Transformasi ke Bentuk Echelon: Eliminasi Gauss mengubah matriks menjadi bentuk echelon. Ini berarti matriks akan memiliki bentuk di mana semua elemen di bawah diagonal utama adalah nol.
- Back Substitution: Setelah matriks berada dalam bentuk echelon, kita lakukan back substitution untuk menemukan solusi persamaan.
Misalkan kita memiliki sistem persamaan linear:
\[ \begin{align*} 3x + y - z &= 5 \\ 4x + 7y - 3z &= 20 \\ 2x - 2y + 5z &= 10 \end{align*} \] Matriks augmented dari sistem ini adalah:
\[ \begin{pmatrix} 3 & 1 & -1 & | & 5 \\ 4 & 7 & -3 & | & 20 \\ 2 & -2 & 5 & | & 10 \end{pmatrix} \]
Implementasi di R
Berikut adalah fungsi Eliminasi Gauss secara general:
# Membuat matriks A (3x3) dan vector b
A <- matrix(c(3, 1, -1,
4, 7, -3,
2, -2, 5), nrow = 3, byrow = TRUE)
b <- c(5, 20, 10)# Fungsi Eliminasi Gauss
gauss_elimination <- function(A,b) {
n <- nrow(A)
A <- cbind(A,b)
for (i in 1:(n-1)) {
for (j in (i+1):n) {
factor <- A[j,i] / A[i,i]
A[j,] <- A[j,] - factor * A[i,]
}
}
return(A)
}Implementasi pada matriks A, didapat:
# Melakukan Eliminasi Gauss
A_gauss <- gauss_elimination(A,b)
A_gauss## b
## [1,] 3 1.000000 -1.000000 5.00000
## [2,] 0 5.666667 -1.666667 13.33333
## [3,] 0 0.000000 4.882353 12.94118Terbentuk matriks A_gauss dengan bentuk matriks segitiga atas.
Kemudian, untuk menemukan solusi, dilakukan back-substitution
# Fungsi Back Substitution
back_substitution <- function(A) {
n <- nrow(A)
x <- numeric(n)
x[n] <- A[n, n+1] / A[n, n]
for (i in (n-1):1) {
x[i] <- (A[i, n+1] - sum(A[i, (i+1):n] * x[(i+1):n])) / A[i, i]
}
return(x)
}
# Menemukan Solusi
solusi_gauss <- back_substitution(A_gauss)
print(solusi_gauss)## [1] 1.506024 3.132530 2.650602Eliminasi Gauss-Jordan
- Transformasi ke Bentuk Echelon Baris Tereduksi: Gauss-Jordan mengambil langkah lebih lanjut dari Eliminasi Gauss dengan tidak hanya membuat elemen di bawah diagonal utama menjadi nol, tetapi juga membuat semua elemen di atas diagonal utama menjadi nol. Akibatnya, matriks yang dihasilkan adalah bentuk echelon baris tereduksi.
- Langsung Dapat Solusi: Setelah matriks berada dalam bentuk ini, solusi langsung dapat dibaca dari matriks tanpa perlu back substitution.
Kita akan menggunakan contoh kasus sebelumnya, yaitu misalkan kita memiliki sistem persamaan linear:
\[ \begin{align*} 3x + y - z &= 5 \\ 4x + 7y - 3z &= 20 \\ 2x - 2y + 5z &= 10 \end{align*} \]
Implementasi di R
gauss_jordan_elimination <- function(A,b) {
n <- nrow(A)
A <- cbind(A,b)
for (k in 1:n) {
A[k,] <- A[k,] / A[k,k]
for (i in 1:n) {
if (i != k) {
factor <- A[i, k] / A[k, k]
A[i,] <- A[i,] - factor * A[k,]
}
}
}
return(A)
}
# Melakukan Eliminasi Gauss-Jordan
A_gauss_jordan <- gauss_jordan_elimination(A,b)
A_gauss_jordan## b
## [1,] 1 0 0 1.506024
## [2,] 0 1 0 3.132530
## [3,] 0 0 1 2.650602Sehingga, bentuk solusi dari x,y,z adalah
A_gauss_jordan[,4]## [1] 1.506024 3.132530 2.650602Iterasi Jacobi
Metode Iterasi Jacobi adalah teknik untuk menemukan solusi dari sistem persamaan linear yang memiliki matriks diagonal dominan. Algoritma ini memisahkan matriks koefisien A menjadi diagonal dan sisa bagian, yaitu D (diagonal) dan R (sisanya), sehingga A = D + R.
Misalkan sistem persamaan linear adalah Ax = b, dengan x sebagai vektor solusi yang ingin ditemukan. Pada setiap iterasi k+1, metode Jacobi menghitung solusi vektor x^{(k+1)} sebagai berikut:
\[ x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j \neq i} a_{ij} x_j^{(k)} \right), \quad i = 1, 2, ..., n \]
dimana \(a_{ij}\) adalah elemen dari matriks A pada baris i dan kolom j, dan \(b_i\) adalah elemen ke-i dari vektor b. Pada setiap iterasi, nilai \(x_i^{(k)}\) dari iterasi sebelumnya digunakan untuk menghitung \(x_i^{(k+1)}\).
Menggunakan sistem persamaan yang diberikan:
\[ \begin{align*} 3x + y - z &= 5 \\ 4x + 7y - 3z &= 20 \\ 2x - 2y + 5z &= 10 \end{align*} \]
Matriks augmented-nya adalah:
\[ \begin{pmatrix} 3 & 1 & -1 & | & 5 \\ 4 & 7 & -3 & | & 20 \\ 2 & -2 & 5 & | & 10 \end{pmatrix} \]
kita akan mengisolasi setiap variabel \(x, y, z\) pada setiap persamaan. Berikut ini adalah formulasi Jacobi yang benar untuk sistem persamaan baru ini:
\[ \begin{align*} x^{(k+1)} &= \frac{1}{3} \left( 5 - y^{(k)} + z^{(k)} \right) \\ y^{(k+1)} &= \frac{1}{7} \left( 20 - 4x^{(k)} + 3z^{(k)} \right) \\ z^{(k+1)} &= \frac{1}{5} \left( 10 - 2x^{(k)} + 2y^{(k)} \right) \end{align*} \]
Dalam persamaan ini, \(x^{(k+1)}, y^{(k+1)}, z^{(k+1)}\) adalah nilai-nilai pada iterasi ke-\((k+1)\), dan \(x^{(k)}, y^{(k)}, z^{(k)}\) adalah nilai-nilai dari iterasi sebelumnya, ke-\(k\). Pada setiap iterasi, nilai-nilai ini diperbarui secara simultan berdasarkan nilai-nilai dari iterasi sebelumnya, sesuai dengan metode Jacobi.
Implementasi di R
# Metode Jacobi
jacobi <- function(A, b, tol = 1e-6, maxiter = 1000) {
n <- length(b)
x <- rep(0, n) # Menginisialisasi solusi dengan 0
for (iter in 1:maxiter) {
x_new <- rep(0, n)
for (i in 1:n) {
s <- sum(A[i, -i] * x[-i])
x_new[i] <- (b[i] - s) / A[i, i]
}
if (all(abs(x - x_new) < tol)) {
return(list(x = x_new, iter = iter))
}
x <- x_new
}
return(list(x = x, iter = maxiter))
}
# Menjalankan metode Jacobi
jacobi_result <- jacobi(A, b)
cat("Jacobi Method: Solusi x y z =", jacobi_result$x, "setelah", jacobi_result$iter, "iterasi.\n")## Jacobi Method: Solusi x y z = 1.506024 3.13253 2.650602 setelah 22 iterasi.Dengan toleransi sebesar 1e-6, didapat solusi di atas.
Iterasi Gauss-Seidel
Metode Gauss-Seidel mirip dengan Jacobi, namun perbedaannya terletak pada pembaruan nilai solusi. Pada metode Gauss-Seidel, saat kita menghitung komponen i dari vektor solusi, kita menggunakan komponen yang sudah diperbarui pada iterasi saat itu jika tersedia.
Untuk sistem persamaan Ax = b, komponen i dari vektor solusi pada iterasi k+1 dihitung menggunakan:
\[ x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j = 1}^{i-1} a_{ij} x_j^{(k+1)} - \sum_{j = i+1}^{n} a_{ij} x_j^{(k)} \right), \quad i = 1, 2, ..., n \]
Dalam kasus ini, kita menggunakan nilai \(x_j^{(k+1)}\) untuk \(j < i\) yang sudah diperbarui pada iterasi saat itu, dan nilai \(x_j^{(k)}\) untuk \(j > i\) dari iterasi sebelumnya.
Untuk menyesuaikan persamaan Gauss-Seidel dengan matriks sebelumnya:
\[ \begin{pmatrix} 3 & 1 & -1 & | & 5 \\ 4 & 7 & -3 & | & 20 \\ 2 & -2 & 5 & | & 10 \end{pmatrix} \]
kita akan mengikuti aturan yang sama seperti sebelumnya, tetapi dengan menggunakan nilai terbaru dari variabel yang diperbarui selama iterasi saat itu. Berikut adalah formulasi yang benar untuk metode Gauss-Seidel pada sistem persamaan baru:
\[ \begin{align*} x^{(k+1)} &= \frac{1}{3} \left( 5 - y^{(k)} + z^{(k)} \right) \\ y^{(k+1)} &= \frac{1}{7} \left( 20 - 4x^{(k+1)} + 3z^{(k)} \right) \\ z^{(k+1)} &= \frac{1}{5} \left( 10 - 2x^{(k+1)} + 2y^{(k+1)} \right) \end{align*} \]
Dalam proses iterasi Gauss-Seidel:
Hitung \(x^{(k+1)}\) menggunakan nilai \(y^{(k)}\) dan \(z^{(k)}\) dari iterasi sebelumnya.
Segera setelah \(x^{(k+1)}\) dihitung, gunakan nilai ini untuk menghitung \(y^{(k+1)}\).
Kemudian, gunakan \(x^{(k+1)}\) dan \(y^{(k+1)}\) yang baru dihitung untuk menghitung \(z^{(k+1)}\).
Langkah-langkah ini diulang hingga solusi konvergen, yaitu ketika perbedaan antara solusi pada iterasi berturut-turut di bawah suatu toleransi yang telah ditentukan.
Perbedaan kunci antara Metode Iterasi Jacobi dan Gauss-Seidel adalah bagaimana nilai-nilai terbaru digunakan. Pada Metode Jacobi, semua nilai dari iterasi sebelumnya digunakan untuk menghitung nilai baru, sedangkan pada Metode Gauss-Seidel, nilai terbaru digunakan segera setelah tersedia, yang dapat menyebabkan konvergensi lebih cepat. Namun, kedua metode ini bergantung pada sifat matriks sistem persamaan linear dan tidak selalu konvergen untuk setiap sistem.
Implementasi di R
# Metode Gauss-Seidel
gaussSeidel <- function(A, b, epsilon, maxIterations) {
# Extract the size of the system
n <- length(b)
# Initialize the solution vector x
x <- numeric(n)
# Iterate until either epsilon is reached or maxIterations is reached
for (iteration in 1:maxIterations) {
x_old <- x
# Update each variable using the Gauss-Seidel formula
for (i in 1:n) {
x[i] <- (b[i] - sum(A[i, -i] * x[-i])) / A[i, i]
}
# Check if epsilon is reached
if (max(abs(x - x_old)) < epsilon) {
break
}
}
# Return the solution and the number of iterations
list(x = x, iterations = iteration)
}
# Menjalankan metode Gauss-Seidel
gauss_seidel_result <- gaussSeidel(A, b, 1e-6, 1000)
cat("Gauss-Seidel Method: Solusi x y z =", gauss_seidel_result$x, "setelah", gauss_seidel_result$iterations, "iterasi.\n")## Gauss-Seidel Method: Solusi x y z = 1.506024 3.13253 2.650602 setelah 15 iterasi.Studi Kasus (Regresi Linear)
Diberikan studi kasus yang bertujuan untuk membuat program R yang
dapat menyelesaikan/mencari koefisien beta bagi permasalahan pada data
trees yang terdapat pada database di R dengan menggunakan
ke-4 metode yang telah digunakan.
Preprocess Data
Input Data Trees
Terdapat 3 peubah dalam data trees, yaitu
Girth: Diameter pohonHeight: Tinggi pohonVolume: Volume pohon (dalam ft kubik)
Ringkasan dari data trees adalah sebagai berikut:
head(trees)str(trees)## 'data.frame': 31 obs. of 3 variables:
## $ Girth : num 8.3 8.6 8.8 10.5 10.7 10.8 11 11 11.1 11.2 ...
## $ Height: num 70 65 63 72 81 83 66 75 80 75 ...
## $ Volume: num 10.3 10.3 10.2 16.4 18.8 19.7 15.6 18.2 22.6 19.9 ...Plot
GGally::ggpairs(trees)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
Dapat terlihat bahwa terdapat hubungan linear antar peubah pada data
trees. Dalam hal ini, kita akan perlakukan Volume sebagai
variabel respons (Y).
Bentuk data dalam sebuah matrix A
pred <- cbind(intercept=1, Girth=trees$Girth, Height=trees$Height)
head(pred)## intercept Girth Height
## [1,] 1 8.3 70
## [2,] 1 8.6 65
## [3,] 1 8.8 63
## [4,] 1 10.5 72
## [5,] 1 10.7 81
## [6,] 1 10.8 83Bentuk Vector b
resp <- trees$Volume
head(resp)## [1] 10.3 10.3 10.2 16.4 18.8 19.7Transformasi:
A <- t(pred) %*% pred #(X'X)
b <- t(pred) %*% resp #(X'y)
Ab <- cbind(A,b)Regresi Linear Sederhana
Sebagai validasi untuk setiap metode yang akan digunakan, kita akan menggunakan base function dari regresi linear sebagai acuan.
reglin <- lm(Volume~.,data=trees)
coef_reglin <- reglin$coefficients
coef_reglin## (Intercept) Girth Height
## -57.9876589 4.7081605 0.3392512Nilai di atas berperan sebagai validasi untuk metode selanjutnya.
Eliminasi Gauss
A_Gauss <- gauss_elimination(A,b)
A_Gauss## intercept Girth Height
## intercept 31 410.7000 2356.0000 935.3000
## Girth 0 295.4374 311.5000 1496.6435
## Height 0 0.0000 889.5641 301.7857coef_gauss <- back_substitution(A_Gauss)
rbind(coef_gauss, coef_reglin)## (Intercept) Girth Height
## coef_gauss -57.98766 4.708161 0.3392512
## coef_reglin -57.98766 4.708161 0.3392512Dapat terlihat bahwa nilai coef gaus dan coef reglin memiliki nilai yang sama persis.
Eliminasi Gauss-Jordan
A_Gauss_jordan <- gauss_jordan_elimination(A,b)
A_Gauss_jordan## intercept Girth Height
## intercept 1 0 0 -57.9876589
## Girth 0 1 0 4.7081605
## Height 0 0 1 0.3392512coef_gauss_jordan <- A_Gauss_jordan[,4]
rbind(coef_gauss_jordan, coef_reglin)## intercept Girth Height
## coef_gauss_jordan -57.98766 4.708161 0.3392512
## coef_reglin -57.98766 4.708161 0.3392512Dapat terlihat bahwa nilai coef gaus jordan dan coef reglin memiliki nilai yang sama persis.
Gauss Seidel
gauss_seidel_result <- gaussSeidel(A, b, 1e-6, 1000)
gauss_seidel_result## $x
## [1] -57.7752385 4.7123396 0.3357443
##
## $iterations
## [1] 1000coef_gauss_seidel <- gauss_seidel_result$x
rbind(coef_gauss_jordan, coef_reglin)## intercept Girth Height
## coef_gauss_jordan -57.98766 4.708161 0.3392512
## coef_reglin -57.98766 4.708161 0.3392512Dapat terlihat bahwa nilai coef gaus jordan dan coef reglin memiliki nilai yang sama persis.
Iterasi Jacobi
# Menjalankan metode Jacobi
jacobi_result <- jacobi(A, b, maxiter=10000)
jacobi_result## $x
## [1] -Inf -Inf -Inf
##
## $iter
## [1] 10000Dari hasil di atas dapat terlihat bahwa setelah 10000 kali iterasi,
nilai \(b_0\), \(b_1\), dan \(b_2\) konvergen ke \(-inf\). Artinya, metode yang digunakan
tidak tepat untuk menganalisis data trees. Untuk lebih
jelasnya, perhatikan pengerjaan iterasi jacobi menggunakan excel pada
gambar di bawah ini:
Dapat terlihat bahwa pada iterasi ke-10 dan seterusnya, nilai error atau selisih antara nilai \(x_{k+1}\) dengan \(x_k\) konstan di angka 151%. Artinya, nilai dari \(x_{k+1}\) (x yang baru) akan selalu 2.5 kali lebih besar daripada nilai \(x_k\) (x sebelumnya). Hal inilah yang membuat solusi yang didapat konvergen ke \(-inf\).
Lantas mengapa nilai konvergensi x terhadap nol tidak dapat ditemukan?
Terdapat beberapa kemungkinan yang dapat terjadi. Namun, hal yang paling mungkin adalah karena diagonal utama pada matrix A (yang dijadikan nilai pembagi untuk setiap iterasi) memiliki rentang cukup jauh seperti yang terlihat pada nilai di bawah:
A## intercept Girth Height
## intercept 31.0 410.70 2356.0
## Girth 410.7 5736.55 31524.7
## Height 2356.0 31524.70 180274.0D <- diag(diag(A))
D## [,1] [,2] [,3]
## [1,] 31 0.00 0
## [2,] 0 5736.55 0
## [3,] 0 0.00 180274Dapat terlihat bahwa nilai untuk pembentuk intersep (\(b_0\)) sangatlah kecil jika dibandingkan dengan nilai yang digunakan untuk membentuk koefisien Height (\(b_1\)). Sehingga, didapatkan invers matriks D:
solve(D)## [,1] [,2] [,3]
## [1,] 0.03225806 0.0000000000 0.000000e+00
## [2,] 0.00000000 0.0001743208 0.000000e+00
## [3,] 0.00000000 0.0000000000 5.547112e-06yang juga memiliki rentang nilai cukup jauh. Nilai \(5.54e-06 \approx 1\) akan membuat matriks \(D^{-1}\) seperti matriks singular (melanggar asumsi metode iterasi jacobi).
Alasan lain yang mungkin (menurut ChatGPT) adalah:
Diagonal Dominance: Agar metode iteratif Jacobi konvergen, matriks harus memiliki sifat dominansi diagonal, artinya elemen diagonal harus lebih besar dari jumlah mutlak dari elemen-elemen lainnya di baris yang sama. Jika matriks tidak memiliki sifat ini, metode Jacobi mungkin tidak akan konvergen.
Pemilihan Nilai Awal (x0): Jika nilai awal yang dipilih jauh dari solusi sebenarnya, metode Jacobi mungkin memerlukan lebih banyak iterasi untuk konvergen atau bahkan tidak akan konvergen sama sekali.
Error Numerik: Pada kasus tertentu, error numerik akibat operasi floating-point yang berulang-ulang dapat menyebabkan divergensi dari metode iteratif.
Matriks Singular atau Hampir Singular: Jika matriks yang digunakan dalam metode Jacobi adalah singular atau hampir singular (determinan mendekati nol), maka metode iteratif mungkin tidak dapat menemukan solusi yang stabil.
Kesalahan dalam Pengimplementasian Algoritma: Terkadang masalah dalam konvergensi bisa berasal dari kesalahan dalam pengkodean algoritma itu sendiri.
Weighted Iterative Jacobi
Adapun solusi yang bisa digunakan adalah dengan menggunakan fungsi
lsolve.jacobi pada package Rlinsolve
berdasarkan buku yang dibuat oleh Demmel JW (1997). Pada fungsi
tersebut, metode yang digunakan adalah Weighted Jacobi atau
pembobotan Jacobi sehingga asumsi nonsingular pada matriks \(D\) dapat terpenuhi kembali. Syntax yang
digunakan dalam fungsi ini adalah sebagai berikut:
library(Rlinsolve)## Warning: package 'Rlinsolve' was built under R version 4.3.2hasil2 <- lsolve.jacobi(A,b,rep(0,3), maxiter=10000, weight = 2/3)
cat("Solusi dari (b0, b1, b2) = (", hasil2$x, ") setelah", hasil2$iter, "iterasi")## Solusi dari (b0, b1, b2) = ( -57.89436 4.709369 0.3378165 ) setelah 3052 iterasiDari hasil di atas, dapat terlihat bahwa nilai koefisien beta konvergen setelah 3052 iterasi dengan menggunakan metode weighted iterative jacobi.
Validasi dengan fungsi base di R
solve(A,b)## [,1]
## intercept -57.9876589
## Girth 4.7081605
## Height 0.3392512Berdasarkan nilai di atas, dapat terlihat bahwa nilai hasil metode
iterasi Jacobi telah cenderung mendekati nilai dari hasil sebenarnya
yang didapat oleh fungsi base R. Oleh karena itu, pencarian koefisien
beta bagi permasalahan pada data trees dengan menggunakan
metode iterasi Jacobi telah berhasil dilakukan.