EricDarve / numerical_linear_algebra

Julia code for the book Numerical Linear Algebra
114 stars 41 forks source link

MINRES algorithm failes to converge and GMRES algorithm syntax error #198

Closed Zan-AA closed 4 years ago

Zan-AA commented 4 years ago

Hi Prof Darve,

I've been playing with the codes in your book for my own understanding and it seems that the MINRES algorithm on page 297 (which is a excerpt from the code available here) has a very poor convergence. I have tried to run the MINRES algorithm on general symmetric matrices and SPD matrices, none of the tests converge in n iterations (n being the size of the matrix). I have tried mostly on small size problems with no preconditioning. To debug, I compared your MINRES algorithm with the one provided from SOL lab (https://web.stanford.edu/group/SOL/software/minres/) and even the first iteration produces different parameter values for the solution update at the end of each iteration. Their version seems to be able to converge in n iteration for the same matrix. I suspect a typo in the algorithm somewhere but the conventions and symbols used in your version and SOL lab's version differ a lot so I couldn't really pin down the exact error. Do you mind double checking the source of the MINRES algorithm? in the meanwhile, I will try to investigate this during my own free time.

As for the GMRES algorithm, in gmres.jl line 53, the Julia version I am using (Julia 1.4.0 in jupyter notebook environment) gives me a syntax error:

"MethodError: Cannot convert an object of type Array{Float64,2} to an object of type Float64"

which bascially says the right hand side produces a 1-element array and it needs to be squeezed down to a scalar in order to be written into the H matrix. An indexed notation seems to solve the problem. H[k,i] = z' q[k] to H[k,i] = (z' q[k])[1]

Other than that, the GMRES algorithm works perfectly fine. For your reference, I have generated a simple code to test minres and gmres algorithm and the test results are shown here.

MINRES_GMRES (1).pdf

Best regards, Zan Xu

Zan-AA commented 4 years ago

Update: the MINRES algorithm works fine. The previously observed problem was attributed to the ill-conditioning of the test matrices. Well-conditioned matrix consistently converges within n iterations.

Please feel free to close this issue once the minor syntax problem in GMRES solver is resolved.

EricDarve commented 4 years ago

👍 I will wait for the GMRES problem to be fixed before closing.

leopoldcambier commented 4 years ago

@Zan-AA can you please provide an example for the gmres issue ? I can't reproduce that Method Error with Julia 1.4.1

EricDarve commented 4 years ago

Line 53 is

H[k,i] = z' * q[k]

with

z = M \ (A * q[i])

I am guessing that your choice of M is the problem. z is probably not a vector. You can check the type of z. z’ * q[k] should be a scalar. It’s a dot product between two vectors.

Eric

Zan-AA commented 4 years ago

Leopold and I discussed about this and it seems the problem lies in the type of vector z as the expression seems to work only if z is of type 1D array.

The reason could be due to precond M, but I also encountered the type problem when I declared RHS vector b as a 2D matrix with only 1 column. Regardless, this can be fixed by just calling the dot function in Julia, which covers the 2D matrix types.

H[k,i] = dot(z, q[k])

(This also explained why I didnt have this issue of using 2D matrix as vector in MINRES as MINRES algorithm uses Julia dot product function in place of explicit multiplications)

Zan-AA commented 4 years ago

I created a PR to fix the minor issue.

EricDarve commented 4 years ago

Great. Thanks.