Implementing a Jacobi iterative method for Ax=b

classic Classic list List threaded Threaded
15 messages Options
Reply | Threaded
Open this post in threaded view
|

Implementing a Jacobi iterative method for Ax=b

Joza
Hi guys. I am trying to implement an iterative method for solving Ax=b using the Jacobi method. When I run the code with any tolerance it doesn't end, however. Any ideas why?

The matrix A is relatively sparse. Should this mean I should do row permutations on it before the iteration begins?
*********************************************************
function [x iter] = Jacobi(A, b, tolerance)

iter = 0;
n = size(A,1);
x_old = zeros(n,1);
converged = 0;

M = diag(diag(A));
N = M - A;

while ~converged

                x_new = M\(N*x_old + b);
                iter = iter + 1;
               
                if norm(b - A*x_new)/norm(b) < tolerance
                        x = x_new;
                        converged = 1;
                else
                        x_old = x_new;
                end
       
end
Reply | Threaded
Open this post in threaded view
|

Re: Implementing a Jacobi iterative method for Ax=b

Joza
I should add that I am using a zero starting guess, and the convergance test is testing the relative residual.
Reply | Threaded
Open this post in threaded view
|

Re: Implementing a Jacobi iterative method for Ax=b

Juan Pablo Carbajal-2
On Sun, Oct 28, 2012 at 5:40 PM, Joza <[hidden email]> wrote:

> I should add that I am using a zero starting guess, and the convergance test
> is testing the relative residual.
>
>
>
> --
> View this message in context: http://octave.1599824.n4.nabble.com/Implementing-a-Jacobi-iterative-method-for-Ax-b-tp4645833p4645834.html
> Sent from the Octave - General mailing list archive at Nabble.com.
> _______________________________________________
> Help-octave mailing list
> [hidden email]
> https://mailman.cae.wisc.edu/listinfo/help-octave

Your code has one error and one unnecessary calculation.

Error:
The step of the Jacobi iteration is
x_new = M\(b - N*x_old);

Unncesessary calculations
M is a diagonal matrix, therefore the inverse is just the inverse of
each diagonal element (take care of zero elements!). Improve the speed
of you code as follows

M = diag(A);

do not build the matrix, just keep the vector

x_new = (b - N*x_old)./M;

The method is not always convergent. Check that the input matrix
guarantees convergence for your tests. For real life usage, limit the
maximum numbe rof interations. When aborting because of max number of
iteration provde a warning and return enough information for the user
to take a decision.
http://en.wikipedia.org/wiki/Jacobi_method#Convergence

Hope that helps
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Implementing a Jacobi iterative method for Ax=b

Joza
Thanks for the input!

My method was correct, but you have informed me that I should be using the "correction form".

How can I permute an arbitrary matrix so that it's diagonal entries are not zero?
this will affect my solution I'm sure!
Reply | Threaded
Open this post in threaded view
|

Re: Implementing a Jacobi iterative method for Ax=b

Joza
Here is my latest code:

function [x, k] = Jacobi2(A, b, tol)

k = 0;
n = size(A,1);
x_old = zeros(n,1);
converged = 0;

for i=1:n
        if A(i,i) == 0.0
                W = A(dmperm(A), :)
                break
        else W = A;
        end
end

for i=1:n
        D_inv(i,i) = 1/W(i,i)
end

while ~converged

                x_new = x_old + D_inv*(b - W*x_old)   % Correction form
                k = k + 1;
               
                if norm(b - W*x_new)/norm(b) < tol
                        x = x_new
                        converged = 1;
                else
                        x_old = x_new;
                end
end
************************************************
I reorder the matrix using dmperm so that no diagonals are zero.

When i run this get the error:

warning: broken pipe -- some output may be lost

Whats going on!!??



Reply | Threaded
Open this post in threaded view
|

Re: Implementing a Jacobi iterative method for Ax=b

Sergei Steshenko





----- Original Message -----

> From: Joza <[hidden email]>
> To: [hidden email]
> Cc:
> Sent: Sunday, October 28, 2012 11:48 PM
> Subject: Re: Implementing a Jacobi iterative method for Ax=b
>
> Here is my latest code:
>
> function [x, k] = Jacobi2(A, b, tol)
>
> k = 0;
> n = size(A,1);
> x_old = zeros(n,1);
> converged = 0;
>
> for i=1:n
>     if A(i,i) == 0.0
>         W = A(dmperm(A), :)
>         break
>     else W = A;
>     end
> end   
>
> for i=1:n
>     D_inv(i,i) = 1/W(i,i)
> end
>
> while ~converged
>
>         x_new = x_old + D_inv*(b - W*x_old)   % Correction form
>         k = k + 1;
>        
>         if norm(b - W*x_new)/norm(b) < tol
>             x = x_new
>             converged = 1;
>         else
>             x_old = x_new;
>         end
> end
> ************************************************
> I reorder the matrix using dmperm so that no diagonals are zero.
>
> When i run this get the error:
>
> warning: broken pipe -- some output may be lost
>
> Whats going on!!??
>

"Whats going on!!??" - you forgot to put a semicolon at the end of

        W = A(dmperm(A), :)
            x = x_new
    D_inv(i,i) = 1/W(i,i)

lines.

Regards,
  Sergei.
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Implementing a Jacobi iterative method for Ax=b

Ozzy Lash
In reply to this post by Joza
On Sun, Oct 28, 2012 at 4:48 PM, Joza <[hidden email]> wrote:

> Here is my latest code:
>
> function [x, k] = Jacobi2(A, b, tol)
>
> k = 0;
> n = size(A,1);
> x_old = zeros(n,1);
> converged = 0;
>
> for i=1:n
>         if A(i,i) == 0.0
>                 W = A(dmperm(A), :)
>                 break
>         else W = A;
>         end
> end
>
> for i=1:n
>         D_inv(i,i) = 1/W(i,i)
> end
>
> while ~converged
>
>                 x_new = x_old + D_inv*(b - W*x_old)   % Correction form
>                 k = k + 1;
>
>                 if norm(b - W*x_new)/norm(b) < tol
>                         x = x_new
>                         converged = 1;
>                 else
>                         x_old = x_new;
>                 end
> end
> ************************************************
> I reorder the matrix using dmperm so that no diagonals are zero.
>
> When i run this get the error:
>
> warning: broken pipe -- some output may be lost
>
> Whats going on!!??
>
>

I'm guessing that the message comes because you had to hit control-c to stop it.

A few things:

The way you have left it, the W matrix still has the diagonal
elements.  Also your update is not quite rite, it isn't

x_new = x_old + D_inv*(b - W*x_old)

it should be:

x_new = D_inv*(b - W*x_old)

assuming you remove the diagonal elements of W.  I changed your script
to the following, and it seems to work for a matrix I tried:

function [x, k] = Jacobi2(A, b, tol)

k = 0;
n = size(A,1);
x_old = zeros(n,1);
converged = 0;

for i=1:n
        if A(i,i) == 0.0
                W = A(dmperm(A), :)
                break
        else W = A;
        end
end
R=W;
for i=1:n
        D_inv(i,i) = 1/W(i,i)
        W(i,i)=0
end

while ~converged

                x_new = D_inv*(b - W*x_old)   % Correction form
                k = k + 1;


                if norm(b - R*x_new)/norm(b) < tol
                        x = x_new
                        converged = 1;
                else
                        x_old = x_new;
                end
end
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Implementing a Jacobi iterative method for Ax=b

c.-2
In reply to this post by Juan Pablo Carbajal-2

On 28 Oct 2012, at 17:59, Juan Pablo Carbajal wrote:

> do not build the matrix, just keep the vector
> x_new = (b - N*x_old)./M;

Octave is clever enough to do that itself, diagonal matrices are a special class:

>> typeinfo (diag (diag (M)))
ans = diagonal matrix


and have their own methods.
c.
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Implementing a Jacobi iterative method for Ax=b

Ozzy Lash
Oh, and if you permute A, I think you will also need to permute b so
that they match, and then permute the result back to the original
form.  I didn't try a matrix with any 0s on the diagonal. I used:

A=[-500, 2, 3; 3, 200, 1; 4,5,600]
b=[1;2;3]
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Implementing a Jacobi iterative method for Ax=b

Joza
Since I permute A's rows using dmperm, how do I also permute b in the same way?
Reply | Threaded
Open this post in threaded view
|

Re: Implementing a Jacobi iterative method for Ax=b

Ozzy Lash
On Sun, Oct 28, 2012 at 6:58 PM, Joza <[hidden email]> wrote:
> Since I permute A's rows using dmperm, how do I also permute b in the same
> way?
>
>
>
> --
> View this message in context: http://octave.1599824.n4.nabble.com/Implementing-a-Jacobi-iterative-method-for-Ax-b-tp4645833p4645857.html
> Sent from the Octave - General mailing list archive at Nabble.com.


It has been a long time since I did anything like this, but since you
are only permuting the rows of A, I think you only need to permute the
rows of b in the same way, and x will still be ok.  If you were
permuting both rows and columns of A, you would need to permute the
rows of b in the same way that you permute the rows of A, and
"unpermute" the x vector to undo the permutation of the columns of A.
That is only from memory and a little Sunday evening thinking, so I
could have that mixed up.

Since you only need to permute b, you could, where you have:

 W = A(dmperm(A), :);

You could put:

P=dmperm(A);
W=A(P,:);
bp=b(p);

And then do your calculations with bp instead of b
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Implementing a Jacobi iterative method for Ax=b

Joza
So I've taken all this on board, and finally got an algorithm which as far as I am concerned SHOULD work...

Please note that I'm using the "correction form" of the Jacobi algorithm: x_new = x_old + D_inv*(b - A*x_old);
where A is from Ax=b (but possibly permuted).

You can derive it from the usual x_new = D_inv*( b - L*x_old - U*x_old ) where A = D + L + U.

So here is the code:
************************************************************************************************************************************************
function [x k] = Jacobi(A, b, tol)

k = 0;
n = size(A,1);
x_old = zeros(n,1);
converged = 0;

W = A(dmperm(A), :); % dmperm permutes rows of A so that no diagonal elements are zero.
v = b(dmperm(A),1);  % must permute b the same way so that x is unchanged.

for i=1:n
        D_inv(i,i) = 1/W(i,i);
end

while ~converged

                x_new = x_old + D_inv*(v - W*x_old);
                k = k + 1;
               
                if norm(v - W*x_new, inf)/norm(v, inf) < tol
                        x = x_new;
                        converged = 1;
                else
                        x_old = x_new;
                end
       
end
*******************************************************************************************************************8
I have been given a starting A and b and told to use the zero starting approximation.

When I run the algorithm, I do not get a return x or k, all the terminal gives me is ans = [NaN, NaN, NaN ...]

So clearly this is running into difficulty somewhere and not finishing correctly, but I really can't see where. Everything should work as far as I can see!!!
Reply | Threaded
Open this post in threaded view
|

Re: Implementing a Jacobi iterative method for Ax=b

Joza
Here are A and b if you want to test it out!

A = [ 3 12 0 -1 0 0; 4 0 31 1 0 0; 2 1 0 0 17 -3; 27 2 0 0 0 1; 0 0 0 -1 1 11; 0 0 0 24 -1 0 ]

(6x6)

b = [ 39 117 12 98 14 55]^T

(6x1)
Reply | Threaded
Open this post in threaded view
|

Re: Implementing a Jacobi iterative method for Ax=b

Juan Pablo Carbajal-2
On Mon, Oct 29, 2012 at 2:16 PM, Joza <[hidden email]> wrote:

> Here are A and b if you want to test it out!
>
> A = [ 3 12 0 -1 0 0; 4 0 31 1 0 0; 2 1 0 0 17 -3; 27 2 0 0 0 1; 0 0 0 -1 1
> 11; 0 0 0 24 -1 0 ]
>
> (6x6)
>
> b = [ 39 117 12 98 14 55]^T
>
> (6x1)
>
>
>
> --
> View this message in context: http://octave.1599824.n4.nabble.com/Implementing-a-Jacobi-iterative-method-for-Ax-b-tp4645833p4645868.html
> Sent from the Octave - General mailing list archive at Nabble.com.
> _______________________________________________
> Help-octave mailing list
> [hidden email]
> https://mailman.cae.wisc.edu/listinfo/help-octave

Hi,

Use Agora to post your code, it is easier to track the changes
http://agora.octave.org/

I posted an improved version of your file.
http://agora.octave.org/snippet/JZaP/

You must always check if the algorithm will converge, in order to know
if you should expect an answer or not. Try
demo jacobi_iter

Cheers
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Implementing a Jacobi iterative method for Ax=b

Ozzy Lash
On Tue, Oct 30, 2012 at 6:10 AM, Juan Pablo Carbajal
<[hidden email]> wrote:

> On Mon, Oct 29, 2012 at 2:16 PM, Joza <[hidden email]> wrote:
>> Here are A and b if you want to test it out!
>>
>> A = [ 3 12 0 -1 0 0; 4 0 31 1 0 0; 2 1 0 0 17 -3; 27 2 0 0 0 1; 0 0 0 -1 1
>> 11; 0 0 0 24 -1 0 ]
>>
>> (6x6)
>>
>> b = [ 39 117 12 98 14 55]^T
>>
>> (6x1)
>>
>>
>>
>> --
>> View this message in context: http://octave.1599824.n4.nabble.com/Implementing-a-Jacobi-iterative-method-for-Ax-b-tp4645833p4645868.html
>> Sent from the Octave - General mailing list archive at Nabble.com.
>> _______________________________________________
>> Help-octave mailing list
>> [hidden email]
>> https://mailman.cae.wisc.edu/listinfo/help-octave
>
> Hi,
>
> Use Agora to post your code, it is easier to track the changes
> http://agora.octave.org/
>
> I posted an improved version of your file.
> http://agora.octave.org/snippet/JZaP/
>
> You must always check if the algorithm will converge, in order to know
> if you should expect an answer or not. Try
> demo jacobi_iter
>
> Cheers
> _______________________________________________
> Help-octave mailing list
> [hidden email]
> https://mailman.cae.wisc.edu/listinfo/help-octave

Juan,

Great pointer to agora.octave.org.  I hadn't heard of it before.

I just posted an update to your script that attempts to create a
better permutation matrix.  It seems to me that to make a matrix
diagonal dominant, it is a necessary condition that the largest
magnitude value in any row be on the diagonal.  I added an initial
permutation "guess" p that is simply the indices of max(abs(A)).  I
then test to see if all of the rows are unique, and, if not fall back
to using dmperm.

This will definitely help with the test case that the initial poster was using.

Bill
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave