# [patch #9282] add minres for sparse linear systems

15 messages
Open this post in threaded view
|

## [patch #9282] add minres for sparse linear systems

 URL:                    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 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 ------------------------------------------------------- Date: Wed 08 Mar 2017 01:13:15 PM UTC  Name: minres.patch  Size: 11kB   By: xierui add minres for sparse linear systems     _______________________________________________________ Reply to this item at:   _______________________________________________   Message sent via/by Savannah   http://savannah.gnu.org/
Open this post in threaded view
|

## [patch #9282] add minres for sparse linear systems

 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:   _______________________________________________   Message sent via/by Savannah   http://savannah.gnu.org/
Open this post in threaded view
|

## [patch #9282] add minres for sparse linear systems

 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:   _______________________________________________   Message sent via/by Savannah   http://savannah.gnu.org/
Open this post in threaded view
|

## [patch #9282] add minres for sparse linear systems

 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:   _______________________________________________   Message sent via/by Savannah   http://savannah.gnu.org/
Open this post in threaded view
|

## [patch #9282] add minres for sparse linear systems

 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:   _______________________________________________   Message sent via/by Savannah   http://savannah.gnu.org/
Open this post in threaded view
|

## [patch #9282] add minres for sparse linear systems

 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:   _______________________________________________   Message sent via/by Savannah   http://savannah.gnu.org/
Open this post in threaded view
|

## [patch #9282] add minres for sparse linear systems

 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:   _______________________________________________   Message sent via/by Savannah   http://savannah.gnu.org/
Open this post in threaded view
|

## [patch #9282] add minres for sparse linear systems

 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:   _______________________________________________   Message sent via/by Savannah   http://savannah.gnu.org/
Open this post in threaded view
|

## [patch #9282] add minres for sparse linear systems

 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:   _______________________________________________   Message sent via/by Savannah   http://savannah.gnu.org/
Open this post in threaded view
|

## [patch #9282] add minres for sparse linear systems

 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:   _______________________________________________   Message sent via/by Savannah   http://savannah.gnu.org/
Open this post in threaded view
|

## [patch #9282] add minres for sparse linear systems

 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:   _______________________________________________   Message sent via/by Savannah   http://savannah.gnu.org/
Open this post in threaded view
|

## [patch #9282] add minres for sparse linear systems

 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:   _______________________________________________   Message sent via/by Savannah   http://savannah.gnu.org/
Open this post in threaded view
|

## [patch #9282] add minres for sparse linear systems

 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:   _______________________________________________   Message sent via/by Savannah   http://savannah.gnu.org/