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 function Priority:  5 - Normal
Status:  None Privacy:  Public
Assigned to:  None Open/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)

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:

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.

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

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)

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.

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.

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.

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

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.

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.

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

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

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 Files:
   
   
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
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by xierui (Submitted the item)
  • -email is unavailable- added by xierui ([GSoC 2017 Application] add minres function for sparse linear systems)
  • -email is 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 add your encouragement to it.
    This task has 0 encouragements so far.

    Only logged-in users can vote.

     

     

     

    Follow 9 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-03-17 xierui Attached File- => Added minres.m, #40018
    2017-03-16 xierui Attached File- => Added minres.m, #40010
    2017-03-15 xierui Attached File- => Added minres.m, #40006
    2017-03-13 xierui Attached File- => Added minres.m, #39988
    2017-03-11 xierui Attached File- => Added minres.m, #39970
    2017-03-08 xierui Attached File- => Added minres.m, #39940
        Attached File- => Added minres.patch, #39941
        Carbon-Copy- => Added marco.caliari@univr.it
        Carbon-Copy- => Added octave-maintainers@gnu.org

    Back to the top


    Powered by Savane 3.5