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

**URL:** <https://discourse.julialang.org/t/unbounded-iterativesolvers-gmres-different-output-to-matlab-slows-convergence-significantly/21499>\
**Category:** Performance\
**Created:** [March 5, 2019, 4:48pm UTC](https://discourse.julialang.org/t/unbounded-iterativesolvers-gmres-different-output-to-matlab-slows-convergence-significantly/21499 "2019-03-05T16:48:13Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![timmm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/timmm/32/6760_2.png) [@timmm](https://discourse.julialang.org/u/timmm)\
**Post date:** [March 5, 2019, 4:48pm UTC](https://discourse.julialang.org/t/unbounded-iterativesolvers-gmres-different-output-to-matlab-slows-convergence-significantly/21499/1 "2019-03-05T16:48:13Z")

</div>

Hi all,

For project here: [GitHub - Timmmdavis/FischerNewton: Julia version of NUM4LCP.m](https://github.com/Timmmdavis/FischerNewton) 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](https://de.mathworks.com/help/matlab/ref/gmres.html)):

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

```

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

```julia
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:

```julia
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](https://github.com/Timmmdavis/CutAndDisplace/blob/master/SharedFunctions/FrictionSolver/Num4LCP/fischer_newton.m)

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

Cheers,

Tim

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [March 5, 2019, 7:11pm UTC](https://discourse.julialang.org/t/unbounded-iterativesolvers-gmres-different-output-to-matlab-slows-convergence-significantly/21499/2 "2019-03-05T19:11:21Z")

</div>

You’re likely to get more answers by posting at [Issues · JuliaLinearAlgebra/IterativeSolvers.jl · GitHub](https://github.com/JuliaMath/IterativeSolvers.jl/issues), and providing a simple reproducer. Did you check the x0 was the same as in matlab?

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [March 5, 2019, 10:58pm UTC](https://discourse.julialang.org/t/unbounded-iterativesolvers-gmres-different-output-to-matlab-slows-convergence-significantly/21499/3 "2019-03-05T22:58:14Z")

</div>

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 :

```julia
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:

```julia
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`,

```julia
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 :

```julia
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

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [March 5, 2019, 10:58pm UTC](https://discourse.julialang.org/t/unbounded-iterativesolvers-gmres-different-output-to-matlab-slows-convergence-significantly/21499/4 "2019-03-05T22:58:42Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [March 6, 2019, 6:22am UTC](https://discourse.julialang.org/t/unbounded-iterativesolvers-gmres-different-output-to-matlab-slows-convergence-significantly/21499/5 "2019-03-06T06:22:31Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [March 6, 2019, 7:06am UTC](https://discourse.julialang.org/t/unbounded-iterativesolvers-gmres-different-output-to-matlab-slows-convergence-significantly/21499/6 "2019-03-06T07:06:08Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![timmm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/timmm/32/6760_2.png) [@timmm](https://discourse.julialang.org/u/timmm)\
**Post date:** [March 6, 2019, 12:25pm UTC](https://discourse.julialang.org/t/unbounded-iterativesolvers-gmres-different-output-to-matlab-slows-convergence-significantly/21499/7 "2019-03-06T12:25:19Z")

</div>

Thanks for the responses,

@[antoine-levitt](https://discourse.julialang.org/u/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.

```julia

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(&quot;##34#37&quot;)){SparseArrays.SparseMatrixCSC{Float64,Int64},Array{Float64,1}},getfield(Krylov, Symbol(&quot;##35#38&quot;)){SparseArrays.SparseMatrixCSC{Float64,Int64},Array{Float64,1}},getfield(Krylov, Symbol(&quot;##36#39&quot;)){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(&quot;#kw##gmres&quot;)))(::NamedTuple{(:rtol, :x0, :itmax),Tuple{Float64,Array{Float64,1},Int64}}, ::typeof(Krylov.gmres), ::LinearOperators.LinearOperator{Complex{Float64},getfield(Krylov, Symbol(&quot;##34#37&quot;)){SparseArrays.SparseMatrixCSC{Float64,Int64},Array{Float64,1}},getfield(Krylov, Symbol(&quot;##35#38&quot;)){SparseArrays.SparseMatrixCSC{Float64,Int64},Array{Float64,1}},getfield(Krylov, Symbol(&quot;##36#39&quot;)){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(&quot;#kw##gmres&quot;)))(::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](https://discourse.julialang.org/u/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`:

```julia

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,&quot;\\test\\Matricies.mat&quot;)

else

G=string(G,&quot;/test/Matricies.mat&quot;)

end

println(G)

#Inputs - Loading influence matricies

using MAT

file = matopen(G)

A=read(file, &quot;A&quot;)

b=read(file, &quot;b&quot;)

#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(&quot;KrylovKit.GMRES&quot;)

for i = 1:10

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

end

println(&quot;Krylov.gmres&quot;)

for i = 1:10

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

end

println(&quot;Krylov.dqgmres&quot;)

for i = 1:10

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

end

```

Results:

```julia

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

```

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [March 6, 2019, 12:46pm UTC](https://discourse.julialang.org/t/unbounded-iterativesolvers-gmres-different-output-to-matlab-slows-convergence-significantly/21499/8 "2019-03-06T12:46:45Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![timmm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/timmm/32/6760_2.png) [@timmm](https://discourse.julialang.org/u/timmm)\
**Post date:** [March 6, 2019, 3:23pm UTC](https://discourse.julialang.org/t/unbounded-iterativesolvers-gmres-different-output-to-matlab-slows-convergence-significantly/21499/9 "2019-03-06T15:23:15Z")

</div>

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

```julia
@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

```

```julia
100.00000000000004
99.99999999999994
99.99999999999997
99.99999999999999
99.99999999999993
100.00000000000016
99.9999999999999
99.99999999999994
99.99999999999993
99.99999999999987

```

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [March 6, 2019, 3:57pm UTC](https://discourse.julialang.org/t/unbounded-iterativesolvers-gmres-different-output-to-matlab-slows-convergence-significantly/21499/10 "2019-03-06T15:57:31Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [March 7, 2019, 12:04am UTC](https://discourse.julialang.org/t/unbounded-iterativesolvers-gmres-different-output-to-matlab-slows-convergence-significantly/21499/11 "2019-03-07T00:04:34Z")

</div>

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 :

```julia
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

```

```julia
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

```

```julia
@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)

```

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [March 7, 2019, 5:58am UTC](https://discourse.julialang.org/t/unbounded-iterativesolvers-gmres-different-output-to-matlab-slows-convergence-significantly/21499/12 "2019-03-07T05:58:50Z")

</div>

> [@timmm](#):
>
> DQ

I am surprised how good Matlab is for allocations!

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [March 7, 2019, 8:17am UTC](https://discourse.julialang.org/t/unbounded-iterativesolvers-gmres-different-output-to-matlab-slows-convergence-significantly/21499/13 "2019-03-07T08:17:32Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [March 7, 2019, 9:06am UTC](https://discourse.julialang.org/t/unbounded-iterativesolvers-gmres-different-output-to-matlab-slows-convergence-significantly/21499/14 "2019-03-07T09:06:55Z")

</div>

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