Unbounded IterativeSolvers.Gmres! different output to MATLAB, slows convergence significantly

General Usage Performance
1

Hi all,

For project here: GitHub - Timmmdavis/FischerNewton: Julia version of NUM4LCP.m I am having issues with the results of IterativeSolvers.Gmres! The main code is under: FischerNewton/src/fischer_newton.jl

If I add the same inputs and settings as MATLABS GMRES solver, I cannot get exactly the same output (I am probably doing something wrong). The current output is different to such an extent it seriously effects the convergence criteria of the rest of my code.

Compare MATLAB vs Julia to the matrices contained in FischerNewton/test/Matricies.mat

Using MATLAB’s GMRES (Solve system of linear equations — generalized minimum residual method - MATLAB gmres - MathWorks Deutschland):

[dx] = gmres( A, b,10, 1e-6,10);

And Julia’s (simple code to load MATs in Julia @ FischerNewton/test/test_Compare2MATLAB.jl)

dx= IterativeSolvers.gmres( A, b ,tol=1e-6,restart=10, maxiter=10);

The results are different. The first 10 julia results as a percentage of the MATLAB result:

95.5277

92.8080

92.1084

95.6201

74.3765

62.9217

69.1012

77.8773

87.8614

92.2786

Are there obvious flags in Julia I can set to get these to match (or converge to a closer solution?)

I am trying to produce a faster example of my MATLAB LCP code here (that calls GMRES): CutAndDisplace/fischer_newton.m at master · Timmmdavis/CutAndDisplace · GitHub

At the moment due to problem described above this converges well over 10* faster than Julia on the large matrix.

Cheers,

Tim

2

You’re likely to get more answers by posting at Issues · JuliaLinearAlgebra/IterativeSolvers.jl · GitHub, and providing a simple reproducer. Did you check the x0 was the same as in matlab?

1 Like
3

Hi Tim,

I think that the problem comes from the implementation of restart in IterativeSolvers.gmres.
Or like you, I don’t understand how works their keyword argument restart.

I tested your code :

using MATLAB
using FischerNewton
import IterativeSolvers

include("test_Compare2MATLAB.jl")
b = b[:,1]

x_matlab = mat"gmres($A, $b, 10, 1e-6, 10)"
x_itsolver = IterativeSolvers.gmres( A, b, tol=1e-6, restart=10, maxiter=10)

And I have the same resultats:

for i = 1:10
       println(x_itsolver[i] / x_matlab[i] * 100)
end
95.52770715094717
92.80799500742638
92.10843453121272
95.62006975408472
74.37649780267112
62.921739325120015
69.10121867573282
77.87726417413496
87.86140357012937
92.27859017728773

I tried with an other implementation of GMRES,

import Krylov
n, m = size(A)
global x_krylov
x_krylov = zeros(n)
nb_restart = 10
for i = 1 : nb_restart
 global x_krylov
  (x_krylov, stats) = Krylov.gmres(A, b, rtol=1e-6, x0=x_krylov, itmax=10)
end

And this time we have the same solution :

for i = 1:10
       global x_krylov
       println(x_krylov[i] / x_matlab[i] * 100)
end
100.00008918762333
100.00015135484914
100.00015326343438
100.00010172365612
100.0003542622564
100.0004725828038
100.00040060447826
100.00031282868942
100.00020481906455
100.00015158974698

I don’t know if IterativeSolvers do the same thing. When x₀ ≠ 0 (restarted version), GMRES implementation in Krylov resolves A * Δx = r₀ with r₀ = b - Ax₀ and returns x = Δx + x₀ .

Cheers,
Alexis

1 Like
4

For the moment GMRES in not in Krylov.jl repository. We are working on adding it.
If you want to test it : julia > ] add https://github.com/amontoison/Krylov.jl.git#gmres

For your problem, did you test with DQGMRES algorithm instead of restarting GMRES? It could be faster too.

5

I don’t have Matlab so I cant test. What about trying KrylovKit?

6

Could it be the results are just shifted by one, so that the first iteration in Matlab corresponds to the second in IterativeSolvers?

7

Thanks for the responses,

@antoine-levitt I don’t think this is the case, otherwise my code would converge as fast as MATLAB did.

Thanks @amontoison, it is good to see a gmres solver that almost reproduces the results of MATLABS solver, having tried this it occasionally throws the following error.


ERROR: LoadError: BoundsError: attempt to access 1-element Array{Float64,1} at index [0]

Stacktrace:

 [1] getindex at .\array.jl:729 [inlined]

 [2] #gmres#15(::Array{Float64,1}, ::LinearOperators.opEye, ::LinearOperators.opEye, ::Float64, ::Float64, ::Int64, ::Bool, ::Bool, ::typeof(Krylov.gmres), ::LinearOperators.LinearOperator{Complex{Float64},getfield(Krylov, Symbol("##34#37")){SparseArrays.SparseMatrixCSC{Float64,Int64},Array{Float64,1}},getfield(Krylov, Symbol("##35#38")){SparseArrays.SparseMatrixCSC{Float64,Int64},Array{Float64,1}},getfield(Krylov, Symbol("##36#39")){SparseArrays.SparseMatrixCSC{Float64,Int64},Array{Float64,1}}}, ::Array{Float64,1}) at C:\Users\timmm\.julia\packages\Krylov\ybdMT\src\gmres.jl:139

 [3] (::getfield(Krylov, Symbol("#kw##gmres")))(::NamedTuple{(:rtol, :x0, :itmax),Tuple{Float64,Array{Float64,1},Int64}}, ::typeof(Krylov.gmres), ::LinearOperators.LinearOperator{Complex{Float64},getfield(Krylov, Symbol("##34#37")){SparseArrays.SparseMatrixCSC{Float64,Int64},Array{Float64,1}},getfield(Krylov, Symbol("##35#38")){SparseArrays.SparseMatrixCSC{Float64,Int64},Array{Float64,1}},getfield(Krylov, Symbol("##36#39")){SparseArrays.SparseMatrixCSC{Float64,Int64},Array{Float64,1}}}, ::Array{Float64,1}) at .\none:0

 [4] #gmres#50(::Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol,Symbol},NamedTuple{(:rtol, :x0, :itmax),Tuple{Float64,Array{Float64,1},Int64}}}, ::Function, ::SparseArrays.SparseMatrixCSC{Float64,Int64}, ::Array{Float64,1}) at C:\Users\timmm\.julia\packages\Krylov\ybdMT\src\variants.jl:15

 [5] (::getfield(Krylov, Symbol("#kw##gmres")))(::NamedTuple{(:rtol, :x0, :itmax),Tuple{Float64,Array{Float64,1},Int64}}, ::typeof(Krylov.gmres), ::SparseArrays.SparseMatrixCSC{Float64,Int64}, ::Array{Float64,1}) at .\none:0

 [6] fischer_newton(::Array{Float64,2}, ::Array{Float64,2}) at C:\Users\timmm\.julia\packages\FischerNewton\myRMj\src\fischer_newton.jl:190

I also tried your DQGMRES and this works well (without errors). I am using this for the moment.

@rveltz I have tried out the KrylovKit.GMRES also and it gives the same (at least to 8 significant figures) error compared to MATLAB on the matrix I was using as Alexis showed for the Krylov.gmres implementation. This also works well.

Testscript comparing KrylovKit and Krylov:


using FischerNewton

using IterativeSolvers

using KrylovKit

using Krylov

using LinearAlgebra

using SparseArrays

using Test

using Profile

using InteractiveUtils

G=pathof(FischerNewton);

G=splitdir(G); #remove file name

G=G[1];

G=splitdir(G); #out of src

G=G[1];

if Sys.iswindows()

G=string(G,"\\test\\Matricies.mat")

else

G=string(G,"/test/Matricies.mat")

end

println(G)

#Inputs - Loading influence matricies

using MAT

file = matopen(G)

A=read(file, "A")

b=read(file, "b")

#Reshape to N length array (not two dims)

b=reshape(b,length(b));

x0=zeros(length(b))

restart=10;

alg = GMRES( krylovdim = restart, maxiter = 10, tol = 1e-6)

@time x_itsolver1, stats1 = @inferred linsolve(A,b,x0,alg)

@info stats1

r1 = norm(b - A * x_itsolver1)/norm(b)

@info r1

x_itsolver2=copy(x0);

@time for i = 1 : restart

(x2,stats2) = Krylov.gmres(A, b, rtol=1e-6, x0=x0, itmax=10)

if i==restart

@info stats2

x_itsolver2.=x2;

end

end

r2 = norm(b - A * x_itsolver2)/norm(b)

@info r2

dqgmres_tol = 1.0e-6

@time (x_itsolver3,stats3) = Krylov.dqgmres(A, b,atol=dqgmres_tol,itmax=restart)

@info stats3

r3 = norm(b - A * x_itsolver3)/norm(b)

@info r3

x_matlab=

[43189395.5369333;

65338300.3629727;

52631061.5842138;

54775644.5591780;

36174129.5940650;

32348224.0211539;

39287541.3646495;

54096387.5483277;

78831593.1510365;

70208640.7061015]

println("KrylovKit.GMRES")

for i = 1:10

println(x_itsolver1[i] / x_matlab[i] * 100)

end

println("Krylov.gmres")

for i = 1:10

println(x_itsolver2[i] / x_matlab[i] * 100)

end

println("Krylov.dqgmres")

for i = 1:10

println(x_itsolver3[i] / x_matlab[i] * 100)

end

Results:


RelativeResidual(res) being: norm(b-A*x)/norm(b)

MATLAB:

[dx,~,res] = gmres( A, b,10, 1e-6,10);

Elapsed time is 0.156341 seconds.

res=2.7714e-06

KrylovKit.GMRES:

0.262323 seconds (304 allocations: 2.472 MiB)

res=8.920428143716319e-9

#Percentage of MATLAB result

100.00008918762342

100.00015135484918

100.00015326343437

100.00010172365623

100.00035426225628

100.00047258280402

100.00040060447817

100.00031282868942

100.0002048190644

100.0001515897469

Krylov.gmres

0.307856 seconds (12.16 k allocations: 1.670 MiB)

res=8.920428259577072e-9

#Percentage of MATLAB result

100.00008918762344

100.00015135484917

100.00015326343438

100.00010172365617

100.00035426225634

100.00047258280395

100.0004006044781

100.00031282868933

100.00020481906455

100.00015158974686

Krylov.dqgmres

0.029348 seconds (56 allocations: 285.250 KiB)

res=0.061037625187066416

#Percentage of MATLAB result

95.52770715094722

92.80799500742638

92.10843453121275

95.62006975408465

74.37649780267117

62.92173932512013

69.1012186757328

77.87726417413496

87.86140357012935

92.27859017728767
8

Matlab has an interface where maxit*restart total iterations (= matvecs) are performed, whereas the maxiter option in IterativeSolvers is directly the total number of iterations (the latter is much saner IMHO). That may be what’s tripping you up. If not, please open an issue over at IterativeSolvers.

9

You are correct, if I do this as in the code in my last post above I get the result:
IterativeSolvers.gmres

@time x_itsolver4=IterativeSolvers.gmres(A, b,tol=1e-6,restart=restart, initially_zero=true,maxiter=10*restart);
r4 = norm(b - A * x_itsolver4)/norm(b)
@info r4

println("IterativeSolvers.gmres")
for i = 1:10
       println(x_itsolver4[i] / x_matlab[i] * 100)
end

100.00000000000004
99.99999999999994
99.99999999999997
99.99999999999999
99.99999999999993
100.00000000000016
99.9999999999999
99.99999999999994
99.99999999999993
99.99999999999987
10

Likewise, I believe you want to run dqgmres with options memory=restart and maxiter=10*restart.

1 Like
11

Thanks @antoine-levitt ! It’s difficult to understand how works keywords arguments sometimes.
@dpo is right about Krylov.dqgmres option.

A little benchmark to see which one is the fastest :

using MATLAB
using FischerNewton
using BenchmarkTools
using LinearOperators

import KrylovKit
import IterativeSolvers
import Krylov

include("test_Compare2MATLAB.jl")
op = PreallocatedLinearOperator(A)
b = b[:,1]

nb_restart = 10
itmax      = 10 # Maximal number of iterations for one pass in restarted GMRES

x_matlab          = mat"gmres($A, $b, 10, 1e-6, 10)"
x_itsolver        = IterativeSolvers.gmres(A, b, tol=1e-6, restart=nb_restart, maxiter=itmax*nb_restart)
(x_krylov, stats) = Krylov.dqgmres(op, b, rtol=1e-6, memory=itmax, itmax=itmax*nb_restart)

println("IterativeSolvers.gmres")
for i = 1:10
	println(x_itsolver[i] / x_matlab[i] * 100)
end
println()
println("Krylov.dqgmres")
for i = 1:10
	println(x_krylov[i] / x_matlab[i] * 100)
end

IterativeSolvers.gmres
100.00000000000009
100.00000000000004
100.00000000000007
99.99999999999993
100.00000000000004
99.99999999999997
100.00000000000003
100.0
99.99999999999997
100.00000000000004

Krylov.dqgmres
99.99996697271747
99.99996297022264
99.99995815977995
99.99996796962814
99.99990375571936
99.99983119752942
99.9998527740654
99.99988897069454
99.9999398964075
99.99994235348828

@btime mat"gmres($A, $b, 10, 1e-6, 10)"
@btime IterativeSolvers.gmres($A, $b, tol=1e-6, restart=nb_restart, maxiter=itmax*nb_restart)
@btime Krylov.dqgmres($op, $b, rtol=1e-6, memory=itmax, itmax=itmax*nb_restart)

110.246 ms (53 allocations: 14.42 KiB)
85.096 ms (992 allocations: 210.73 KiB)
78.740 ms (49 allocations: 262.34 KiB)
1 Like
12

I am surprised how good Matlab is for allocations!

13

I think, Julia is not able to measure memory allocated in MATLAB functions.
If we do 10 iterations for each inner GMRES, we need to allocate at least 10 vectors of dimension n for bulding the orthogonal basis V during the Arnoldi process. n=1560 with @timmm problem. Each coefficient is stored with 8 bytes (64 bits).

So we need (at least) 1560 * 8 * 10 / 1024 = 121.875 KiB.

The solution x of the problem is a n-vector : 12.1875 KiB, so I suppose that Julia measure only that when the solution is returned.

1 Like
14

The speed is not too bad as well… I am guessing Matlab is multithreaded