Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix GMRES memory usage #131

Merged

Conversation

@haampie
Copy link
Member

@haampie haampie commented Jun 28, 2017

This PR contains a replacement for the current GMRES implementation aimed to fix #127. Feedback is welcome :)

What's different w.r.t. the current implementation:

  1. Pre-allocates (almost) all vectors. The only problem is how to deal with preconditioners efficiently
  2. Computes the residual without computing the QR decomposition of the Hessenberg matrix. This way you don't have to store the Given's rotations + the inner loop of GMRES gets cleaner
  3. Adds a special Hessenberg type (different from LinAlg.Hessenberg) that takes care of the small least-squares problem via Given's rotations + UpperTriangular.
  4. Uses UnsafeVectorView for non-allocating, small views of vectors. Bug prone, so removed.
  5. Only append to the ConvergenceHistory when logging is turned on
  6. Adds an Identity type as preconditioner so that without preconditioning we can do a non-allocating expansion of the Krylov subspace (i.e. V[k + 1] = A * V[k] will not allocate`)
  7. Adds some benchmarks
  8. Changes V to a Matrix rather than Vector{Vector} so that BLAS-2 can be exploited in the update of the solution x += Vy.

Todo:

  • Rename to gmres
  • Add different orthogonalization strategies~

Running LinearSystemsBench.gmres() gives something like:

julia> improved
BenchmarkTools.Trial:
  memory estimate:  13.75 MiB
  allocs estimate:  313
  --------------
  minimum time:     11.454 s (0.00% GC)

julia> old
BenchmarkTools.Trial:
  memory estimate:  1.26 GiB
  allocs estimate:  13378
  --------------
  minimum time:     13.568 s (8.59% GC)
Copy link
Member

@andreasnoack andreasnoack left a comment

Great work


@inline Base.size(H::Hessenberg, args...) = size(H.H, args...)

function solve!(H::Hessenberg, rhs)
Copy link
Member

@andreasnoack andreasnoack Jun 28, 2017

We usually call these A_ldiv_B!

Copy link
Member Author

@haampie haampie Jun 29, 2017

I wasn't completely sure, because this method mutates the H matrix as well

Copy link
Member

@andreasnoack andreasnoack Jul 6, 2017

I think we should stick with A_ldiv_B! here even though it invalidates the factorization. Probably a good idea to add a comment, though.


# Solve the upper triangular problem.
U = UpperTriangular(view(H.H, 1 : width, 1 : width))
U \ UnsafeVectorView(rhs, 1 : width)
Copy link
Member

@andreasnoack andreasnoack Jun 28, 2017

Why not in-place solve here?

Copy link
Member Author

@haampie haampie Jun 29, 2017

I'm doing that now in the latest commit, but it's a bit awkward because it's a least-squares problem. The Hessenberg matrix is not square, so it comes down to having the solution in the first k elements of rhs and having a last element at index k + 1 not being part of the solution.

@@ -0,0 +1,39 @@
import Base.LinAlg: Givens, givensAlgorithm

mutable struct Hessenberg{T<:AbstractMatrix}
Copy link
Member

@andreasnoack andreasnoack Jun 28, 2017

We should consider upstreaming this to Base

r.β = β
end

function solve_least_squares!(arnoldi::ArnoldiDecomp{T}, β, k::Int) where {T}
Copy link
Member

@andreasnoack andreasnoack Jun 28, 2017

Maybe a more precise name here since it is not solving arbitrary least squares problems

@codecov-io
Copy link

@codecov-io codecov-io commented Jun 29, 2017

Codecov Report

Merging #131 into master will decrease coverage by 1.39%.
The diff coverage is 78.68%.

Impacted file tree graph

@@           Coverage Diff            @@
##           master    #131     +/-   ##
========================================
- Coverage   86.29%   84.9%   -1.4%     
========================================
  Files          18      19      +1     
  Lines        1394    1444     +50     
========================================
+ Hits         1203    1226     +23     
- Misses        191     218     +27
Impacted Files Coverage Δ
src/hessenberg.jl 100% <100%> (ø)
src/common.jl 82.6% <100%> (+3.66%) ⬆️
src/orthogonalize.jl 25% <25%> (ø)
src/gmres.jl 97.36% <97.14%> (+6.45%) ⬆️
src/krylov.jl 50.98% <0%> (-15.69%) ⬇️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 4e92e2d...24e5bbd. Read the comment docs.

@haampie haampie force-pushed the fix-gmres-memory-usage branch 3 times, most recently from 15dd13c to c17265f Jul 2, 2017
@haampie haampie force-pushed the fix-gmres-memory-usage branch from c17265f to a86d093 Jul 6, 2017
@andreasnoack
Copy link
Member

@andreasnoack andreasnoack commented Jul 6, 2017

I think we a ready for this one so if you finish it up we can go ahead and merge

@haampie haampie force-pushed the fix-gmres-memory-usage branch from a86d093 to d97462a Jul 6, 2017
@haampie haampie force-pushed the fix-gmres-memory-usage branch from d92d69a to 24e5bbd Jul 6, 2017
@andreasnoack andreasnoack merged commit 81e7872 into JuliaLinearAlgebra:master Jul 6, 2017
1 of 2 checks passed
haampie added a commit to haampie/IterativeSolvers.jl that referenced this issue Jul 10, 2017
* Add improved versions of GMRES and CG

* Add preconditioners

* Add Hessenberg type

* Hessenberg solver

* Fix incorrect comment

* Benchmark Givens rotation + UpperTriangular vs \

* Rename to Hessenberg

* Use UnsafeVectorView when solving the small system

* Clean up some

* Add Identity preconditioner type

* Improve GMRES with specialized functions for preconditioners, unsafe views, and conditional logging

* Don't run by default

* Some default values for gmres bench

* Reduce the number of allocs slightly by solving the least-squares-problem with the Hessenberg matrix in-place

* Fix for complex numbers

* Move from Vector{Vector{T}} to Matrix{T} using views and make orthogonalization routine reusable

* Make the comparison fair and the matrix sparser to ~10 nonzeros per row

* Fix a bug for complex numbers

* Temporarily remove tests with matrices as preconditioner, since we don't have A_ldiv_B for them at the moment

* Add tests for various orthogonalization methods

* Fix a problem for complex numbers again

* Export Identity

* Add tests to runtests.jl

* Fix compatibility with Julia 0.5

* Override the current gmres method and move solve! to A_ldiv_B!
haampie added a commit to haampie/IterativeSolvers.jl that referenced this issue Jul 11, 2017
* Add improved versions of GMRES and CG

* Add preconditioners

* Add Hessenberg type

* Hessenberg solver

* Fix incorrect comment

* Benchmark Givens rotation + UpperTriangular vs \

* Rename to Hessenberg

* Use UnsafeVectorView when solving the small system

* Clean up some

* Add Identity preconditioner type

* Improve GMRES with specialized functions for preconditioners, unsafe views, and conditional logging

* Don't run by default

* Some default values for gmres bench

* Reduce the number of allocs slightly by solving the least-squares-problem with the Hessenberg matrix in-place

* Fix for complex numbers

* Move from Vector{Vector{T}} to Matrix{T} using views and make orthogonalization routine reusable

* Make the comparison fair and the matrix sparser to ~10 nonzeros per row

* Fix a bug for complex numbers

* Temporarily remove tests with matrices as preconditioner, since we don't have A_ldiv_B for them at the moment

* Add tests for various orthogonalization methods

* Fix a problem for complex numbers again

* Export Identity

* Add tests to runtests.jl

* Fix compatibility with Julia 0.5

* Override the current gmres method and move solve! to A_ldiv_B!
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Linked issues

Successfully merging this pull request may close these issues.

3 participants