Optimisasi Statistika - Akar Persamaan SPL

Video Pembelajaran - P2

Video Pembelajaran dapat diakses melalui link berikut : https://ipb.link/materiopstat

Pengantar

Dalam praktikum ini, akan dibahas berbagai metode penyelesaian sistem persamaan linear (SPL) menggunakan pendekatan:

  1. Eliminasi Gauss
  2. Eliminasi Gauss-Jordan
  3. Iterasi Jacobi
  4. Iterasi Gauss-Seidel

Setiap metode memiliki karakteristik dan pendekatan tersendiri, mulai dari transformasi matriks hingga metode iterasi untuk mendekati solusi. Studi kasus pada data trees digunakan untuk mengilustrasikan aplikasi praktis metode ini dalam konteks regresi linear.


Teori dan Metode

1. Eliminasi Gauss

Metode ini bertujuan untuk mengubah matriks augmented SPL ke bentuk segitiga atas (row echelon form) melalui operasi baris elementer. Setelah bentuk ini tercapai, solusi diperoleh melalui substitusi mundur (back substitution).

Langkah Utama:

  1. Transformasi matriks ke bentuk segitiga atas dengan menggunakan operasi baris elementer: \[ R_i \to R_i - \frac{a_{i,j}}{a_{j,j}} R_j, \quad i > j \] dimana \(R_i\) adalah baris ke-\(i\) dalam matriks.

  2. Substitusi mundur untuk menemukan solusi: \[ x_i = \frac{b_i - \sum_{j=i+1}^n a_{i,j}x_j}{a_{i,i}}. \]

Contoh SPL: \[ \begin{aligned} 3x + y - z &= 5 \\ 4x + 7y - 3z &= 20 \\ 2x - 2y + 5z &= 10 \end{aligned} \]

Matriks augmented SPL: \[ \begin{pmatrix} 3 & 1 & -1 & | & 5 \\ 4 & 7 & -3 & | & 20 \\ 2 & -2 & 5 & | & 10 \end{pmatrix} \]

Implementasi dalam R:

# Membuat matriks A dan vektor 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)
}

# Menjalankan eliminasi
A_gauss <- gauss_elimination(A, b)
print(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.94118

Terbentuk 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.650602

Kelebihan dan Kekurangan:

  • Kelebihan:

    • Mudah dipahami dan diimplementasikan untuk SPL kecil.

    • Dapat digunakan pada sistem persamaan dengan koefisien real.

  • Kekurangan:

    • Rentan terhadap error pembulatan pada matriks besar.

    • Membutuhkan operasi komputasi yang signifikan untuk matriks berukuran besar.

2. Eliminasi Gauss-Jordan

Metode ini memperluas Eliminasi Gauss dengan mengubah matriks ke bentuk reduced row echelon form (RREF), di mana elemen diagonal bernilai 1 dan elemen lainnya bernilai 0. Dengan bentuk ini, solusi langsung dapat dibaca tanpa substitusi.

Langkah Utama:

  1. Mengubah elemen diagonal menjadi 1: \[ R_i \to \frac{R_i}{a_{i,i}} \]

  2. Mengeliminasi elemen di atas dan di bawah diagonal: \[ R_j \to R_j - a_{j,i}R_i, \quad j \neq i \]

Implementasi dalam 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)
}

# Menjalankan eliminasi
A_gauss_jordan <- gauss_jordan_elimination(A, b)
print(A_gauss_jordan)
##                   b
## [1,] 1 0 0 1.506024
## [2,] 0 1 0 3.132530
## [3,] 0 0 1 2.650602

Sehingga, bentuk solusi dari x,y,z adalah

A_gauss_jordan[,4]
## [1] 1.506024 3.132530 2.650602

Kelebihan dan Kekurangan:

  • Kelebihan:

    • Memberikan solusi langsung tanpa substitusi mundur.

    • Berguna untuk analisis data atau aljabar linier simbolik.

  • Kekurangan:

    • Lebih kompleks dibandingkan Eliminasi Gauss.

    • Membutuhkan lebih banyak langkah komputasi.

3. Iterasi Jacobi

Metode Jacobi merupakan teknik iteratif yang menggunakan pendekatan berulang untuk menghitung solusi. Matriks A dipecah menjadi:

\[ A = D + R \]

dengan \(D\) adalah matriks diagonal, dan \(R\) adalah matriks residu. Persamaan iteratifnya:

\[ x_i^{(k+1)} = \frac{1}{a_{ii}} \left(b_i - \sum_{j \neq i} a_{ij} x_j^{(k)}\right). \]

Kriteria Konvergensi: Metode Jacobi akan konvergen jika matriks \(A\) memiliki dominansi diagonal: \[ |a_{ii}| > \sum_{j \neq i} |a_{ij}| \quad \text{untuk setiap } i. \]

Implementasi dalam R:

jacobi <- function(A, b, tol = 1e-6, maxiter = 1000) {
  n <- length(b)
  x <- rep(0, n) # Inisialisasi solusi
  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 iterasi Jacobi
jacobi_result <- jacobi(A, b)
print(jacobi_result)
## $x
## [1] 1.506024 3.132530 2.650602
## 
## $iter
## [1] 22

Kelebihan dan Kekurangan:

  • Kelebihan:

    • Cocok untuk matriks besar yang jarang (sparse).

    • Proses iterasi paralel lebih mudah diimplementasikan.

  • Kekurangan:

    • Tidak selalu konvergen, terutama jika tidak ada dominansi diagonal.

    • Membutuhkan iterasi lebih banyak dibandingkan Gauss-Seidel.

4. Iterasi Gauss-Seidel

Berbeda dari Jacobi, Gauss-Seidel menggunakan nilai terbaru dari iterasi saat itu untuk memperbarui solusi:

\[ 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). \]

Kriteria Konvergensi: Sama seperti Jacobi, metode ini membutuhkan matriks \(A\) dengan dominansi diagonal untuk menjamin konvergensi.

Implementasi dalam 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.

Kelebihan dan Kekurangan:

  • Kelebihan:

    • Konvergensi lebih cepat dibandingkan Jacobi.

    • Efisien untuk matriks besar yang jarang.

  • Kekurangan:

    • Tidak selalu konvergen jika matriks tidak memenuhi kriteria dominansi diagonal.

Studi Kasus: Regresi Linear

Data trees digunakan untuk memperkirakan volume pohon (variabel respons) berdasarkan diameter dan tinggi pohon. Matriks prediktor \(X\) dan vektor respons \(y\) diformulasikan sebagai berikut:

Bentuk Matriks: \[ X = \begin{pmatrix} 1 & x_{11} & x_{12} \\ 1 & x_{21} & x_{22} \\ \vdots & \vdots & \vdots \end{pmatrix}, \quad y = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \end{pmatrix} \]

Persamaan Normal: \[ \beta = (X^T X)^{-1} X^T y \]

Metode Eliminasi dan Iterasi digunakan untuk menghitung \(\beta\) dengan validasi dari fungsi lm di R. Hasil regresi diverifikasi dengan berbagai metode untuk memastikan kesesuaian solusi.

head(trees)
##   Girth Height Volume
## 1   8.3     70   10.3
## 2   8.6     65   10.3
## 3   8.8     63   10.2
## 4  10.5     72   16.4
## 5  10.7     81   18.8
## 6  10.8     83   19.7
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 ...

Visualisasi: Grafik hubungan antara diameter, tinggi, dan volume pohon dapat dibuat untuk mendukung analisis:

library(ggplot2)
ggplot(trees, aes(x = Girth, y = Volume)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Hubungan Diameter dan Volume Pohon",
       x = "Diameter Pohon (Girth)",
       y = "Volume Pohon")

GGally::ggpairs(trees)

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     83

Bentuk Vector b

resp <- trees$Volume
head(resp)
## [1] 10.3 10.3 10.2 16.4 18.8 19.7

Transformasi:

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.3392512

Nilai 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.7857
coef_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.3392512

Dapat 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.3392512
coef_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.3392512

Dapat 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] 1000
coef_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.3392512