Quantcast

[patch #9282] add minres for sparse linear systems

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
15 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[patch #9282] add minres for sparse linear systems

Marco Caliari-5
URL:
  <http://savannah.gnu.org/patch/?9282>

                 Summary: add minres for sparse linear systems
                 Project: GNU Octave
            Submitted by: xierui
            Submitted on: Wed 08 Mar 2017 01:13:15 PM UTC
                Category: Core : new function
                Priority: 5 - Normal
                  Status: None
                 Privacy: Public
             Assigned to: None
        Originator Email:
             Open/Closed: Open
         Discussion Lock: Any

    _______________________________________________________

Details:

It solves the linear system of equations A*x = b by means of the minimum
residual method.
Refenrences:
1. Matlab documentation of minres
<https://www.mathworks.com/help/matlab/ref/minres.html>
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)



    _______________________________________________________

File Attachments:


-------------------------------------------------------
Date: Wed 08 Mar 2017 01:13:15 PM UTC  Name: minres.m  Size: 10kB   By: xierui
add minres for sparse linear systems
<http://savannah.gnu.org/patch/download.php?file_id=39940>
-------------------------------------------------------
Date: Wed 08 Mar 2017 01:13:15 PM UTC  Name: minres.patch  Size: 11kB   By:
xierui
add minres for sparse linear systems
<http://savannah.gnu.org/patch/download.php?file_id=39941>

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/patch/?9282>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[patch #9282] add minres for sparse linear systems

Marco Caliari-5
Follow-up Comment #1, patch #9282 (project octave):

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

n = 100; on = ones(n,1);
A = spdiags([-2*on 4*on -2*on],-1:1,n,n);
b = sum(A,2);
tol = 1e-10;
maxit = 50;
M1 = spdiags(4*on,0,n,n);

x = minres(A,b,tol,maxit,M1);
minres converged at iteration 49 to a solution with relative
residual 4.7e-014


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:

>> x'
ans =
 Columns 1 through 8:
   0.97030   0.94062   0.91098   0.88142   0.85195   0.82260   0.79338  
0.76433
 Columns 9 through 16:
   0.73547   0.70681   0.67839   0.65023   0.62234   0.59476   0.56750  
0.54059
 Columns 17 through 24:
   0.51406   0.48792   0.46220   0.43692   0.41211   0.38779   0.36398  
0.34071
 Columns 25 through 32:
   0.31800   0.29586   0.27434   0.25344   0.23320   0.21363   0.19476  
0.17661
 Columns 33 through 40:
   0.15921   0.14257   0.12673   0.11171   0.09752   0.08419   0.07175  
0.06022
 Columns 41 through 48:
   0.04962   0.03998   0.03131   0.02365   0.01701   0.01142   0.00690  
0.00347
 Columns 49 through 56:
   0.00116   0.00000   0.00000   0.00116   0.00347   0.00690   0.01142  
0.01701
 Columns 57 through 64:
   0.02365   0.03131   0.03998   0.04962   0.06022   0.07175   0.08419  
0.09752
 Columns 65 through 72:
   0.11171   0.12673   0.14257   0.15921   0.17661   0.19476   0.21363  
0.23320
 Columns 73 through 80:
   0.25344   0.27434   0.29586   0.31800   0.34071   0.36398   0.38779  
0.41211
 Columns 81 through 88:
   0.43692   0.46220   0.48792   0.51406   0.54059   0.56750   0.59476  
0.62234
 Columns 89 through 96:
   0.65023   0.67839   0.70681   0.73547   0.76433   0.79338   0.82260  
0.85195
 Columns 97 through 100:
   0.88142   0.91098   0.94062   0.97030


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.


## test for algorithm accuracy and compatibility from matlab doc example
%!test
%! n = 100;
%! on = ones (n, 1);
%! A = spdiags ([-2*on 4*on -2*on], -1:1, n, n);
%! b = sum (A, 2);
%! tol = 1e-10;
%! maxit = 50;
%! M1 = spdiags (4*on, 0, n, n);
%! x = minres (A, b, tol, maxit, M1);
%! assert (size (x), [100, 1]);
%! assert (x,ones(100,1),1e-13);


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


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/patch/?9282>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[patch #9282] add minres for sparse linear systems

Marco Caliari-5
Follow-up Comment #2, patch #9282 (project octave):

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.

n = 100; on = ones(n,1);
A = spdiags([-2*on 4*on -2*on],-1:1,n,n);
b = sum(A,2);
tol = 1e-10;
maxit = 50;
M1 = spdiags(4*on,0,n,n);

[x,flag,relres,iter,resvec,resveccg] = minres(A,b,tol, maxit,M1)

The result of relres (relative residual) and iter (iteration) is

relres =

   4.6537e-14
iter =

    50

not 49. It is strange.

For my m-file, I try to set maxit as 100. That is

n = 100; on = ones(n,1);
A = spdiags([-2*on 4*on -2*on],-1:1,n,n);
b = sum(A,2);
tol = 1e-10;
maxit = 100;
M1 = spdiags(4*on,0,n,n);

[x,flag,relres,iter,resvec] = minres(A,b,tol,maxit,M1)

Then the result is

>> x'
ans =

 Columns 1 through 13:

   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000  
1.00000   1.00000   1.00000   1.00000   1.00000   1.00000

 Columns 14 through 26:

   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000  
1.00000   1.00000   1.00000   1.00000   1.00000   1.00000

 Columns 27 through 39:

   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000  
1.00000   1.00000   1.00000   1.00000   1.00000   1.00000

 Columns 40 through 52:

   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000  
1.00000   1.00000   1.00000   1.00000   1.00000   1.00000

 Columns 53 through 65:

   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000  
1.00000   1.00000   1.00000   1.00000   1.00000   1.00000

 Columns 66 through 78:

   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000  
1.00000   1.00000   1.00000   1.00000   1.00000   1.00000

 Columns 79 through 91:

   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000  
1.00000   1.00000   1.00000   1.00000   1.00000   1.00000

 Columns 92 through 100:

   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000  
1.00000   1.00000

and

flag = 0
relres =   4.4991e-014
iter =  50

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.

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/patch/?9282>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[patch #9282] add minres for sparse linear systems

Marco Caliari-5
Follow-up Comment #3, patch #9282 (project octave):

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

0 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.
0 I fixed some bugs, such as maxit and the using of preconditioner matrix.
0 I deleted the errors and call print_usage(). (I am not sure  whether I have
done it correctly.)
0 I added the description of usage.

(file #39970)
    _______________________________________________________

Additional Item Attachment:

File name: minres.m                       Size:21 KB


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/patch/?9282>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[patch #9282] add minres for sparse linear systems

Marco Caliari-5
Follow-up Comment #4, patch #9282 (project octave):

There are still a lot of problems now:

0 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.
0 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?
0 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.
0 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.





    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/patch/?9282>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[patch #9282] add minres for sparse linear systems

Marco Caliari-5
Follow-up Comment #5, patch #9282 (project octave):

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?

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/patch/?9282>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[patch #9282] add minres for sparse linear systems

Marco Caliari-5
Follow-up Comment #6, patch #9282 (project octave):

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.

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/patch/?9282>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[patch #9282] add minres for sparse linear systems

Marco Caliari-5
Follow-up Comment #7, patch #9282 (project octave):

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)
    _______________________________________________________

Additional Item Attachment:

File name: minres.m                       Size:20 KB


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/patch/?9282>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[patch #9282] add minres for sparse linear systems

Marco Caliari-5
Follow-up Comment #8, patch #9282 (project octave):

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.


>> A = A+A*1i;

>> x = minres(A,b,tol,maxit,M);
minres stopped at iteration 50 without converging to the desired tolerance
1e-10
because the maximum number of iterations was reached.
The iterate returned (number 2) has relative residual 0.7.


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/patch/?9282>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[patch #9282] add minres for sparse linear systems

Marco Caliari-5
Follow-up Comment #9, patch #9282 (project octave):

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.


n = 100; on = ones(n,1);
A = spdiags ([-2*on 4*on -2*on], -1:1, n, n);
b = sum(A,2);
tol = 1e-10;
maxit = 50;
M1 = spdiags(4*on,0,n,n);

A = spdiags([(-2 - 2 * 1i)*on (4 + 4)*on (-2 + 2 * 1i)*on],-1:1,n,n); %I added
this line

x = minres(A,b,tol,maxit,M1);

% [x,flag,relres,iter,resvec] = minres(A,b,tol,maxit,M1);
% use this if you want the full outputs.


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.

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/patch/?9282>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[patch #9282] add minres for sparse linear systems

Marco Caliari-5
Follow-up Comment #10, patch #9282 (project octave):

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),


N = 3;
A = rand(3);
A = A*A';
M = [1 0 0; 0 1 0; 0 0 0]; % the last rows is zero
[x,flag] = minres (A, ones(3,1), [], [], M);


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:


% Use many random hermitian A and complex b to test

N = 20; %Number of tests
n = 100; %Size of matrices
tol = 1e-10;
maxit = 2 * n;

flag_set = zeros(N, 1);
relres_set = zeros(N, 1);
iter_set = zeros(N, 1);

A_set = rand(N * n, n) + 1i * rand(N * n, n);
b_set = rand(N * n, 1) + 1i * rand(N * n, 1);

for j = 1: N
    A = A_set(((j - 1) * n + 1): j * n, :);
    A = (A + A') / 2; %Make A hermitian
    b = b_set(((j - 1) * n + 1): j * n, :);
   
    [x, flag, relres, iter, resvec] = minres (A, b, tol, maxit);
   
    flag_set(j) = flag;
    relres_set(j) = relres;
    iter_set(j) = iter;
end

allconverge = all(flag_set == zeros(N, 1));
max_converge_iter = max(iter_set);
min_converge_iter = min(iter_set);
max_relres = max(relres_set);


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




(file #40006)
    _______________________________________________________

Additional Item Attachment:

File name: minres.m                       Size:20 KB


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/patch/?9282>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[patch #9282] add minres for sparse linear systems

Marco Caliari-5
Follow-up Comment #11, patch #9282 (project octave):

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.


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/patch/?9282>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[patch #9282] add minres for sparse linear systems

Marco Caliari-5
Follow-up Comment #12, patch #9282 (project octave):

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)
    _______________________________________________________

Additional Item Attachment:

File name: minres.m                       Size:20 KB


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/patch/?9282>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[patch #9282] add minres for sparse linear systems

Marco Caliari-5
Follow-up Comment #13, patch #9282 (project octave):

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:


A = toeplitz (1:4);
M = toeplitz ([2,1,0,0]);
b = A * ones (4, 1);
M2 = chol (M);
M1 = M2';
minres (A ,b ,0 ,1, M1, M2)
minres (A, b, 0, 1, M) % it should be exactly the same
M2 \ minres (M1 \ A / M2, M1 \ b, 0, 1) % equivalent formulation


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


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/patch/?9282>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[patch #9282] add minres for sparse linear systems

Marco Caliari-5
Follow-up Comment #14, patch #9282 (project octave):

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

(file #40018)
    _______________________________________________________

Additional Item Attachment:

File name: minres.m                       Size:20 KB


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/patch/?9282>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/


Loading...