patchGNU Octave - Patches: patch #9282, add minres for sparse linear...

 
 

patch #9282: add minres for sparse linear systems

Submitted by:  Xie Rui <xierui>
Submitted on:  Wed 08 Mar 2017 01:13:15 PM UTC  
 
Category: Core : new functionPriority: 5 - Normal
Status: NonePrivacy: Public
Assigned to: NoneOpen/Closed: Open

Add a New Comment(Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

(Jump to the original submission Jump to the original submission)

Fri 17 Mar 2017 03:10:41 AM UTC, comment #14:

Now it always uses a split preconditioner, and the warning about too fast convergence was removed.

(file #40018)

Xie Rui <xierui>
Thu 16 Mar 2017 09:13:04 AM UTC, comment #13:

If you look at algorithm 9.1 on Saad's book, you see that the split preconditioning for the conjugate gradient is realised at step 6. And it does not matter whether the preconditioner M is M1*M2. The only thing to do is to multiply the residual by inv(M). I see that Matlab always use a split preconditioner:

I would remove the warning about too fast convergence. And some tests do not work for me.

Marco Caliari <caliari>
Project Member
Thu 16 Mar 2017 03:47:12 AM UTC, comment #12:

I removed the try...catch inside the iteration and the check of ishermitian.

I changed the stagnation criterion to norm(resvec(k+1)-resvec(k)) <= eps * norm(resvec(k+1)).

I changed the use of preconditioner. Now it uses m1 and m2 as a split preconditioner. However, it still uses m1 as a left preconditioner when m2 does not exist.

(file #40010)

Xie Rui <xierui>
Wed 15 Mar 2017 07:59:22 AM UTC, comment #11:

You can remove the try..catch inside the iteration: if the preconditioner is singular in the ## Initiation stage, the code will return. Cristiano used ishermitian only in some tests, I think you can remove from the code.

I think you should change the stagnation criterion. It is almost impossible that two consecutive iterations coincide (with ==). You should use something like norm(resvec(k+1)-resvec(k)) <= eps * norm(resvec(k+1)).

I see you use a left preconditioner, while Matlab a split one.

Marco Caliari <caliari>
Project Member
Wed 15 Mar 2017 02:00:25 AM UTC, comment #10:

Re: Comment #5

1. Solved by Comment #7.

2. I learned from Cristiano, and added try...catch to detect singular preconditioner matrix. Matlab checks it, too. For example (it is similar to a test from pcg of Cristiano),

Then flag = 2, which means that Matlab detected singular preconditioner matrix.

However, Matlab does not check whether A is symmetric or whether A is hermitian. Because pcg of Cristiano checks whether A is hermitian, my minres checks it.

3. I fixed it.

4. Matlab accepts complex matrix and works on it. Test code:

If allconverge = 1, then Matlab minres works for all the A and b above.

(file #40006)

Xie Rui <xierui>
Tue 14 Mar 2017 01:40:22 AM UTC, comment #9:

Hi Nicholas,

Thank you for the test.
In Comment #8, A is symmetric but not hermitian. I think we should test for hermitian A as well. For example, I substituted A into the matlab example in Comment #1.

I tried it with my minres. It converged after 26 iterations, and the relative residual is 8.4181e-011.
In addition, it failed to converge as for the test in comment #8.

Xie Rui <xierui>
Mon 13 Mar 2017 03:27:21 PM UTC, comment #8:

re: questions from Comment #5, I'm not sure how to test 2, but I made a lousy complex A to substitute into the matlab example in Comment #1. It failed to converge, but it didn't produce any errors about complex input. Barring any other test, I'll assume A should be allowed to be complex.

Nicholas Jankowski <nrjank>
Mon 13 Mar 2017 10:55:08 AM UTC, comment #7:

Dear Marco,

I looked at formulas (5.13-5.16) and (7.4), but I think it is difficult to use them when there is a preconditioner matrix. Therefore, I adopted another method. I accumulated the residual in every iteration just like x, to avoid matrix-vector multiplication per step.
Moreover, I used this method and formula (7.1) to calculate resveccg, i.e. the norm of r^c in Paige&Saunders. By the way, my x is x^M, and my resvec is the norm of r^M in Paige&Saunders.

(file #39988)

Xie Rui <xierui>
Sun 12 Mar 2017 07:51:57 PM UTC, comment #6:

If you give me a few test scripts I can reply with the MATLAB output. If you brainstorm a large number of tests I could turn them into some test code blocks for you to debug against.

Nicholas Jankowski <nrjank>
Sat 11 Mar 2017 04:28:01 PM UTC, comment #5:

1. Could it be what given in formula (5.17) of Paige&Saunders?
Moreover, it seems that you compute the true residual (b-A*x), which requires an additional matrix-vector per step. I think you should use formulas (5.13)--(5.16).

2. If by "ill-conditioned" you mean that octave gives the warning "matrix is singular...", I think Cristiano did something with try... catch.
You could use issymmetric. Does matlab check it? For instance, there is not such a check in pcg.

4. Does matlab accept complex matrices? Does it work with complex matrices?

Marco Caliari <caliari>
Project Member
Sat 11 Mar 2017 12:52:26 PM UTC, comment #4:

There are still a lot of problems now:

  1. There is an output that matlab minres has but I do not know its meaning. This is resveccg. Therefore, I didn't add it as an output.
  2. My minres has 5 kinds of flags now, but I doubt its performance. Is there good methods to detect ill-conditioned matrix, when the information of this matrix is given by the form of matrix or functions which return the results of applying the inverse? And what about detecting matrix that is not sysmmetric or hermitian? Need tolerance here?
  3. My minres requires that if m1 and m2 both exist, then they should be both matrices or both function handles. However, in function pcg, m1 and m2 may be of different types.
  4. The paper I based on (matlab bases on it too according to the docummentation) says the matrix A should be real. I have not checked carefully whether it works with complex matrices, although it passed some tests.

I think maybe it is far from being a new function in Octave now, and I will continue to do the remaining work.

Xie Rui <xierui>
Sat 11 Mar 2017 09:37:50 AM UTC, comment #3:

I uploaded a new m-file. It is different from the old one in these aspects:

  1. I added some tests, including the suggestion of Nicholas Jankowski, modified tests from pcg.m and some tests for hermitian but NOT positive definite matrices.
  2. I fixed some bugs, such as maxit and the using of preconditioner matrix.
  3. I deleted the errors and call print_usage(). (I am not sure whether I have done it correctly.)
  4. I added the description of usage.

(file #39970)

Xie Rui <xierui>
Thu 09 Mar 2017 11:14:11 AM UTC, comment #2:

Thank you very much for your advices!
I ran the example with a slightly different code in matlab. That is, use the complete form of the outputs.
The result of relres (relative residual) and iter (iteration) is
not 49. It is strange.

For my m-file, I try to set maxit as 100. That is
Then the result is
and
Note that it converged at iteration 50 to a solution with relative residul 4.5e-014.
There is a bug about using maxit in my code. I am fixing it, and I will do others that you mentioned.

Xie Rui <xierui>
Wed 08 Mar 2017 04:28:41 PM UTC, comment #1:

Nice work. some comments:

are there any tests you could add to the end of the m-file, perhaps to verify expected input/output handling, errors, maybe verifying results as expected from the matlab examples?

I did run the example given on the referenced matlab help page

the output x is more or less equal to ones(100,1), with error on the order of 1e-14.

running the same on your minres produces:

Since it's an example we can use to test algorithm compatibility, I would suggest adding the following test block to the end of the file in addition to any others to catch input/output, format, etc.

Some of your errors related to providing incorrect inputs (nargin<2 test for example) should call print_usage()

Regarding usage, I don't know how mandatory this is, but usually there's a separate line for each valid calling method. so a line for minres(A,b), then minres(A,b,tol), etc.

Last, a complete patch should also update the scripts/sparse/module.mk file, the documentation file (I guess under doc/interpreter/sparse.txi), and add the new function to /NEWS

Nicholas Jankowski <nrjank>
Wed 08 Mar 2017 01:13:15 PM UTC, original submission:

It solves the linear system of equations A*x = b by means of the minimum residual method.
Refenrences:
1. Matlab documentation of minres
2. C. C. PAIGE and M. A. SAUNDERS, @cite{Solution of Sparse Indefinite Systems of Linear Equations},
SIAM J. Numer. Anal., 1975. (the minimum residual method)

Xie Rui <xierui>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach File(s):
   
   
Comment:
   

Attached Files
file #40018:  minres.m added by xierui (20KiB - text/plain)
file #40010:  minres.m added by xierui (21KiB - text/plain)
file #40006:  minres.m added by xierui (21KiB - text/plain)
file #39988:  minres.m added by xierui (21KiB - text/plain)
file #39970:  minres.m added by xierui (21KiB - text/plain)
file #39940:  minres.m added by xierui (10KiB - text/plain - add minres for sparse linear systems)
file #39941:  minres.patch added by xierui (11KiB - application/octet-stream - add minres for sparse linear systems)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by caliari (Posted a comment)
  • -unavailable- added by nrjank (Posted a comment)
  • -unavailable- added by xierui (Submitted the item)
  • -unavailable- added by xierui ([GSoC 2017 Application] add minres function for sparse linear systems)
  • -unavailable- added by xierui ([GSoC 2017 Application] add minres function for sparse linear systems)
  •  

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only logged-in users can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 9 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Fri 17 Mar 2017 03:10:41 AM UTCxieruiAttached File-=>Added minres.m, #40018
    Thu 16 Mar 2017 03:47:12 AM UTCxieruiAttached File-=>Added minres.m, #40010
    Wed 15 Mar 2017 02:00:25 AM UTCxieruiAttached File-=>Added minres.m, #40006
    Mon 13 Mar 2017 10:55:08 AM UTCxieruiAttached File-=>Added minres.m, #39988
    Sat 11 Mar 2017 09:37:50 AM UTCxieruiAttached File-=>Added minres.m, #39970
    Wed 08 Mar 2017 01:13:15 PM UTCxieruiAttached File-=>Added minres.m, #39940
      Attached File-=>Added minres.patch, #39941
      Carbon-Copy-=>Added -unavailable-
      Carbon-Copy-=>Added -unavailable-

    Back to the top


    Powered by Savane 3.1-cleanup1