Gaussian Elimination in Matlab

Question

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.

Answer 1

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.

Answer 2

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(n³) algorithms. If any one approach is better than another depends on your particular situation and is something you would need to investigate more.