How to Do Gaussian Elimination in R (Do Not Use "Solve")

How to do Gaussian elimination in R (do not use solve)

I had to reorder you matrix A, don't know how do a generic code:

A <- matrix(c(2,-5,4,1,-4,6,1,-2.5,1),byrow=T,nrow=3,ncol=3)
b <- matrix(c(-3,5,10),nrow=3,ncol=1)
p <- nrow(A)
(U.pls <- cbind(A,b))

U.pls[1,] <- U.pls[1,]/U.pls[1,1]

for (i in 2:p){
 for (j in i:p) {
  U.pls[j, ] <- U.pls[j, ] - U.pls[i-1, ] * U.pls[j, i-1]
 }
 U.pls[i,] <- U.pls[i,]/U.pls[i,i]
}

for (i in p:2){
 for (j in i:2-1) {
  U.pls[j, ] <- U.pls[j, ] - U.pls[i, ] * U.pls[j, i]
 }
}
U.pls

EDIT:

A <- matrix(c(2,-5,4,1,-2.5,1,1,-4,6),byrow=T,nrow=3,ncol=3)
b <- matrix(c(-3,5,10),nrow=3,ncol=1)
p <- nrow(A)
(U.pls <- cbind(A,b))

U.pls[1,] <- U.pls[1,]/U.pls[1,1]

i <- 2
while (i < p+1) {
 j <- i
 while (j < p+1) {
  U.pls[j, ] <- U.pls[j, ] - U.pls[i-1, ] * U.pls[j, i-1]
  j <- j+1
 }
 while (U.pls[i,i] == 0) {
  U.pls <- rbind(U.pls[-i,],U.pls[i,])
 }
 U.pls[i,] <- U.pls[i,]/U.pls[i,i]
 i <- i+1
}
for (i in p:2){
 for (j in i:2-1) {
  U.pls[j, ] <- U.pls[j, ] - U.pls[i, ] * U.pls[j, i]
 }
}
U.pls

Gaussian Elimination not working

You say that your code creates the upper triangular matrix of A correctly, but it doesn't. Let me show you an example.

Let A and b be

A =

 3     2     3     4
 4     3     2     1
 1     0     4     0
 0     5     0     3

b =

 2
 4
 6
 7

If we run your code as it is and look at A and b we get

A =

3.0000    2.0000    3.0000    4.0000
1.3333    0.3333   -2.0000   -4.3333
0.3333   -2.0000   -1.0000  -10.0000
     0   15.0000  -30.0000 -232.0000

b =

7.0000
-203.0000
187.0000
7.0000

Which is neither a triangular matrix, nor the b that we expected. But if we modify slightly your program to be:

function [ x ] = Gauss( A,b )
n=length(b);
for j=1:n-1
    if A(j,j)==0
        break;
    end;

    for i=j+1:n
       f=A(i,j)/A(j,j);  %%Save the proportion between the rows in a
                         %%different variable outside the matrix, or
                         %%you will loose the value that was originally there

       b(i)=b(i)-f*b(j);   %%The operation has to be done in the row you are currently working

       for k=1:n    %%You have to make the operation in the full row,
                    %%not only in the remaining columns, also you can 
                    %%make this without a for loop using `:`
                    %%indexing, but if you dont know about it, 
                    %%leave as it is, it works
          A(i,k)=A(i,k)-f*A(j,k);
       end
    end
end
A
b
end

You get this result

A =

3.0000    2.0000    3.0000    4.0000
     0    0.3333   -2.0000   -4.3333
     0         0   -1.0000  -10.0000
     0         0         0 -232.0000

b =

2.0000
1.3333
8.0000
227.0000

Which is a upper triangular matrix, and the b we want. Hopefully you can take it from here, just as a reference the next steps should look like this

A =

1.0000    0.6667    1.0000    1.3333
     0    1.0000   -6.0000  -13.0000
     0         0    1.0000   10.0000
     0         0         0    1.0000

b =

0.6667
4.0000
-8.0000
-0.9784

then

 A =

 1     0     0     0
 0     1     0     0
 0     0     1     0
 0     0     0     1

b =

-1.1379
1.9871
1.7845
-0.9784

In which A is already an identity matrix, which mean that b is already our answer, that we can corroborate doing

A\b

ans=

-1.1379
1.9871
1.7845
-0.9784

J: Gauss-Jordan elimination

Thanks to Eelvex, who advised me to look in addons/math/misc/linear.ijs, I've concluded the task with this nice code:

   gj=: monad :0
   I=. i.#y
   for_i. I do. y=. y - (col - i=I) */ (i{y) % i{col=. i{"1 y end.
   )
   gj Ab
1 0 0 0 _1
0 1 0 0  3
0 0 1 0 _2
0 0 0 1  2

It has taken some time to understand verb pivot in linear.ijs - but pencil-paper method helps.

R Gaussian Elimination and qr factorization

When you read the help for qr you see that R uses a pivoted QR-decomposition.
So

str(b) 

gives

List of 4
 $ qr   : num [1:30, 1:6] -3.2292 0.218 0.0623 0.0371 0.302 ...
 $ rank : int 4
 $ qraux: num [1:6] 1.05 1.11 1.04 1.22 0 ...
 $ pivot: int [1:6] 1 3 5 6 2 4
 - attr(*, "class")= chr "qr"

Thus you need to apply pivot to a or the inverse of pivot to d to line the matrices up correctly. So

pivots <- b$pivot
d.ok <- d[,order(pivots)]
all.equal(a,d.ok)

gives

[1] TRUE

You can also do

a.p <- a[,pivots]
all.equal(a.p,d)

which also results in TRUE.

What is wrong with my Gauss-Jordan elimination?

Problem is solved. The question has been updated, among other with the new code (the old one is still available, to allow comparison). There were two bugs (the below "STEP XYZ" references the corresponding source code's STEP, not the steps mentionned on this StackOverflow question, which are presented a bit differently) :

1. The subtraction concerning the identity matrix didn't use identity matrix's coefficient (STEP 4). Bug fix : identity_matrix(id_line)(id_column) -= identity_matrix(id_last_pivot)(id_column) * tmp

2. Second, in STEP 3, I forgot to store the pivot in a temporary variable, in order to divide the both matrices (original and identity ones) with it. Without storing it, the value of the pivot changed after the division on the original matrix. Bug fix :

       val tmp = mutable_being_inversed_matrix(id_last_pivot)(general_id_column)
       mutable_being_inversed_matrix(id_last_pivot) = mutable_being_inversed_matrix(id_last_pivot).map(coefficient => coefficient / tmp)
       identity_matrix(id_last_pivot) = identity_matrix(id_last_pivot).map(coefficient => coefficient / tmp)
   

---

Related Topics