• Vektor dan matriks Pada Chapter 2.7 dan Chapter 2.8 telah dijelaskan sekilas bagaimana cara melakukan operasi pada vektor dan matriks. Pada chapter ini, penulis akan menambahkan operasi-operasi lain yang dapat dilakukan pada vektor dan matriks. Dasar-dasar operasi ini selanjutnya akan digunakan sebagai dasar menyusun algoritma penyelesaian sistem persamaan linier.
    1. Operasi Vektor Untuk lebih memahami operasi tersebut, berikut penulis berikan contoh penerapannya pada R: 
    u <- seq(1,5)
v <- seq(6,10)

# penjumlahan
u+v
    ## [1]  7  9 11 13 15
    # penguranga
u-v
    ## [1] -5 -5 -5 -5 -5

    Bagaimana jika kita melakukan operasi dua vektor, dimaana salah satu vektor memiliki penjang yang berbeda?. Untuk memnjawab hal tersebut, perhatikan sintaks berikut: 

    x <- seq(1,2)
u+x
    ## Warning in u + x: longer object length is not a multiple of shorter object
## length
    ## [1] 2 4 4 6 6

    Berdasarkan contoh tersebut, R akan mengeluarkan peringatan yang menunjukkan operasi dilakukan pada vektor dengan panjang berbeda. R akan tetap melakukan perhitungan dengan menjumlahkan kembali vektor u yang belum dijumlahkan dengan vektor x sampai seluruh elemen vektor u dilakukan operasi penjumlahan.

    Berikut adalah contoh bagaimana cara menghitung inner product dan panjang vektor menggunakan R: 

    # inner product
u%*%v
    ##      [,1]
## [1,]  130
    # panjang vektor u
sqrt(sum(u*u))
    ## [1] 7.416198
    1. Operasi matriks berikut disajikan contoh operasi penjumlahan pada matriks: 
    A <- matrix(1:9,3)
B <- matrix(10:18,3)
C <- matrix(1:6,3)
# penjumlahan dengan skalar
A+1
    ##      [,1] [,2] [,3]
## [1,]    2    5    8
## [2,]    3    6    9
## [3,]    4    7   10
    # penjumlahan A+B
A+B
    ##      [,1] [,2] [,3]
## [1,]   11   17   23
## [2,]   13   19   25
## [3,]   15   21   27
    # Perkalian matriks
A%*%B
    ##      [,1] [,2] [,3]
## [1,]  138  174  210
## [2,]  171  216  261
## [3,]  204  258  312
    1. Operasi Baris Elementer Terdapat tiga buah operasi dasar pada baris matriksoperasi baris elementer. Ketiga operasi ini akan menjadi dasar operasi sub-chapter selanjutnya. Ketiga operasi dasar tersebut antara lain:

    A. Row Scalling. Mengalikan baris matriks dengan konstanta bukan nol. B. Row Swaping. Menukar urutan baris pada sebuah matriks (contoh: menukar baris 1 dengan baris 2 dan sebaliknya). c. Row Replacement. Baris matriks diganti dengan hasil penjumlahan atau pengurangan baris matriks tersebut dengan baris matriks lainnya, dimana baris matriks lainnya yang akan dijumlahkan/dikurangkan dengan matriks tersebut telah dilakukan proses row scalling. Luaran yang diperoleh pada umumnya adalah nilai nol pada baris matriks awal atau akhir.

    Ketiga proses tersebut akan terjadi secara berulang, khusunya jika kita hendak mengerjakan sistem persamaan linier menggunakan algoritma eliminasi Gauss. Untuk mempermudah proses tersebut, kita dapat membuat masing-masing fungsi untuk masing-masing operasi tersebut. Algoritma fungsi-fungsi tersebut selanjutnya menjadi dasar penyusunan algoritma fungsi-fungsi eliminasi Gauss dan dekomposisi matriks yang akan dijelaskan pada chapter selanjutnya.

    Fungsi row scalling pada R dapat dituliskan pada sintaks berikut: 

    scale_row <- function(m, row, k){
 m[row, ] <- m[row, ]*k
 return(m)
}
# membuat matriks A
(A <- matrix(1:15, nrow=5))
    ##      [,1] [,2] [,3]
## [1,]    1    6   11
## [2,]    2    7   12
## [3,]    3    8   13
## [4,]    4    9   14
## [5,]    5   10   15
    # lakukan scaling pada row 2 dengan nilai 10
scale_row(m=A, row=2, 10)
    ##      [,1] [,2] [,3]
## [1,]    1    6   11
## [2,]   20   70  120
## [3,]    3    8   13
## [4,]    4    9   14
## [5,]    5   10   15

    Row swapping merupakan proses yang berulang, kita perlu menyimpan terlebih dahulu baris matriks pertama kedalam sebuah objek. Baris matriks pertama selanjutnya diganti dengan baris matriks kedua, sedangkan baris matriks kedua selanjutnya akan diganti dengan baris matriks pertama yang telah terlebih dahulu disimpan dalam sebuah objek. Fungsi row swapping pada R dapat dituliskan pada sintaks berikut: 

    swap_row <- function(m, row1, row2){
  row_tmp <- m[row1, ]
  m[row1, ] <- m[row2, ]
  m[row2, ] <- row_tmp
  return(m)
}
# pertukarkan baris 2 dengan baris 5
swap_row(m=A, row1=2, row2=5)
    ##      [,1] [,2] [,3]
## [1,]    1    6   11
## [2,]    5   10   15
## [3,]    3    8   13
## [4,]    4    9   14
## [5,]    2    7   12

    Pada proses row replacement, proses perhitungan dilakukan dengan melakukan penjumlahan suatu baris matriks dengan baris matriks lainnya dengan terlebih dahulu melakukan row scalling terhadap matriks lainnya. Berikut adalah fungsi replace_row() yang ditulis pada R: 

    replace_row <- function(m, row1, row2, k){
  m[row2, ] <- m[row2, ] + m[row1, ]*k
  return(m)
}

    Berikut adalah contoh penerapan fungsi replace_row(): 

    replace_row(m=A, row1=1, row2=3, k=-3)
    ##      [,1] [,2] [,3]
## [1,]    1    6   11
## [2,]    2    7   12
## [3,]    0  -10  -20
## [4,]    4    9   14
## [5,]    5   10   15
    1. Eliminasi Gauss Pada sub-chapter ini kita akan menggunakan operasi baris elementer yang telah dijelaskan pada Chapter 2.5. Terdapat dua topik yang akan dibahas pada sub-chapter ini, yaitu: row echelon form termasuk reduced row echelon form dan matriks tridiagonal.

    Eliminasi Gauss merupakan sebuah cara untuk mencari penyelesaian sistem persamaan linier. Ide dasar dari eliminasi Gauss adalah melakukan operasi matematika pada baris matriks (lihat Chapter 2.5) dan melanjutkannya sampai hanya tersisa satu variabel saja. Kita dapat melakukan lebih dari satu operasi baris elementer pada proses elmininasi ini (contoh: mengalikan sebuah baris dengan konstanta dan menjumlahkan hasilnya pada baris lain).

    1. Row Echelon Form Sebuah matriks merupakan row echelon form jika matriks tersebut memenuhi beberapa kondisi:
    1. Angka bukan nol pertama dari kiri (leading coefficient) selalu di sebelah kanan angka bukan nol pertama pada baris di atasnya.
    2. Baris yang terdiri dari semua nol ada di bagian bawah matriks.

    Teorema 6.1 (spltheorem) Suatu sistem persamaan linier mempunyai penyelesaian tunggal bila memenuhi syarat-syarat sebagai berikut: 1. ukuran persamaan linier simultan bujursangkar (jumlah persamaan sama dengan jumlah variabel bebas). 2. sistem persamaan linier non-homogen di mana minimal ada satu nilai vektor konstanta B tidak nol atau terdapat bn≠0 . 3. Determinan dari matriks koefisiensistem persamaan linier tidak sama dengan nol.

    Algoritma Row Echelon Form 1. Masukkan matriks A ,dan vektor B beserta ukurannya n 2. Buat augmented matrix [A|B] namakan dengan A Untuk baris ke- i dimana i = 1 s/d n ,perhatikan apakah nilai ai,jsama dengan nol. a) Bila iya, lakukan row swapping antara baris ke-i dan baris ke-i+k≤n, dimana ai+k,j tidak sama dengan nol. Bila tidak ada berarti perhitungan tidak bisa dilanjutkan dan proses dihentikan dengan tanpa penyelesaian, b) Bila tidak, lanjutkan. 4. Untuk baris ke- j, dimana j=i+1s/d n,lakukan operasi baris elementer:a) Hitung c=aj,iai,i, b) untuk kolom k, dimana k=1 s/d n+1, hitung aj,k=aj,k−c.ai,k. 5. Hitung akar, untuk i=n s/d 1 (bergerak dari baris pertama) menggunakan Persamaan (6.16).

    Berdasarkan algoritma tersebut, kita dapat menyusun fungsi pada R untuk menyelesaikan sistem persamaan linier menggunakan matriks row echelon form. Fungsi yang akan dibentuk hanya sampai pada algoritma ke-4. Proses substitusi akan dilakukan secara manual. Berikut adalah sintaks yang digunakan: 

    ref_matrix <- function(a){
  m <- nrow(a)
  n <- ncol(a)
  piv <- 1
  
# cek elemen diagonal apakah bernilai nol
  for(row_curr in 1:m){
    if(piv <= n){
      i <- row_curr
      while(a[i, piv] == 0 && i < m){
        i <- i+1
        if(i > m){
          i <- row_curr
          piv <- piv+1
          if(piv > n)
            return(a)
        }
      }
      
# jika diagonal bernilai nol, lakukan row swapping
    if(i != row_curr)
      a <- swap_row(a, i, row_curr)
    
# proses triangulasi untuk membentuk matriks segitiga atas
    for(j in row_curr:m)
      if(j != row_curr){
        c <- a[j, piv]/a[row_curr, piv]
        a <- replace_row(a, row_curr, j, -c)
      }
    piv <- piv+1
    }
  }
  return(a)
}

    Dengan menggunakan fungsi ref_matrix(), kita dapat membentuk matriks row echelon form pada Contoh 6.1. 

    am <- c(1,1,2,
        1,2,1,
        1,-1,2,
        6,2,10)
(m <- matrix(am, nrow=3))
    ##      [,1] [,2] [,3] [,4]
## [1,]    1    1    1    6
## [2,]    1    2   -1    2
## [3,]    2    1    2   10
    ref_matrix(m)
    ##      [,1] [,2] [,3] [,4]
## [1,]    1    1    1    6
## [2,]    0    1   -2   -4
## [3,]    0    0   -2   -6

    Contoh 6.2 Dengan menggunakan fungsi ref_matrix(), buatlah matriks row echelon form dari sistem persamaan linier berikut: 2x1+x2−x3=1 3x1+2x2−2x3=1 x1−5x2+4x3=3 Penyelesaian matriks tersebut adalah sebagai berikut: 

    (m <- matrix(c(2,3,1,
              1,2,-5,
              -1,-2,4,
              1,1,3), nrow=3))
    ##      [,1] [,2] [,3] [,4]
## [1,]    2    1   -1    1
## [2,]    3    2   -2    1
## [3,]    1   -5    4    3
    ref_matrix(m)
    ##      [,1] [,2] [,3] [,4]
## [1,]    2  1.0 -1.0  1.0
## [2,]    0  0.5 -0.5 -0.5
## [3,]    0  0.0 -1.0 -3.0

    Proses lebih lanjut akan menghasilkan penyelesaian sebagai berikut: x1=1 x2=2 x3=3

    1. Eliminasi Gauss-Jordan Berbeda dengan metode eliminasi Gauss yang telah dijelaskan pada Chapter 6.3.1, metode eliminasi Gauss-Jordan membentuk matriks menjadi bentuk reduced row echelon form. Metode ini merupakan pengembangan metode eliminasi Gauss, dimana matriks sebelah kiri augmented matrix diubah menjadi matriks diagonal algorritma eliminasi gouss-jordan kita dapat membangun sebuah fungsi menggunakan R. Fungsi tersebut adalah sebagai berikut: 
    gauss_jordan <- function (a){
    m <- nrow (a)
    n <- ncol (a)
    piv <- 1
    
# cek elemen diagonal utama apakah bernilai nol
    for(row_curr in 1:m){
        if(piv <= n){
            i <- row_curr
            while(a[i, piv] == 0 && i < m){
                i <- i + 1
                if(i > m){
                    i <- row_curr
                    piv <- piv + 1
                    if(piv > n)
                        return (a)
                }
            }

# jika diagonal utama bernilai nol,lakukan row swapping
            if(i != row_curr)
                a <- swap_row(a, i, row_curr)
            
# proses pembentukan matriks reduced row echelon form
            piv_val <- a[row_curr , piv]
            a <- scale_row (a, row_curr , 1 / piv_val)
            for(j in 1: m){
                if(j != row_curr){
                    k <- a[j, piv]/a[row_curr, piv]
                    a <- replace_row (a, row_curr, j, -k)
                }
            }
            piv <- piv + 1
        }
    }
    return (a)
}

    Dengan menggunakan fungsi gauss_jordan(), sistem persamaan linier : 

    (m <- matrix(c(1,2,1,4,3,8), nrow=2))
    ##      [,1] [,2] [,3]
## [1,]    1    1    3
## [2,]    2    4    8
    gauss_jordan(m)
    ##      [,1] [,2] [,3]
## [1,]    1    0    2
## [2,]    0    1    1

    Contoh 6.4 Dengan menggunakan fungsi gauss_jordan(), carilah penyelesaian sistem persamaan linier pada Contoh 6.1 dan Contoh 6.2:

    Jawab: Untuk Contoh 6.1: 

    am <- c(1,1,2,
        1,2,1,
        1,-1,2,
        6,2,10)
m <- matrix(am, nrow=3)

gauss_jordan(m)
    ##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0    1
## [2,]    0    1    0    2
## [3,]    0    0    1    3

    Untuk Contoh 6.2: 

    m <- matrix(c(2,3,1,1,2,-5,
              -1,-2,4,1,1,3), 
            nrow=3)
gauss_jordan(m)
    ##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0    1
## [2,]    0    1    0    2
## [3,]    0    0    1    3
    1. Matrik Tridiagonal Metode eliminasi Gauss merupakan metode yang sederhana untuk digunakan khususnya jika semua koefisien bukan nol berkumpul pada diagonal utama dan beberapa diagonal sekitarnya. Suatu sistem yang bersifat demikian disebut sebagai banded dan banyaknya diagonal yang memuat koefisien bukan nol disebut sebagai bandwidth. Contoh khusus yang sering dijumpai adalah matriks tridiagonal yang memiliki bandwidth tiga. Proses eliminasi untuk matriks tridiagonal bersifat trivial karena dengan membentuk sebuah subdiagonal tambahan, proses substitusi mundur segera dapat dilakukan.

    fungsi pada R yang dapat dibangun untuk menyelesaikan permasalahan dari algoritma matriks tradiagonal sebagai berikut: 

    tridiagmatrix <- function (L, D, U, b){
  n <- length (D)
  L <- c(NA , L)
  
  ## forward sweep
  U[1] <- U[1] / D[1]
  b[1] <- b[1] / D[1]
  for(i in 2:(n - 1)){
      U[i] <- U[i] / (D[i] - L[i] * U[i - 1])
      b[i] <- (b[i] - L[i] * b[i - 1]) /
      (D[i] - L[i] * U[i - 1])
  }
  b[n] <- (b[n] - L[n] * b[n - 1])/(D[n] - L[n] * U[n - 1])
  
  ## backward sweep
  x <- rep.int (0, n)
  x[n] <- b[n]
  for(i in (n - 1) :1)
      x[i] <- b[i] - U[i] * x[i + 1]
  return (x)
}

    Contoh 6.5 Selesaikan sistem persamaan berikut menggunakan fungsi tridiagmatrix() dan fungsi gauss_jordan()! 3x1+4x2=20 4x1+5x2−2x3=28 2x2+5x3−3x4=18 3x3+5x4=18 Untuk menyelesaikan persamaan tersebut menggunakan fungsi tridiagmatrix(), kita perlu membentuk vektor diagonal
    l, d, u, dan b. 

    l <- u <- c(4, 2, 3); d <- c(3, 5, 5, 5)
b <- c(20, 28, 18, 18)
#Setelah terbentuk, vektor tersebut dapat langsung dimasukkan ke dalam fungsi tridiagmatrix().
tridiagmatrix(L=l, D=d, U=u, b=b)
    ## [1] 4 2 1 3

    Untuk menyelesaikannya menggunakan fungsi gauss_jordan(), kita perlu membentuk augmented matrix-nya terlebih dahulu. 

    m <- matrix(c(3,4,0,0,4,5,2,0,
              0,2,5,3,0,0,3,5,
              20,28,18,18), nrow=4)
gauss_jordan(m)
    ##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    0    0    0    4
## [2,]    0    1    0    0    2
## [3,]    0    0    1    0    1
## [4,]    0    0    0    1    3
    1. Penyelesaian Sistem Persamaan Linier Menggunakan Fungsi solve() R menyediakan fungsi bawaan solve() untuk menyelesaiakan sistem persamaan linier. Format fungsi solve() adalah sebagai berikut: 
    #solve(a,b)

#Catatan:
#a: matriks koefisien atau matriks segiempat
#b: vektor konstanta

    Berikut adalah contoh penerapan fungsi solve() pada sistem persamaan linier yang disajikan pada Contoh 6.2: 

    # memecah matriks m menjadi matriks koefisien dan vektor konstanta
a <- matrix(c(2,3,1,1,2,-5,-1,-2,4),nrow=3)
b <- c(1,1,3)

solve(a,b)
    ## [1] 1 2 3

    Jika kita hanya memasukkan matriks persegi, maka output yang akan dihasilkan adalah invers dari matriks yang kita masukkan. 

    solve(a)
    ##      [,1] [,2]          [,3]
## [1,]    2   -1  7.401487e-17
## [2,]   14   -9 -1.000000e+00
## [3,]   17  -11 -1.000000e+00

    Jika kita mengalikan invers dengan matriks semula, maka akan dihasilkan output berupa matriks identitas. 

    a%*%solve(a)
    ##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1
    1. Penyelesaian Sistem Persamaan Linier Menggunakan Fungsi ’Solve.tridiag()` Penyelesaian matriks tridiagonal selain menggunakan fungsi solve(), juga dapat menggunakan fungsi Solve.tridiag() dari Paket limSolve. Untuk menginstall dan mengaktifkan Paket tersebut, jalankan sintaks berikut:

    Fungsi Solve.tridiag() memiliki format sebagai berikut: 

    library(limSolve)
    ## Warning: package 'limSolve' was built under R version 4.2.2
    #Solve.tridiag ( diam1, dia, diap1, B=rep(0,times=length(dia)))


#Catatan:
#diam1: vektor bukan nol di bawah diagonal matriks
#dia: vektor bukan nol pada diagonal matriks
#diap1: vektor bukan nol di atas diagonal matriks
#B: vektor konstanta

    Untuk memahami penerapannya, kita akan menggunakan kembali matriks yang ada pada Contoh 6.5. 

    l <- u <- c(4, 2, 3); d <- c(3, 5, 5, 5)
b <- c(20, 28, 18, 18)
Solve.tridiag(diam1=l, dia=d, diap1=u, B=b)
    ##      [,1]
## [1,]    4
## [2,]    2
## [3,]    1
## [4,]    3

    1. Dekomposisi Matriks Seringkali kita diminta untuk memperoleh nilai penyelesaian suatu persamaan linier Ax=B, dimana nilai vektor B yang selalu berubah-ubah. Penggunaan metode eliminasi Gauss mengharuskan untuk menyelesaikan sistem persamaan linier Ax=B secara terpisah untuk setiap perubahan vektor B. Untuk menghindari pekerjaan eliminasi yang selalu berulang-ulang, faktorisasi menjadi suatu hal yang dapat dilakukan untuk mempersingkat prosesnya. Faktorisasi atau dekomposisi matriks merupakan suatu algoritma untuk memecah matriks A, hasil pemecahan ini selanjutnya digunakan untuk memperoleh penyelesaian sistem persamaan linier melalui perkalian antara vektor B dan hasil faktorisasi matriks A.

    2. Dekomposisi LU Misalkan kita memiliki persamaan linier seperti yang ditunjukkan oleh Persamaan (6.12). Pada metode dekomposisi LU, matriks A difaktorkan menjadi matriks L dan matriks U, dimana ukuran kedua matriks tersebut harus sama dengan ukuran matriks A atau dapat kita tuliskan bahwa hasil perkalian kedua matriks tersebut akan menghasilkan matriks A.

    A=LU (6.22)

    Sehingga Persamaan (6.12) akan menjadi Persamaan (6.23).

    LUx=b 6.23)

    Langkah penyelesaian sistem persamaan linier, diawali dengan menghadirkan vektor t yang ditunjukkan pada Persamaan (6.24).

    Ux=t (6.24)

    Langkah pada Persamaan (6.24) tidak dimaksudkan untuk menghitung vektor t, melainkan untuk menghitung vektor x. Vektor t diperoleh dengan menggunakan Persamaan (6.25).

    Lx=t 6.25)

    Kita dapat menyelesaikan sistem persamaan yang ditunjukkan pada Persamaan (6.24) dan Persamaan (6.25) menggunakan berbagai algoritma penyelesaian yang telah dibahas sebelumnya. Namun, karena matriks L merupakan matriks segitiga bawah dengan nilai nol berada pada bagian atas diagonal utama, penyelesaian t mengambil langkah yang lebih sedikit. Kondisi ini sama dengan kondisi penyelesaian matriks tridiagonal, dimana kita memanfaatkan sejumlah jalan pintas penyelesaiaannya guna mempercepat komputasi. Matriks segitia bawah L akan berupa matriks persegi dengan ukuran m, di mana m merupakan jumlah baris matriks A. Persamaan (6.25) dalam bentuk matriks akan terlihat seperti Persamaan (6.26).

    Berdasarkan algoritma Dekomposisi LU, kita dapat menyusun algoritma faktorisasi LU menggunakan R. Berikut adalah sintaks yang digunakan: 

    lu_solve <- function(a, b=NULL){
    m <- nrow(a)
    n <- ncol(a)
    piv <- 1

# membentuk matriks identitas P dan L
    P <- L <- diag(n)

# cek elemen diagonal utama apakah bernilai nol
    for(row_curr in 1:m){
        if(piv <= n){
            i <- row_curr
            while(a[i, piv] == 0 && i < m){
                i <- i + 1
                if(i > m){
                    i <- row_curr
                    piv <- piv + 1
                    if(piv > n)
                        return(list(P = P, L = L, U = a))
                }
            }
            
# jika elemen diagonal utama bernilai nol,lakukan row swapping
            if(i != row_curr){
                a <- swap_row(a, i, row_curr)
                P <- swap_row(P, i, row_curr)
            }
            
  # pembentukan matriks L dan U
            for(j in row_curr:m)
                if(j != row_curr){
                    k <- a[j, piv]/a[row_curr, piv]
                    
  # matriks U
                    a <- replace_row(a, row_curr, j, -k)
                    
  # pengisian elemen matriks L
                    L[j, piv] <- k
                }
            piv <- piv + 1
        }
    }
    
# penyelesaian persamaan linier
    if(is.null(b)){
      return(list(P = P, L = L, U = a))
    }else{
      
      # forward substitution
      t <- forwardsolve(L, b)
      
      # backward substitution
      x <- backsolve(a, t)
      return(list(P = P, L = L, U = a, result=x))
     }
}

    Kita dapat menyelesaikan sistem persamaan linier pada Contoh 6.6 menggunakan fungsi yang telah kita buat. 

    # membuat matriks a dan vektor b
a <- matrix(c(1,2,3,-1,1,1,-1,2,
              0,-1,-1,3,3,1,2,-1),
            nrow=4)
b <- c(4,1,-3,4)

# penyelesaian
decomp<-lu_solve(a,b)

    Untuk membentuk kembali matriks A, kita dapat mengalikan matriks L,U, dan P. 

    decomp$L%*%decomp$U%*%decomp$P
    ##      [,1] [,2] [,3] [,4]
## [1,]    1    1    0    3
## [2,]    2    1   -1    1
## [3,]    3   -1   -1    2
## [4,]   -1    2    3   -1

    Contoh 6.7 Lakukan dekomposisi LU pada matriks berikut dan lakukan pengecekan apakah perkalian hasil dekomposisi matriks akan menghasilkan matriks semula!

    ⎡0 1 -1⎤ ⎢1 5 9⎥ ⎣7 -1 -5⎦ Jawab: Lakukan proses dekomposisi menggunakan fungsi lu_solve(). 

    library(limSolve)
# membentuk matriks a
(A <- matrix(c(0, 1, 7, 1, 5, -1, -2, 9, -5), 3))
    ##      [,1] [,2] [,3]
## [1,]    0    1   -2
## [2,]    1    5    9
## [3,]    7   -1   -5
    # dekomposisi lu
#decomp <-lu_solve(A)

    Lakukan pengecekan apakah matriks hasil dekomposisi akan menghasilkan matriks A. 

    decomp$P %*% decomp$L %*% decomp$U
    ##      [,1] [,2] [,3] [,4]
## [1,]    1    1    0    3
## [2,]    2    1   -1    1
## [3,]    3   -1   -1    2
## [4,]   -1    2    3   -1

    Fungsi lu() pada Paket Matrix dapat digunakan untuk melakukan dekomposisi LU. Untuk meggunakan fungsi tersebut, kita harus menginstall dan mengaktifkan Paket Matrix. 

    library(Matrix)
    ## Warning: package 'Matrix' was built under R version 4.2.2
    #Untuk dapat menggunakannya kita hanya perlu menginputkan matriks kedalam fungsi tersebut. Berikut adalah contoh penerapannya:

# membuat matriks a 
a <- Matrix::Matrix(round(rnorm(9),2), nrow=3)

# dekomposisi
lum <- Matrix::lu(a)
lum
    ## 'MatrixFactorization' of Formal class 'denseLU' [package "Matrix"] with 4 slots
##   ..@ Dimnames:List of 2
##   .. ..$ : NULL
##   .. ..$ : NULL
##   ..@ x       : num [1:9] -1.69 -0.604 0.112 -0.66 1.752 ...
##   ..@ perm    : int [1:3] 2 3 3
##   ..@ Dim     : int [1:2] 3 3

    Untuk menampilkan hasil dekomposisi, jalankan fungsi expand(). 

    decomp <- Matrix::expand(lum)
decomp
    ## $L
## 3 x 3 Matrix of class "dtrMatrix" (unitriangular)
##      [,1]       [,2]       [,3]      
## [1,]  1.0000000          .          .
## [2,] -0.6035503  1.0000000          .
## [3,]  0.1124260 -0.1859947  1.0000000
## 
## $U
## 3 x 3 Matrix of class "dtrMatrix"
##      [,1]       [,2]       [,3]      
## [1,] -1.6900000 -0.6600000 -0.4500000
## [2,]          .  1.7516568 -1.1715976
## [3,]          .          . -0.3173192
## 
## $P
## 3 x 3 sparse Matrix of class "pMatrix"
##           
## [1,] . . |
## [2,] | . .
## [3,] . | .
    1. Dekomposisi Cholesky Dekomposisi Cholesky memberikan faktorisasi matriks alternatif sehingga A=LLTA=LLT, di mana LTLT merupakan transpose konjugat dari matriks LL. Dalam kasus ini, penulis hanya bekerja dengan matriks rill dengan nilai rill dan bagian imajiner nol. Jadi untuk tujuan sub-chapter ini, matriks LTLT hanyalah transpose dari matriks LL. Seperti dekomposisi LU, dekomposisi Cholesky dapat digunakan untuk menyelesaikan sistem persamaan linier. Kelebihannya, Menemukan dekomposisi Cholesky jauh lebih cepat daripada dekomposisi LU. Namun, dekomposisi ini hanya terbatas pada matriks tertentu saja. Dekomposisi Cholesky hanya dapat digunakan pada matriks definit positif dan simetris. Matriks simeteris merupakan matriks yang nilai di atas dan di bawah diagonalnya simetris atau sama; secara matematis, untuk semua ii dan jj pada matriks AA, ai;j=aj;iai;j=aj;i. Definit positif berarti bahwa setiap entri pivot (nilai elemen diagonal utama) selelu bernilai positif. Selain itu, untuk matriks definit positif, hubungan xAx>0xAx>0 untuk semua vektor, xx. Karena L∗L∗ transpose dari matriks LL, maka lTi,j=lj,ili,jT=lj,i untuk semua nilai ii dan jj. Tanpa kendala (constraint) ini, dekomposisi Cholesky akan mirip dekomposisi LU. Tetapi dengan kendala ini, nilai elemen matriks LL dan LTLT harus dipilih dengan cermat sehingga hubungan A=LLTA=LLT berlaku.

    Berdasarkan algoritma Dekomposisi Cholesky, kita dapat menyusun fungsi pada R untuk melakukan dekomposisi Cholesky. Fungsi tersebut disajikan pada sintaks berikut: 

    cholesky_solve <- function(a, b=NULL){
    m <- nrow(a)
   
# membentuk matriks L dengan elemen nol
    L = diag(0,m)

# Perhitungan elemen matriks L
    for(i in 1:m){
        for(k in 1:i){
            p_sum <- 0
            for(j in 1:k)
                p_sum <- p_sum + L[j,i]*L[j,k]
                
# Pehitungan elemen diagonal utama
            if(i==k)
                L[k,i]<-sqrt(a[i,i]-p_sum)
            else
                L[k,i]<-(a[k,i]-p_sum)/L[k,k]
        }
    }
    
# Perhitungan elemn matriks L*
    tL <- t(L)
    
# penyelesaian persamaan linier
    if(is.null(b)){
      return(list(L = L, tL = tL, a = a))
    }else{
      
      # forward substitution
      t <- forwardsolve(L, b)
      
      # backward substitution
      x <- backsolve(tL, t)
      return(list(L = L, tL = tL, a = a, result=x))
     }
}

    Contoh 6.8 Dengan menggunakan fungsi cholesky_solve(), lakukan dekomposisi pada matriks berikut! Lakukan pengecekan pada hasil dekomposisi apakah hasil kali matriks dekomposisi akan menghasilkan matriks semula!

    ⎡9 -3 6⎤ ⎢-3 17-10⎥ ⎣6 -10 12⎦

    Jawab: Dekomposisi Cholesky menggunakan fungsi cholesky_solve(), disajikan pada sintaks berikut: 

    library(MatrixModels)
    ## Warning: package 'MatrixModels' was built under R version 4.2.2
    a <- matrix(c(9,-3,6,-3,17,-10,6,-10,12),3)

# dekomposisi Cholesky
#(decomp<-cholesky_solve(a))
# mengecek hasil dekomposisi
#decomp$tL %*% decomp$L

    Fungsi lain yang dapat digunakan untuk melakukan dekomposisi Cholesky adalah menggunakan fungsi chol() pada Paket Matrix. Pada fungsi tersebut, kita hanya perlu menginputkan objek matrik kedalamnya. Berikut adalah contoh penerapan fungsi tersebut menggunakan matriks pada Contoh 6.8. 

    chol(a)
    ##      [,1] [,2] [,3]
## [1,]    3   -1    2
## [2,]    0    4   -2
## [3,]    0    0    2
    #Penting!!!
#Fungsi chol() hanya menampilkan matriks LT. Untuk menampilkan matriks L, kita perlu melakukan transpose
    1. Dekomposisi Lainnya Terdapat beberapa algoritma lain yang telah dikembangkan untuk melakukan dekomposisi matriks. Pada buku ini hanya akan dijelaskan secara singkat terkait fungsi yang digunakan dalam melakukan dekomposisi matriks. Algoritma yang akan dijelaskan pada sub-chapter ini antara lain: QR, singular value decomposition (SVD), dan dekomposisi eigen. Untuk algoritma lainnya, pembaca dapat membaca buku terkait atau mengecek dokumentasinya pada Paket base.
    1. Dekomposisi QR Dekomposisi QR merupakan dekomposisi yang penting dalam menyelesaikan sistem persamaan linier. Dekomposisi ini juga berperan penting untuk menghitung koefisien regresi dan pengaplikasian algoritma Newton-Raphson. Untuk memperoleh informasi terkait dekomposisi ini, pembaca dapat mengetikkan sintaks berikut pada R: 
    ?qr
    ## starting httpd help server ... done
    #Berikut merupakan contoh penerapan fungsi qr() untuk menyelesaikan sistem persamaan linier:

# membuat matriks A dan B
set.seed(123)
A <- matrix((1:12)+rnorm(12), nrow=4)
b <- 2:5

# dekomposisi matriks A
qr(A)
    ## $qr
##            [,1]        [,2]       [,3]
## [1,] -6.3777985 -12.1257372 -19.850120
## [2,]  0.2774974  -6.3105459  -7.939245
## [3,]  0.7147777  -0.6461294   2.351193
## [4,]  0.6382310  -0.5653624   0.276672
## 
## $rank
## [1] 3
## 
## $qraux
## [1] 1.068915 1.512720 1.960964
## 
## $pivot
## [1] 1 2 3
## 
## attr(,"class")
## [1] "qr"
    # memperoleh penyelesaian SPL
qr.solve(A,b)
    ## [1]  0.3045952 -0.1111081  0.3236862
    1. Singular Value Decomposition Singular value decomposition (SVD) merupakan algoritma faktorisasi matriks yang mendekomposisi matriks segiempat menjadi matriks UDVHUDVH, dimana DD merupakan mmatriks diagonal non negatif, UU dan VV merupakan matriks unitary, dan VHVH merupakan matriks tanspose konjugat dari matriks VV. Algoritma ini banyak digunakan dalam analisis principal component. Pada R, SVD dapat dilakukan menggunakan fungsi svd() dari Paket base. Berikut adalah sintaks untuk memperoleh informasi terkait fungsi tersebut: 
    ?svd

#Berikut adalah contoh penerapan fungsi svd():

# dekomposisi matriks A
svd(A)
    ## $d
## [1] 26.094305  2.727495  1.329585
## 
## $u
##            [,1]       [,2]       [,3]
## [1,] -0.3684647 -0.5661362  0.6650563
## [2,] -0.4706958 -0.5703414 -0.6113237
## [3,] -0.5740273  0.4324050 -0.2589679
## [4,] -0.5596176  0.4089332  0.3419343
## 
## $v
##            [,1]        [,2]       [,3]
## [1,] -0.2257101  0.87169376 -0.4349769
## [2,] -0.5201583 -0.48536069 -0.7027520
## [3,] -0.8237052  0.06763865  0.5629696
    1. Dekomposisi Eigen Proses umum yang digunakan untuk menemukan nilai eigen dan vektor eigen suatu matriks segiempat dapat dilihat sebagai proses dari dekomposisi eigen. Proses ini akan mendekomposisi matriks menjadi VDV−1VDV−1, dimana DD merupakan matriks diagonal yang terbentuk dari nilai eigen, dan VV merupakan vektor eigen. Proses dekomposisi ini akan berguna bagi pembaca yang ingin mempelajari principal component analysis. Fungsi eigen() pada Paket base dapat digunakan untuk melakukan dekomposisi eigen. Untuk mempelajari lebih jauh terkait fungsi ini, pambaca dapat menjalankan sintaks berikut: 
    ?eigen

#Berikut adalah contoh sintaks untuk melakukan dekomposisi eigen:

A <- matrix(c(2,-1,0,-1,2,-1,0,-1,2), nrow=3)

# dekomposisi matriks A
eigen(A)
    ## eigen() decomposition
## $values
## [1] 3.4142136 2.0000000 0.5857864
## 
## $vectors
##            [,1]          [,2]      [,3]
## [1,] -0.5000000 -7.071068e-01 0.5000000
## [2,]  0.7071068  1.099065e-15 0.7071068
## [3,] -0.5000000  7.071068e-01 0.5000000
    1. Metode Iterasi Pada Chapter 6.5 kita akan membahas penyelesaian persamaan linier dengan menggunakan metode iterasi. Terdapat dua metode iterasi yang akan dibahas yaitu iterasi Jacobi dan Gauss-Seidel.

    Metode iterasi dimulai dengan estimasi nilai akhir. Setelah menerapkan beberapa perlakuan pada nilai estimasi, hasil perlakuan selanjutnya menjadi nilai estimasi untuk iterasi berikutnya. Proses tersebut akan berlangsung secara terus-menerus hingga ambang batas dipenuhi. Nilai ambang batas dapat berupa jumlah iterasi maksimum atau selisih antara nilai estimasi baru dan estimasi semula lebih kecil dari suatu nilai toleransi yang ditetapkan.

    Jumlah kuadrat merupakan metode yang sering digunakan untuk mengecek apakah selisih nilai estimasi baru terhadap estimasi lama lebih kecil dari nilai toleransi yang ditetapkan.

    1. Iterasi Jacobi Untuk menyelesaikan matriks menggunakan metode iterasi, kita dapat mulai dengan premis terdapat matriks AA dan vektor xx dan b, sehingga Ax=bAx=b. Dengan menggunakan metode Jacobi, pertama-tama kita dapat amati bahwa terdapat matriks RR dan DD yang memiliki hubungan A=R+DA=R+D. Berdasarkan kedua hubungan tersebut, dapat diturunkan operasi matriks melalui persamaan berikut: Ax=b (6.34)
      Rx+Dx=b (6.35) Dx=b−Rx (6.36) x=D−1(b−Rx) (6.37) Persamaan (6.37) merupakan persamaan yang dapat kita gunakan untuk memperoleh nilai x. Jika kita menulis kembali persamaan tersebut, maka kita akan memperoleh persamaan yang digunakan sebagai acuan iterasi Jacobi.

    xn+1=D−1(b−Rxn) (6.38)

    dimana DD merupakan matriks diagonal dengan nilai elemen diagonal berupa diagonal utama matriks AA. Invers dari matriks DD secara sederhana sebagai matriks diagonal sama dengan satu dibagi dengan elemen diagonal utama matriks AA. Matriks RR identik dengan matriks AA. Namun, diagonal utamanya bernilai nol. Suatu iterasi dikatakan konvergen jika jumlah kuadrat dari vektor x(n+1)x(n+1) dan vektor x(n)x(n) semakin mengecil. Suatu persamaan linier yang hendak diselesaikan dengan menggunakan metode iterasi Jacobi harus memenuhi syarat nilai elemen diagonal utama matriks harus lebih dominan. Maksudnya adalah nilai absolut diagonal utama matriks harus lebih besar dari jumlah nilai absolut elemen matriks lainnya pada satu kolom.

    Contoh 6.9 Selesaikan sistem persamaan berikut menggunakan iterasi Jacobi!

    ⎡5 2 3⎤ 40 ⎢2 7 4⎥x= 39
    ⎣1 3 8⎦ 55

    Jawab: Berdasarkan matriks AA (matriks koefisien), kita dapat memastikan bahwa matriks tersebut memiliki nilai dominan pada elemen diagonal utama. Sebagai contoh: |5|>|2|+|1| (kolom 1)|5|>|2|+|1| (kolom 1) |7|>|2|+|3| (kolom 2)|7|>|2|+|3| (kolom 2) Untuk mempermudah proses iterasi, kita akan menggunakan bantuan R untuk melakukan komputasi. Langkah pertama yang perlu dilakukan adalah menyiapkan matriks AA, vektor bb, dan vektor xx (nilai taksiran awal). 

    (A <- matrix(c(5,2,1,2,7,3,3,4,8), 3))
    ##      [,1] [,2] [,3]
## [1,]    5    2    3
## [2,]    2    7    4
## [3,]    1    3    8
    (b <- c(40,39,55))
    ## [1] 40 39 55
    (x <- rep(0,3))
    ## [1] 0 0 0
    #Langkah selanjutnya adalah memperoleh invers matriks D

(Dinv <- diag(1/diag(A)))
    ##      [,1]      [,2]  [,3]
## [1,]  0.2 0.0000000 0.000
## [2,]  0.0 0.1428571 0.000
## [3,]  0.0 0.0000000 0.125
    #Persiapan terakhir sebelum iterasi dilakukan adalah menyiapkan matriks R

(R<-A-diag(diag(A)))
    ##      [,1] [,2] [,3]
## [1,]    0    2    3
## [2,]    2    0    4
## [3,]    1    3    0

    Iterasi selanjutnya dilakukan menggunakan Persamaan (6.38). 

    #literasi 1
(x1 <- Dinv %*% (b-R%*%x))
    ##          [,1]
## [1,] 8.000000
## [2,] 5.571429
## [3,] 6.875000
    #iterasi 2
(x2 <- Dinv %*% (b-R%*%x1))
    ##            [,1]
## [1,]  1.6464286
## [2,] -0.6428571
## [3,]  3.7857143
    #iterasi 3
(x3 <- Dinv %*% (b-R%*%x2))
    ##          [,1]
## [1,] 5.985714
## [2,] 2.937755
## [3,] 6.910268

    Selama proses iterasi,jumlah akar jumlah kuadrat dihitung. Sebagai contoh berikut disajikan akar jumlah kuadrat pada iterasi ke-3: 

    sqrt(sum(x3-x2)^2)
    ## [1] 11.04445

    Berdasarkan algoritma Iterasi Jacobi, kita dapat menyusun fungsi sebuah fungsi untuk melakukan iterasi Jacobi. Berikut sintaks yang digunakan: 

    jacobi <- function(a, b, tol=1e-7, maxiter=100){
  n <- length(b)
  iter <- 0
  
  Dinv <- diag(1/diag(a))
  R <- a-diag(diag(a))
  x <- rep(0,n)
  x_new <- rep(tol, n)
  
  while(sqrt(sum(x_new-x)^2)>tol){
            if(iter>maxiter){
              warning("iterasi maksimum tercapai")
              break
            }
            x <- x_new
            x_new <- Dinv %*% (b - R %*% x)
            iter <- iter+1
  }
  return(list(X = x_new, iter=iter))
  
}

    Berikut adalah penerpan fungsi jacobi() tersebut: 

    jacobi(A,b)
    ## $X
##      [,1]
## [1,]    4
## [2,]    1
## [3,]    6
## 
## $iter
## [1] 62

    Contoh 6.10 Selesaikan sistem persamaan berikut menggunakan fungsi jacobi()

    ⎡27 6 -1⎤ 85 ⎢6 15 2⎥x= 72
    ⎣1 1 54⎦ 110

    Jawab: Matriks A (matriks koefisien) berdasarkan sistem persamaan linier tersebut telah memenuhi syarat dari algoritma Jacobi (nilai diagonal utama dominan dibanding nilai lainnya pada satu kolom). Penyelesaian sistem persamaan tersebut, sebagai berikut: 

    A <- matrix(c(27,6,1,6,15,1,-1,2,54), 3)
b <- c(85,72,110)

jacobi(A,b)
    ## $X
##          [,1]
## [1,] 2.425476
## [2,] 3.573016
## [3,] 1.925954
## 
## $iter
## [1] 17
    #Nilai vektor x sesungguhnya dapat diperoleh menggunakan fungsi solve().

solve(A,b)
    ## [1] 2.425476 3.573016 1.925954
    1. Iterasi Gauss-Seidel Metode iterasi Gauss-Seidel melakukan dekomposisi pada matriks AA menjadi matriks segitiga atas UU dan matriks segitiga bawah LL. Dekomposisi ini tidak sama dengan dekomposisi LU pada Chapter 6.4.1. Matriks UU pada metode Gauss-Seidel merupakan elemen (entri) matriks AA pada bagian atas diagonal utama, sedangkan matriks LL merupakan elemen diagonal utama dan bagian bawah diagonal utama matriks AA. Elemen selain yang penulis sebutkan pada kedua matriks tersebut akan bernilai nol. Persamaan iterasi Gauss-Seidel ditampilkan pada Persamaan (6.39).

    xn+1=L−1(b−Uxn)(6.39)

    Syarat agar suatu sistem persamaan linier dapat diselesaikan menggunakan metode Gauss-Seidel adalah matriks harus memiliki nilai diagonal utama yang dominan. Maksudnya, nilai absolut diagonal utama lebih besar dari jumlah nilai absolut elemen lainnya dalam satu kolom. Jika syarat ini tidak terpenuhi maka metode ini tidak akan memperoleh penyelesaian yang konvergen.

    Contoh 6.11 Selesaikan sistem persamaan pada Contoh 6.10 menggunakan iterasi Gauss-Seidel! Jawab: Kita akan kembali menggunakan bantuan R untuk melakukan kalkulasi pada proses iterasi Gauss-Seidel. Kita telah melakukan pengecekan pada sistem persamaan linier pada contoh tersebut dan menghasilkan kesimpulan bahwa persamaan linier tersebut dapat diselesaikan dengan metode Gauss-Seidel. Langkah selanjutnya adalah membentuk matriks LL dan UU. 

    # membentuk matriks U dan L dari matriks A
(L <- U <- A)
    ##      [,1] [,2] [,3]
## [1,]   27    6   -1
## [2,]    6   15    2
## [3,]    1    1   54
    # membentuk matriks L dari entri bagian bawah diagonal utama matriks A
L[upper.tri(A, diag=FALSE)]<-0
L
    ##      [,1] [,2] [,3]
## [1,]   27    0    0
## [2,]    6   15    0
## [3,]    1    1   54
    # membentuk matriks L dari entri bagian bawah diagonal utama matriks A
L[upper.tri(A, diag=FALSE)]<-0
L
    ##      [,1] [,2] [,3]
## [1,]   27    0    0
## [2,]    6   15    0
## [3,]    1    1   54

    Selanjutya lakukan invers terhadap matriks L menggunakn fungsi solve(). 

    (Linv <- solve(L))
    ##               [,1]         [,2]       [,3]
## [1,]  0.0370370370  0.000000000 0.00000000
## [2,] -0.0148148148  0.066666667 0.00000000
## [3,] -0.0004115226 -0.001234568 0.01851852

    tetapkan nilai estimasi awal dan nilai toleransi yang dikehendaki. Nilai toleransi pada proses ini ditetapkan sebesar 10−7. 

    # tebakan awal nilai x
(x <- rep(0, length(b)))
    ## [1] 0 0 0
    #Lakukan iterasi menggunakan Persamaan (6.39).

#Iterasi 1
(x1 <- Linv %*% (b - U %*% x))
    ##          [,1]
## [1,] 3.148148
## [2,] 3.540741
## [3,] 1.913169
      # akar jumlah kuadrat
sqrt(sum(x1-x)^2)
    ## [1] 8.602058
    #Iterasi 2
(x2 <- Linv %*% (b - U %*% x1))
    ##             [,1]
## [1,] -0.71597317
## [2,]  0.03130011
## [3,]  0.01267913
      # akar jumlah kuadrat
sqrt(sum(x2-x1)^2)
    ## [1] 9.274052

    Iterasi terus dilakukan sampai dengan nilai akar jumlah kuadrat lebih kecil dari nilai toleransi. Setelah iterasi ke-7 diperoleh nilai vektor x sebesar: 

    ⎡2,425476⎤

    x= ⎢3,573016⎥
    ⎣1,925954⎦

    Berdasarkan algoritma Iterasi Gauss-Seidel, kita dapat menyusun fungsi sebuah fungsi untuk melakukan iterasi Gauss-Seidel. Berikut sintaks yang digunakan: 

    gauss_seidel <- function(a, b, tol=1e-7, maxiter=100){
  n <- length(b)
  iter <- 0
  
  
  L <- U <- a
  L[upper.tri(a, diag=FALSE)] <- 0
  U[lower.tri(a, diag=TRUE)] <- 0
  Linv <- solve(L)
  
  x <- rep(0,n)
  x_new <- rep(tol, n)
  
  while(sqrt(sum(x_new-x)^2)>tol){
            if(iter>maxiter){
              warning("iterasi maksimum tercapai")
              break
            }
            x <- x_new
            x_new <- Linv %*% (b - U %*% x)
            iter <- iter+1
  }
      return(list(X = x_new, iter=iter))
}

    Contoh 6.12 Selesaikan sistem persamaan pada Contoh 6.10 menggunakan fungsi gauss_seidel()! Jawab: Penyelesaiansistem persamaan linier tersebut menggunakan fungsi gauss_seidel() disajikan pada sintaks berikut: 

    gauss_seidel(A,b)
    ## $X
##          [,1]
## [1,] 2.425476
## [2,] 3.573016
## [3,] 1.925954
## 
## $iter
## [1] 7
    1. Studi Kasus Aljabar linier banyak diaplikasikan baik dalam bidang engineering, fisika, sampai dengan statistika. Pada sub-chapter ini penulis akan menjelaskan penerapan aljabar linier pada metode kuadrat terkecil dan aliran massa dalam reaktor. Untuk penerapan lainnya pembaca dapat membaca buku lainnya terkait aljabar linier.
    1. Metode Kuadrat Terkecil Metode kuadrat terkecil merupakan salah satu aplikasi penerapan aljabar linier yang paling populer. Intuisi dibalik metode ini adalah bagaimana kita meminimalkan jarak antara sejumlah titik dengan garis regresi. Misalkan kita menggambarkan scatterplot antara dua buah variabel. Pola yang terbentuk dari plot tersebut adalah terjadi korelasi positif antara variabel pada sumbu xx dan sumbu yy. Kita ingin menggambarkan garis regresi terbaik yang dapat menangkap seluruh pola tersebut. Garis regresi terbaik terjadi ketika jumlah kuadrat jarak antara titik observasi dan garis regresi yang terbentuk seminimal mungkin. Untuk lebih memahaminya kita akan melakukan latihan menggunakan dataset trees yang berisi data hasil pengukuran kayu dari pohon yang ditebang. Pada dataset ini terdapat 31 observasi dan 3 buah kolom. Keterangan dari ketiga buah kolom tersebut adalah sebagai berikut: • Girth: diameter pohon dalam satuan inch. • Height: tinggi pohon dalam satuan feet. • Volume: volume kayu dalam satuan cubic feet. Untuk mengecek 6 observasi pertama dan struktur data, jalankan sintaks berikut: 
    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 ...

    Kita ingin membuat sebuah model linier untuk memprediksi Volume kayu berdasarkan variabel Girth dan Heiht atau volume sebagia fungsi dari variabel Girth dan Heiht. Kita dapat menuliskan relasi antara variabel volume sebagai fungsi dari variabel Girth dan Heiht menggunakan Persamaan (6.40). Volume=βgirthGirth+βheightHeight+β0(6.40 dimana β0β0 merupakan intersep persamaan regresi linier dan nilai ββ lainnya merupakan koefisien dari variabel Girth dan Heiht. Variabel Volume disebut sebagai variabel respon, sedangkan variabel Girth dan Heiht disebut sebagai variabel prediktor. Metode kuadrat terkecil berusaha memperoleh seluruh koefisien variabel dan intersep dari persamaan regresi linier. Berdasarkan yang telah penulis jelaskan garis regresi terbaik adalah garis yang memiliki nilai kuadrat terkecil jarak antara titik observasi dan garis regresi. Dasar dari metode kuadrat terkecil merupakan persamaan yang relatif sederhana yang ditunjukkan pada Persamaan (6.41). ATA=ATb(6.41) dimana bb merupakan vektor dari variabel respon (Volume) dan matrik AA merupakan matriks variabel prediktor (variabel Girth dan Heiht).

    Untuk menginputkan intercept kedalam persamaan linier kita perlu menmabhakan satu kolom di awal matriks AA yang berisi nilai 1. Berikut adalah sintaks yang digunakan untuk membentuk matriks AA: 

    # membentuk matriks A
pred <- cbind(intercept=1, Girth=trees$Girth, Height=trees$Height)
head(A)
    ##      [,1] [,2] [,3]
## [1,]   27    6   -1
## [2,]    6   15    2
## [3,]    1    1   54
    #Langkah selanjutnya adalah membentuk matriks b. Berikut adalah sintaks yang digunakan:
resp<- trees$Volume
head(resp)
    ## [1] 10.3 10.3 10.2 16.4 18.8 19.7
    #Untuk memperoleh koefisien β , kita dapat mencarinya dengan cara menyelesaikan Persamaan (6.41). Berikut adalah sintaks yang digunakan:
A <- t(pred) %*% pred
b <- t(pred) %*% resp

Ab <- cbind(A,b)
(x <- gauss_jordan(Ab))
    ##           intercept Girth Height            
## intercept         1     0      0 -57.9876589
## Girth             0     1      0   4.7081605
## Height            0     0      1   0.3392512

    Berdasarkan hasil yang diperoleh, persamaan linier yang terbentuk disajikan pada Persamaan (6.42).

    Volume=4.7081605Girth+0.3392512Height−57.9876589(6.42)

    Pembaca juga dapat menggunakan fungsi lain untuk memperoleh nilai koefisien tersebut, seperti: lu_solve()dansolve(). untuk fungsi jacobi() dan gauss_seidel(), kita harus pastikan syarat-syarat terkait metode tersebut. Berikut adalah contoh penyelesaian menggunakan sintaks lainnya: 

    # metode LU
lu_solve(A,b)
    ## $P
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1
## 
## $L
##          [,1]     [,2] [,3]
## [1,]  1.00000 0.000000    0
## [2,] 13.24839 1.000000    0
## [3,] 76.00000 1.054369    1
## 
## $U
##           intercept    Girth    Height
## intercept        31 410.7000 2356.0000
## Girth             0 295.4374  311.5000
## Height            0   0.0000  889.5641
## 
## $result
##             [,1]
## [1,] -57.9876589
## [2,]   4.7081605
## [3,]   0.3392512
    # fungsi solve()
solve(A,b)
    ##                  [,1]
## intercept -57.9876589
## Girth       4.7081605
## Height      0.3392512

    R juga menyediakan fungsi untuk membentuk model regresi linier. Fungsi yang digunakan adalah lm(). Berikut sintaks yang digunakan untuk membentuk model linier menggunakan fungsi lm(): 

    lm(Volume~Girth+Height, data=trees)
    ## 
## Call:
## lm(formula = Volume ~ Girth + Height, data = trees)
## 
## Coefficients:
## (Intercept)        Girth       Height  
##    -57.9877       4.7082       0.3393
    1. Aliran Massa Dalam Reaktor Pada sub-chapter ini penulis akan memberikan penerapan aljabar linier untuk menghitung konsentrasi suatu zat atau parameter lingkungan dalam reaktor yang saling terhubung. Pada contoh kasus kali ini diasumsikan terdapat lima buah reaktor yang saling terhubung satu sama lain sesuai Gambar 6.2. Debit air (m3detikm3detik) dan konsentrasi zat pencemar (mgm3mgm3) disajikan pula diagram alir tersebut. Diasumsikan kelima buah reaktor tersebut dalam kondisi steady dan volume reaktor diasumsikan sama. Kesetimbangan massa persatuan waktu dalam kondisi steady disajikan pada Persamaan (6.43).

    min=mout(6.43) QinCin=QoutCout(6.44)

    Berdasarkan Gambar 6.2, dapat dibentuk lima buah sistem persamaan linier. Persamaan linier yang terbentuk disajikan sebagai berikut: 6c1−c3=50 −3c1+3c2=0 −c2+9c3=160 −c2−8c3+11c4−2c5=0 −3c1−c2+4c5=06c1 Untuk menyelesaiakan sistem persamaan linier tersebut dan memperoleh nilai cc dari masing-masing reaktor, kiat perlu mengubahnya dulu kedalam bentuk matriks Ax=bAx=b. Berikut adalah matriks yang terbentuk:

    ⎡ 6 0 -1 0 0⎤ 50
    ⎢-3 3 0 0 0⎥ 0 ⎢ 0 -1 9 0 0⎥c= 160 ⎢ 0 -1 -8 11 -2⎥ 0 ⎢-3 -1 0 0 4⎥ 0 ⎣ ⎦

    Kita akan menyelesaikannya dengan menggunakan metode elminasi Gauss-Jordan, dekomposisi LU, iterasi Jacobi, dan iterasi Gauss-Seidel. Untuk dapat menyelesaikannya menggunakan metode-metode tersebut pada R, kita perlu membentuk matriksnya terlebih dahulu: 

    (A <- matrix(c(6,-3,0,0,-3,
              0,3,-1,-1,-1,
              -1,0,9,-8,0,
              0,0,0,11,0,
              0,0,0,-2,4),nrow=5))
    ##      [,1] [,2] [,3] [,4] [,5]
## [1,]    6    0   -1    0    0
## [2,]   -3    3    0    0    0
## [3,]    0   -1    9    0    0
## [4,]    0   -1   -8   11   -2
## [5,]   -3   -1    0    0    4
    (b <- c(50,0,160,0,0))
    ## [1]  50   0 160   0   0
    #Metode Eliminasi Gauss-Jordan
gauss_jordan(cbind(A,b))
    ##                       b
## [1,] 1 0 0 0 0 11.50943
## [2,] 0 1 0 0 0 11.50943
## [3,] 0 0 1 0 0 19.05660
## [4,] 0 0 0 1 0 16.99828
## [5,] 0 0 0 0 1 11.50943
    #Metode Dekomposisi LU
lu_solve(A,b)
    ## $P
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    0    0    0    0
## [2,]    0    1    0    0    0
## [3,]    0    0    1    0    0
## [4,]    0    0    0    1    0
## [5,]    0    0    0    0    1
## 
## $L
##      [,1]       [,2]       [,3] [,4] [,5]
## [1,]  1.0  0.0000000  0.0000000    0    0
## [2,] -0.5  1.0000000  0.0000000    0    0
## [3,]  0.0 -0.3333333  1.0000000    0    0
## [4,]  0.0 -0.3333333 -0.9245283    1    0
## [5,] -0.5 -0.3333333 -0.0754717    0    1
## 
## $U
##      [,1] [,2]          [,3] [,4] [,5]
## [1,]    6    0 -1.000000e+00    0    0
## [2,]    0    3 -5.000000e-01    0    0
## [3,]    0    0  8.833333e+00    0    0
## [4,]    0    0  0.000000e+00   11   -2
## [5,]    0    0  1.110223e-16    0    4
## 
## $result
## [1] 11.50943 11.50943 19.05660 16.99828 11.50943
    #Metode Iterasi Jacobi
jacobi(A,b, maxiter=100)
    ## $X
##          [,1]
## [1,] 11.50943
## [2,] 11.50943
## [3,] 19.05660
## [4,] 16.99828
## [5,] 11.50943
## 
## $iter
## [1] 17
    #Metode Iterasi Gauss-Seidel
gauss_seidel(A,b, maxiter=200)
    ## $X
##          [,1]
## [1,] 11.50943
## [2,] 11.50943
## [3,] 19.05660
## [4,] 16.99828
## [5,] 11.50943
## 
## $iter
## [1] 7

    Berdasarkan seluruh metode tersebut, diperoleh konsentrasi zat pencemar pada masing-masing reaktor adalah sebagai berikut: c1=11,50943mgm3 c2=11,50943mgm3 c3=19,05660mgm3 c4=16,99828mgm3 c5=11.50943mgm3

    REFERENSI: https://bookdown.org/moh_rosidi2610/Metode_Numerik/linearaljabar.html#jacobiiter