2019年9月11日星期三

A function for Gauss-Seidel Iterative methods

Gauss-Seidel Iterative method is essentially the back fitting method for GAM model.

Such as to solve equations

 5x1- 2x2+3x3=-1
-3x1+9x2+  x3= 2
 2x1 - x2 - 7x3= 3


library(MASS) # package needed for ginv() function
Gauss_Seidel<-function(A,b,x0){
 
B<-t(A)%*%A
C<-t(A)%*%b

iter<-100          # number of Gauss-Seidel iterations to run,100 iterations
L<-lower.tri(B)*B# lower triang. B,lower.tri() Returns a matrix of logicals,therefore need to *B
U<-upper.tri(B)*B  # upper triang. B
D<-diag(diag(B))   # diag of A
n<-1
while(n<=iter){
  print(x0)
 
  # Gauss-Seidel formula
  x1<-(ginv((L+D)))%*%((-U%*%x0)+C)
  x0<-x1
  n=n+1
 
}

return(x0)

}

A<-matrix(c(5,-3,2,-2,9,-1,3,1,-7),ncol=3,nrow=3)
b<-matrix(c(-1,2,3),ncol=1,nrow=3)
x0<-t(t(c(0,0,0)))

Gauss_Seidel(A,b,x0)

After 10 iterations:


> Gauss_Seidel(A,b,x0)
     [,1]
[1,]    0
[2,]    0
[3,]    0
           [,1]
[1,] -0.1315789
[2,]  0.1380049
[3,] -0.4007323
            [,1]
[1,] -0.01103352
[2,]  0.23926763
[3,] -0.41380921
            [,1]
[1,]  0.09220576
[2,]  0.28760601
[3,] -0.41850252
           [,1]
[1,]  0.1415692
[2,]  0.3105375
[3,] -0.4207159
           [,1]
[1,]  0.1649876
[2,]  0.3214148
[3,] -0.4217656
           [,1]
[1,]  0.1760960
[2,]  0.3265744
[3,] -0.4222636
           [,1]
[1,]  0.1813651
[2,]  0.3290218
[3,] -0.4224998
           [,1]
[1,]  0.1838645
[2,]  0.3301827
[3,] -0.4226118
           [,1]
[1,]  0.1850501
[2,]  0.3307334
[3,] -0.4226650
           [,1]
[1,]  0.1856124
[2,]  0.3309946
[3,] -0.4226902

#Check with solve function.
solve(A,b)