Jacobi Iteration With Ruby
Wed 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.
For example, let’s assume we’ve been given this matrix A:
We will then zero-out everything above and including the diagonal entries. This will give us our new L matrix.
Likewise, we’ll remove all entries below, and including the diagonal from matrix A to give us our matrix U.
Finally, we’ll zero-out all non-diagonal entries, to yield our new matrix D.
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:
And since D is a diagonal matrix, the above equation can be rewritten in terms of the components of A and b.
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.