CODEWITHSUNDEEP
matlabgaussian
sosruko
sosruko

Reputation: 1079

Gaussian Elimination in Matlab

I am using the matlab code from this book: http://books.google.com/books/about/Probability_Markov_chains_queues_and_sim.html?id=HdAQdzAjl60C Here is the Code:

    function [pi] = GE(Q)

    A = Q';
    n = size(A);
    for i=1:n-1
      for j=i+1:n
         A(j,i) = -A(j,i)/A(i,i);
      end
         for j =i+1:n
            for k=i+1:n
        A(j,k) = A(j,k)+ A(j,i) * A(i,k);
         end
      end
      end

     x(n) = 1;
      for i = n-1:-1:1
        for j= i+1:n
          x(i) = x(i) + A(i,j)*x(j);
        end
       x(i) = -x(i)/A(i,i);
      end

      pi = x/norm(x,1);

Is there a faster code that I am not aware of? I am calling this functions millions of times, and it takes too much time.

Upvotes: 5

Views: 97029

Answers (7)

Dominiksr
Dominiksr

Reputation: 78

I think you can use the matlab function rref:

[R,jb] = rref(A,tol)

It produces a matrix in reduced row echelon form. In my case it wasn't the fastest solution. The solution below was faster in my case by about 30 percent.

function C = gauss_elimination(A,B)
i = 1; % loop variable
X = [ A B ]; 
[ nX mX ] = size( X); % determining the size of matrix  

while i <= nX % start of loop 
    if X(i,i) == 0 % checking if the diagonal elements are zero or not
        disp('Diagonal element zero') % displaying the result if there exists zero 
        return
    end
    X = elimination(X,i,i); % proceeding forward if diagonal elements are non-zero
    i = i +1;
end

C = X(:,mX);

function X = elimination(X,i,j)
% Pivoting (i,j) element of matrix X and eliminating other column
% elements to zero
[ nX mX ] = size( X); 
a = X(i,j);
X(i,:) = X(i,:)/a;

for k =  1:nX % loop to find triangular form 
    if k == i
        continue
    end
    X(k,:) = X(k,:) - X(i,:)*X(k,j); 
end

Upvotes: 0

Feras Alsheet
Feras Alsheet

Reputation: 11

Let's assume Ax=d Where A and d are known matrices. We want to represent "A" as "LU" using "LU decomposition" function embedded in matlab thus: LUx = d This can be done in matlab following: [L,U] = lu(A) which in terms returns an upper triangular matrix in U and a permuted lower triangular matrix in L such that A = LU. Return value L is a product of lower triangular and permutation matrices. (https://www.mathworks.com/help/matlab/ref/lu.html)

Then if we assume Ly = d where y=Ux. Since x is Unknown, thus y is unknown too, by knowing y we find x as follows: y=L\d; x=U\y

and the solution is stored in x.

This is the simplest way to solve system of linear equations providing that the matrices are not singular (i.e. the determinant of matrix A and d is not zero), otherwise, the quality of the solution would not be as good as expected and might yield wrong results.

if the matrices are singular thus cannot be inversed, another method should be used to solve the system of the linear equations.

Upvotes: 1

achini
achini

Reputation: 449

function x = naiv_gauss(A,b);
n = length(b); x = zeros(n,1);
for k=1:n-1 % forward elimination
      for i=k+1:n
           xmult = A(i,k)/A(k,k);
           for j=k+1:n
             A(i,j) = A(i,j)-xmult*A(k,j);
           end
           b(i) = b(i)-xmult*b(k);
      end
end
 % back substitution
x(n) = b(n)/A(n,n);
for i=n-1:-1:1
   sum = b(i);
   for j=i+1:n
     sum = sum-A(i,j)*x(j);
   end
   x(i) = sum/A(i,i);
end
end

Upvotes: 5

ntox
ntox

Reputation: 11

function Sol = GaussianElimination(A,b)


[i,j] = size(A);

for j = 1:i-1


    for i = j+1:i


        Sol(i,j) = Sol(i,:) -( Sol(i,j)/(Sol(j,j)*Sol(j,:)));


    end


end


disp(Sol);


end

Upvotes: 0

Sean
Sean

Reputation: 1637

For the naive approach (aka without row swapping) for an n by n matrix:

function A = naiveGauss(A)

% find's the size

n = size(A);
n = n(1);

B = zeros(n,1);

% We have 3 steps for a 4x4 matrix so we have
% n-1 steps for an nxn matrix
for k = 1 : n-1     
    for i = k+1 : n
        % step 1: Create multiples that would make the top left 1
        % printf("multi = %d / %d\n", A(i,k), A(k,k), A(i,k)/A(k,k) )
        for j = k : n
            A(i,j) = A(i,j) - (A(i,k)/A(k,k)) *  A(k,j);
        end
        B(i) = B(i) - (A(i,k)/A(k,k))  * B(k);
    end
end

Upvotes: 0

David Alber
David Alber

Reputation: 18111

Unless you are specifically looking to implement your own, you should use Matlab's backslash operator (mldivide) or, if you want the factors, lu. Note that mldivide can do more than Gaussian elimination (e.g., it does linear least squares, when appropriate).

The algorithms used by mldivide and lu are from C and Fortran libraries, and your own implementation in Matlab will never be as fast. If, however, you are determined to use your own implementation and want it to be faster, one option is to look for ways to vectorize your implementation (maybe start here).

One other thing to note: the implementation from the question does not do any pivoting, so its numerical stability will generally be worse than an implementation that does pivoting, and it will even fail for some nonsingular matrices.

Different variants of Gaussian elimination exist, but they are all O(n3) algorithms. If any one approach is better than another depends on your particular situation and is something you would need to investigate more.

Upvotes: 8

Darren Engwirda
Darren Engwirda

Reputation: 7015

MATLAB has a whole set of built-in linear algebra routines - type help slash, help lu or help chol to get started with a few of the common ways to efficiently solve linear equations in MATLAB.

Under the hood these functions are generally calling optimised LAPACK/BLAS library routines, which are generally the fastest way to do linear algebra in any programming language. Compared with a "slow" language like MATLAB it would not be unexpected if they were orders of magnitude faster than an m-file implementation.

Hope this helps.

Upvotes: 9

Related Questions