profile image

Jacobi Iteration With Ruby

April 28 2010

When solving linear algebraic systems, it is often beneficial to use iterative methods. Iterative solvers are particularly useful when a matrix is diagonally dominant. One such iterative method that is particularly useful is Jacobi Iteration.

Given a matrix A, you first decompose it into 3 separate matrices: L, D, and U. These are the lower, diagonal, and upper sections of the matrix, respectively.
Jacobi Decomposition

For example, let’s assume we’ve been given this matrix A:
Example Matrix

We will then zero-out everything above and including the diagonal entries. This will give us our new L matrix.
Example L Matrix

Likewise, we’ll remove all entries below, and including the diagonal from matrix A to give us our matrix U.
Example U Matrix

Finally, we’ll zero-out all non-diagonal entries, to yield our new matrix D.
Example D Matrix

As stated previously, since A=L+D+U, we can deduce that Ax=Dx+(L+U)x. And since Ax=b, we can conclude that Dx=b-(L+U)x, which is the basis of Jacobi’s Method.

Given an estimate vector (x) (which is generally initialized to 0 if no educated guess can be provided), the next estimate is found by:
Basic Jacobi Form

And since D is a diagonal matrix, the above equation can be rewritten in terms of the components of A and b.
Jacobi Equation

So now that we understand how Jacobi Iteration works, we can write a program to do it for us! Ruby is not necessarily the best language to do this in, as these systems can often times be quite large. Compiled languages such as C++ are going to converge much faster.

To start out, we’re going to write a method that takes in our Matrix A, our right hand size b, the number of times we’d like to iterate, and our initial guess.

def jacobi_iteration(A, b, iterations, guess |= 0)

From here, we will need to fire up our estimate vector x. We’ll initialize its first element to be the initial “guess,” and the remaining entries to be 0.

  # Initialize
  x = Array.new(A.size)
  x[0] = guess
  1.upto(A.size) do |i|
    x[i] = 0
  end

Finally, we implement the equation above as nested loops to iteratively try to converge to a solution.

  0.upto(iterations-1) do |its|
    0.upto(A.size) do |i|
      sum = 0
      0.upto(A.size) do |j|
        unless j==i
          sum = A[i][j] * x[j]
        end
      end
      x[i] = 1.0/A[i][i] * (b[i] - sum)
    end
  end
end

And there you have it! If you're doing this on large matrices, I would definitely suggest using a different language, such as C++ to handle this, since Ruby can not only be a bit slower, but it's also quite a memory hog when dealing with so many objects.