diff --git a/src/BitStringAddresses/occupationnumberfs.jl b/src/BitStringAddresses/occupationnumberfs.jl
index 6cad03a3a..5383c3300 100644
--- a/src/BitStringAddresses/occupationnumberfs.jl
+++ b/src/BitStringAddresses/occupationnumberfs.jl
@@ -66,7 +66,14 @@ end
 Base.reverse(ofs::OccupationNumberFS) = OccupationNumberFS(reverse(ofs.onr))
 
 onr(ofs::OccupationNumberFS) = ofs.onr
-Base.isless(a::OccupationNumberFS, b::OccupationNumberFS) = isless(b.onr, a.onr)
+function Base.isless(a::OccupationNumberFS{M}, b::OccupationNumberFS{M}) where {M}
+    # equivalent to `isless(reverse(a.onr), reverse(b.onr))`
+    i = length(a.onr)
+    while i > 1 && a.onr[i] == b.onr[i]
+        i -= 1
+    end
+    return isless(a.onr[i], b.onr[i])
+end
 # reversing the order here to make it consistent with BoseFS
 Base.:(==)(a::OccupationNumberFS, b::OccupationNumberFS) = a.onr == b.onr
 Base.hash(ofs::OccupationNumberFS, h::UInt) = hash(ofs.onr, h)
diff --git a/src/ExactDiagonalization/ExactDiagonalization.jl b/src/ExactDiagonalization/ExactDiagonalization.jl
index e54f9681d..805b4121b 100644
--- a/src/ExactDiagonalization/ExactDiagonalization.jl
+++ b/src/ExactDiagonalization/ExactDiagonalization.jl
@@ -25,14 +25,16 @@ using CommonSolve: CommonSolve, solve, init
 using VectorInterface: VectorInterface, add
 using OrderedCollections: freeze
 using NamedTupleTools: delete
+using StaticArrays: setindex
 
 using Rimu: Rimu, DictVectors, Hamiltonians, Interfaces, BitStringAddresses, replace_keys,
     clean_and_warn_if_others_present
 using ..Interfaces: AbstractDVec, AbstractHamiltonian, AdjointUnknown,
     diagonal_element, offdiagonals, starting_address, LOStructure, IsHermitian
-using ..BitStringAddresses: AbstractFockAddress
+using ..BitStringAddresses: AbstractFockAddress, BoseFS, FermiFS, CompositeFS, near_uniform
 using ..DictVectors: FrozenDVec, PDVec, DVec
-using ..Hamiltonians: check_address_type, dimension, ParitySymmetry, TimeReversalSymmetry
+using ..Hamiltonians: allows_address_type, check_address_type, dimension,
+    ParitySymmetry, TimeReversalSymmetry
 
 
 export ExactDiagonalizationProblem, KrylovKitSolver, LinearAlgebraSolver
@@ -42,6 +44,8 @@ export BasisSetRepresentation, build_basis
 export sparse # from SparseArrays
 
 
+include("basis_breadth_first_search.jl")
+include("basis_fock.jl")
 include("basis_set_representation.jl")
 include("algorithms.jl")
 include("exact_diagonalization_problem.jl")
diff --git a/src/ExactDiagonalization/basis_breadth_first_search.jl b/src/ExactDiagonalization/basis_breadth_first_search.jl
new file mode 100644
index 000000000..952706183
--- /dev/null
+++ b/src/ExactDiagonalization/basis_breadth_first_search.jl
@@ -0,0 +1,399 @@
+const SubVector{T} = SubArray{T,1,Vector{T},Tuple{UnitRange{Int64}},true}
+"""
+    UniformSplit(array::Vector, min_chunk_size, max_chunks)
+
+Split the array into at most `max_chunks` subarrays, with each at least `min_chunk_size` long.
+If the number of elements in `array` is not divisible by the determined chunk size, the leftover elements are placed in the last chunk.
+
+```jldoctest
+julia> using Rimu.ExactDiagonalization: UniformSplit
+
+julia> UniformSplit(collect(1:10), 3, 3)
+3-element UniformSplit{Int64}:
+ [1, 2, 3]
+ [4, 5, 6]
+ [7, 8, 9, 10]
+
+julia> UniformSplit(collect(1:13), 5, 3)
+2-element UniformSplit{Int64}:
+ [1, 2, 3, 4, 5, 6]
+ [7, 8, 9, 10, 11, 12, 13]
+```
+"""
+struct UniformSplit{T} <: AbstractVector{SubVector{T}}
+    array::Vector{T}
+    n_chunks::Int
+    chunk_size::Int
+
+    function UniformSplit(array::Vector{T}, min_chunk_size, max_chunks) where {T}
+        chunk_size = fld(length(array), max_chunks)
+        if min_chunk_size > chunk_size
+            n_chunks = max(fld(length(array), min_chunk_size), 1)
+            chunk_size = fld(length(array), n_chunks)
+        else
+            n_chunks = max(fld(length(array), chunk_size), 1)
+        end
+        return new{T}(array, n_chunks, chunk_size)
+    end
+end
+
+Base.size(split::UniformSplit) = (split.n_chunks,)
+
+function Base.getindex(split::UniformSplit, i)
+    index_start = (i - 1) * split.chunk_size + 1
+    if i == split.n_chunks
+        index_end = length(split.array)
+    elseif 0 < i < split.n_chunks
+        index_end = i * split.chunk_size
+    else
+        throw(BoundsError(split, i))
+    end
+    return view(split.array, index_start:index_end)
+end
+
+"""
+    bb! = BasisBuilder{A::Type}(; col_hint=0)
+    bb!(operator::AbstractOperator, pairs, seen)
+
+Functor used with [`basis_breadth_first_search`](@ref) to build a basis of addresses of type
+`A` from an operator. It contains a set of addresses (`bb!.frontier`) that is collected from
+the `offdiagonals` of `operator` with all addresses contained the list of address-value
+pairs `pairs` that are not element of `seen`.
+"""
+struct BasisBuilder{A}
+    frontier::Set{A}
+
+    function BasisBuilder{A}(; col_hint=0) where {A}
+        return new{A}(sizehint!(Set{A}(), col_hint))
+    end
+end
+init_accumulator(::Type{<:BasisBuilder}) = nothing
+update_accumulator!(_, ::BasisBuilder, _) = nothing
+function finalize_accumulator!(::Nothing, basis, sort; kwargs...)
+    if sort
+        return sort!(basis; kwargs...)
+    else
+        return basis
+    end
+end
+
+function (bb::BasisBuilder)(operator, ks, seen)
+    empty!(bb.frontier)
+    for (k1, _) in ks
+        for (k2, val) in offdiagonals(operator, k1)
+            if val ≠ 0 && k2 ∉ seen
+                push!(bb.frontier, k2)
+            end
+        end
+    end
+end
+
+"""
+    mb! = MatrixBuilder{A::Type}(; col_hint=0)
+    mb!(operator::AbstractOperator, pairs, seen)
+
+Functor used with [`basis_breadth_first_search`](@ref) to build a matrix from an
+operator. It contains a set of addresses (`mb!.frontier`) that is collected from the
+`offdiagonals` of `operator` with all addresses contained the list of address-value pairs
+`pairs` that are not element of `seen`.
+
+It also collects `mb!.is`, `mb!.js`, and `mb!.vs` which are used to build the sparse matrix
+via `sparse!`.
+"""
+struct MatrixBuilder{A,T,I}
+    js::Vector{I}
+    is::Vector{A}
+    vs::Vector{T}
+    column_buffer::Dict{A,T}
+    frontier::Set{A}
+
+    function MatrixBuilder{A,T,I}(; col_hint=0) where {A,T,I}
+        return new{A,T,I}(
+            sizehint!(I[], col_hint),
+            sizehint!(A[], col_hint),
+            sizehint!(T[], col_hint),
+            sizehint!(Dict{A,T}(), col_hint),
+            sizehint!(Set{A}(), col_hint),
+        )
+    end
+end
+
+function init_accumulator(::Type{<:MatrixBuilder{<:Any,T,I}}) where {I,T}
+    return MatrixBuilderAccumulator{T,I}()
+end
+
+function (builder::MatrixBuilder{<:Any,T})(operator, columns, seen) where {T}
+    empty!(builder.is)
+    empty!(builder.js)
+    empty!(builder.vs)
+    empty!(builder.frontier)
+    for (col_key, col_index) in columns
+        empty!(builder.column_buffer)
+        diag = diagonal_element(operator, col_key)
+        if !iszero(diag)
+            builder.column_buffer[col_key] = diag
+        end
+        for (row_key, value) in offdiagonals(operator, col_key)
+            if !iszero(value)
+                new_value = get(builder.column_buffer, row_key, zero(T)) + value
+                builder.column_buffer[row_key] = new_value
+                row_key ∉ seen && push!(builder.frontier, row_key)
+            end
+        end
+
+        old_len = length(builder.is)
+        new_len = length(builder.is) + length(builder.column_buffer)
+        resize!(builder.is, new_len)
+        resize!(builder.js, new_len)
+        resize!(builder.vs, new_len)
+        @inbounds for (i, (row_key, value)) in enumerate(builder.column_buffer)
+            builder.is[old_len + i] = row_key
+            builder.js[old_len + i] = col_index
+            builder.vs[old_len + i] = value
+        end
+    end
+end
+
+"""
+    struct MatrixBuilderAccumulator
+
+Used in conjunction with [`basis_breadth_first_search`](@ref) and
+[`MatrixBuilder`](@ref). It is used to combine the sparse matrix as it is being built by
+multiple threads.
+"""
+struct MatrixBuilderAccumulator{T,I}
+    is::Vector{I}
+    js::Vector{I}
+    vs::Vector{T}
+
+    MatrixBuilderAccumulator{T,I}() where {T,I} = new{T,I}(I[], I[], T[])
+end
+
+function update_accumulator!(acc::MatrixBuilderAccumulator, mb::MatrixBuilder, mapping)
+    for i in eachindex(mb.is)
+        if haskey(mapping, mb.is[i]) # skip keys that were filtered
+            push!(acc.is, mapping[mb.is[i]])
+            push!(acc.js, mb.js[i])
+            push!(acc.vs, mb.vs[i])
+        end
+    end
+end
+
+function finalize_accumulator!(
+    acc::MatrixBuilderAccumulator{T,I}, basis, sort; kwargs...
+) where {T,I}
+    n = length(basis)
+
+    # see docstring of `sparse!` for an explanation for what this is
+    klasttouch = Vector{I}(undef, n)
+    csrrowptr = Vector{I}(undef, n + 1)
+    csrcolval = Vector{I}(undef, length(acc.is))
+    csrnzval = Vector{T}(undef, length(acc.is))
+
+    matrix = SparseArrays.sparse!(
+        acc.is, acc.js, acc.vs, n, n, +,
+        klasttouch, csrrowptr, csrcolval, csrnzval,
+        acc.is, acc.js, acc.vs,
+    )
+    if sort
+        # `csrcolval` is no longer needed, so we can reuse it
+        perm = resize!(csrcolval, length(basis))
+        sortperm!(perm, basis; kwargs...)
+        permute!(matrix, perm, perm), permute!(basis, perm)
+    else
+        return matrix, basis
+    end
+end
+
+"""
+    basis_breadth_first_search(::Type{Builder}, operator, starting_basis)
+
+Internal function that performs breadth-first search (BFS) on an operator.
+
+`Builder` is either [`MatrixBuilder`](@ref) or [`BasisBuilder`](@ref), which triggers
+building a matrix and basis, or only the basis of addresses, respectively.
+"""
+function basis_breadth_first_search(
+    ::Type{Builder}, operator, basis::Vector{A};
+    min_batch_size=100,
+    max_tasks=4 * Threads.nthreads(),
+
+    max_depth=Inf,
+    minimum_size=Inf,
+
+    cutoff=nothing,
+    filter=isnothing(cutoff) ? Returns(true) : a -> diagonal_element(operator, a) ≤ cutoff,
+
+    sizelim=Inf,
+    nnzs=0,
+    col_hint=0,
+
+    sort=false,
+    by=identity,
+    rev=false,
+    lt=isless,
+    order=Base.Forward,
+) where {Builder,A}
+
+    dim = dimension(operator, basis[1])
+    if dim > sizelim
+        throw(ArgumentError("dimension $dim larger than sizelim $sizelim"))
+    end
+    sizehint!(basis, nnzs)
+    sort_kwargs = (; by, rev, lt, order)
+
+    # addresses already visited and addresses visited in the last layer
+    seen = Set(basis)
+    last_seen = empty(seen)
+
+    # addresses to visit and their index
+    curr_frontier = collect(zip(basis, eachindex(basis)))
+    next_frontier = empty(curr_frontier)
+
+    # map from address to index
+    mapping = Dict(zip(basis, eachindex(basis)))
+
+    # builders store task-local storage and choose whether a matrix is to be built or not
+    builders = [Builder(; col_hint) for _ in 1:max_tasks]
+    result_accumulator = init_accumulator(Builder)
+
+    depth = 0
+    early_stop = false
+    while !isempty(curr_frontier)
+        depth += 1
+        early_stop = length(basis) > minimum_size || depth > max_depth
+
+        # We can stop here when not constructing the matrix. If we are, we still need
+        # to do another round to evaluate the columns.
+        if Builder <: BasisBuilder && early_stop
+            break
+        end
+
+        # Split the workload into chunks and spawn a task for each. These now run
+        # asynchronously while the main thread continues.
+        split_frontier = UniformSplit(curr_frontier, min_batch_size, max_tasks)
+        tasks = map(enumerate(split_frontier)) do (i, sub_frontier)
+            Threads.@spawn builders[$i]($operator, $sub_frontier, seen)
+        end
+
+        # Procesess the tasks in order they were spawned. This is fine as we've used more
+        # tasks than threads.
+        for (task, builder) in zip(tasks, builders)
+            wait(task)
+
+            result = builder.frontier
+            for k in result
+                if k ∉ last_seen
+                    push!(last_seen, k)
+                    if !early_stop && (isnothing(filter) || filter(k))
+                        push!(basis, k)
+                        push!(next_frontier, (k, length(basis)))
+                        mapping[k] = length(basis)
+                    end
+                end
+            end
+
+            # update matrix if needed
+            update_accumulator!(result_accumulator, builder, mapping)
+        end
+
+        # This had to wait until now so we don't get into data races
+        union!(seen, last_seen)
+        curr_frontier, next_frontier = next_frontier, curr_frontier
+        empty!(next_frontier)
+        empty!(last_seen)
+    end
+
+    return finalize_accumulator!(result_accumulator, basis, sort; sort_kwargs...)
+end
+
+function _address_to_basis(operator, addr_or_basis)
+    if addr_or_basis isa AbstractVector || addr_or_basis isa Tuple
+        check_address_type(operator, eltype(addr_or_basis))
+        return collect(addr_or_basis)
+    else
+        check_address_type(operator, typeof(addr_or_basis))
+        return [addr_or_basis]
+    end
+end
+
+"""
+    build_basis(
+        ham, address=starting_address(ham);
+        cutoff, filter, sizelim, sort=false, kwargs...
+    ) -> basis
+    build_basis(ham, addresses::AbstractVector; kwargs...)
+
+Get all basis element of a linear operator `ham` that are reachable (via
+non-zero matrix elements) from the address `address`, returned as a vector.
+Instead of a single address, a vector of `addresses` can be passed.
+Does not return the matrix, for that purpose use [`BasisSetRepresentation`](@ref).
+
+Providing an energy cutoff will skip addresses with diagonal elements greater
+than `cutoff`. Alternatively, an arbitrary `filter` function can be used instead.
+Addresses passed as arguments are not filtered.
+
+Providing a `max_depth` will limit the size of the basis by only visiting addresses that are
+connected to the `starting_address` through `max_depth` hops through the
+Hamiltonian. Similarly, providing `minimum_size` will stop the bulding process after the
+basis reaches a length of at least `minimum_size`.
+
+A maximum basis size `sizelim` can be set which will throw an error if the expected
+dimension of `ham` is larger than `sizelim`. This may be useful when memory may be a
+concern. These options are disabled by default.
+
+!!! warning
+        The order the basis is returned in is arbitrary and non-deterministic. Use
+        `sort=true` if the ordering matters.
+
+"""
+function build_basis(operator, addr=starting_address(operator); sizelim=Inf, kwargs...)
+    basis = _address_to_basis(operator, addr)
+    return basis_breadth_first_search(BasisBuilder{eltype(basis)}, operator, basis; sizelim, kwargs...)
+end
+
+"""
+    build_sparse_matrix_from_LO(
+        ham, address=starting_address(ham);
+        cutoff, filter=nothing, nnzs, col_hint, sort=false, kwargs...
+    ) -> sparse_matrix, basis
+    build_sparse_matrix_from_LO(ham, addresses::AbstractVector; kwargs...)
+
+Create a sparse matrix `sparse_matrix` of all reachable matrix elements of a linear operator
+`ham` starting from `address`. Instead of a single address, a vector of `addresses` can be
+passed.  The vector `basis` contains the addresses of basis configurations.
+
+Providing the number `nnzs` of expected calculated matrix elements and `col_hint` for the
+estimated number of nonzero off-diagonal matrix elements in each matrix column may improve
+performance.
+
+Providing an energy cutoff will skip the columns and rows with diagonal elements greater
+than `cutoff`. Alternatively, an arbitrary `filter` function can be used instead. These are
+not enabled by default. To generate the matrix truncated to the subspace spanned by the
+`addresses`, use `filter = Returns(false)`.
+
+Providing a `max_depth` will limit the size of the matrix by only visiting addresses that
+are connected to the `starting_address` through `max_depth` hops through the
+Hamiltonian. Similarly, providing `minimum_size` will stop the bulding process after the
+basis reaches a length of at least `minimum_size`.
+
+Setting `sort` to `true` will sort the `basis` and order the matrix rows and columns
+accordingly. This is useful when the order of the columns matters, e.g. when comparing
+matrices. Any additional keyword arguments are passed on to `Base.sortperm`.
+
+!!! warning
+        The order of the returned rows and columns is arbitrary and non-deterministic. Use
+        `sort=true` if the ordering matters.
+
+See [`BasisSetRepresentation`](@ref).
+"""
+function build_sparse_matrix_from_LO(
+    operator, addr=starting_address(operator); sizelim=1e7, kwargs...
+)
+    basis = _address_to_basis(operator, addr)
+    T = eltype(operator)
+    return basis_breadth_first_search(
+        MatrixBuilder{eltype(basis),T,Int32}, operator, basis;
+        sizelim, kwargs...
+    )
+end
diff --git a/src/ExactDiagonalization/basis_fock.jl b/src/ExactDiagonalization/basis_fock.jl
new file mode 100644
index 000000000..648588918
--- /dev/null
+++ b/src/ExactDiagonalization/basis_fock.jl
@@ -0,0 +1,193 @@
+"""
+    build_basis(addr::AbstractFockAddress)
+    build_basis(::Type{<:AbstractFockAddress}) -> basis
+
+Return all possible Fock states of a given type as a vector. This method is _much_ faster
+than `build_basis(::AbstractHamiltonian, ...)`, but does not take matrix blocking into
+account. This version of `build_basis` accepts no additional arguments.
+
+All address types except [`OccupationNumberFS`](@ref Main.Rimu.OccupationNumberFS) are
+supported.
+
+Returns a sorted vector of length equal to the [`dimension`](@ref) of `addr`.
+
+See also [`AbstractFockAddress`](@ref).
+"""
+function build_basis(addr::AbstractFockAddress)
+    return build_basis(typeof(addr))
+end
+
+###
+### BoseFS
+###
+# Algorithm:
+# Given an address type of N particles in M modes,
+# For all n_1 ∈ (0 ... N)
+# - put n_1 particles in the first mode
+# - recursively call the function to put the remaining (N-n_1) particles into (M-1) modes.
+# The recursion is parallelised through spawning tasks.
+
+# this is equivalent to `dimension(BoseFS{n,m})`, but does not return a `BigInt`.
+_bose_dimension(n, m) = binomial(n + m - 1, n)
+
+function build_basis(::Type{<:BoseFS{N,M}}) where {N,M}
+    T = typeof(near_uniform(BoseFS{N,M}))
+    result = Vector{T}(undef, _bose_dimension(N, M))
+    _bose_basis!(result, (), 1, N, Val(M), 2 * Threads.nthreads())
+    return result
+end
+
+# result: the basis that is eventually returned
+# postfix: partially build ONR
+# index: the position the built address is written to in result
+# remaining_n: the number of particles to be placed in the remaining part of the ONR
+# M: the number of modes in the ONR left to fill
+# n_tasks: number of tasks to spawn to build the basis in parallel
+@inline function _bose_basis!(
+    result::Vector, postfix, index, remaining_n, ::Val{M}, n_tasks::Int
+) where {M}
+    @sync if M < 5 || remaining_n ≤ 1 || n_tasks ≤ 0
+        _bose_basis!(result, postfix, index, remaining_n, Val(M))
+    else
+        n_tasks ÷= 2
+        Threads.@spawn begin
+            _bose_basis!(result, $(0, postfix...), $index, $remaining_n, Val(M-1), $n_tasks)
+        end
+        index += _bose_dimension(remaining_n, M - 1)
+        Threads.@spawn begin
+            _bose_basis!(
+                result, $(1, postfix...), $index, $(remaining_n - 1), Val(M-1), $n_tasks
+            )
+        end
+        index += _bose_dimension(remaining_n - 1, M - 1)
+        for n in 2:remaining_n
+            _bose_basis!(result, (n, postfix...), index, remaining_n - n, Val(M-1))
+            index += _bose_dimension(remaining_n - n, M - 1)
+        end
+    end
+end
+
+@inline function _bose_basis!(
+    result::Vector{T}, postfix, index, remaining_n, ::Val{M}
+) where {M,T}
+    if remaining_n == 0
+        @inbounds result[index] = T((ntuple(Returns(0), Val(M))..., postfix...))
+    elseif remaining_n == 1
+        _basis_basecase_N1!(result, postfix, index, Val(M))
+    elseif M == 1
+        @inbounds result[index] = T((remaining_n, postfix...))
+    elseif M == 2
+        _bose_basis_basecase_M2!(result, postfix, index, remaining_n)
+    elseif M == 3
+        _bose_basis_basecase_M3!(result, postfix, index, remaining_n)
+    else
+        for n in 0:remaining_n
+            _bose_basis!(result, (n, postfix...), index, remaining_n - n, Val(M - 1))
+            index += _bose_dimension(remaining_n - n, M - 1)
+        end
+    end
+    return nothing
+end
+
+###
+### FermiFS
+###
+# The algorithm is similar to the BoseFS one.
+
+# this is equivalent to `dimension(FermiFS{n,m})`, but does not return a `BigInt`.
+_fermi_dimension(n, m) = binomial(m, n)
+
+function build_basis(::Type{<:FermiFS{N,M}}) where {N,M}
+    T = typeof(near_uniform(FermiFS{N,M}))
+    result = Vector{T}(undef, _fermi_dimension(N, M))
+    _fermi_basis!(result, (), 1, N, Val(M), 2 * Threads.nthreads())
+    return result
+end
+
+# result: the basis that is eventually returned
+# postfix: partially build ONR
+# index: the position the built address is written to in result
+# remaining_n: the number of particles to be placed in the remaining part of the ONR
+# M: the number of modes in the ONR left to fill
+# n_tasks: number of tasks to spawn to build the basis in parallel
+@inline function _fermi_basis!(
+    result::Vector, postfix, index, remaining_n, ::Val{M}, n_tasks
+) where {M}
+    @sync if M < 5 || remaining_n ≤ M || remaining_n == 1 || n_tasks ≤ 0
+        _fermi_basis!(result, postfix, index, remaining_n, Val(M))
+    else
+        n_tasks ÷= 2
+        Threads.@spawn begin
+            _fermi_basis!(result, $(0, postfix...), $index, $remaining_n, Val(M-1), $n_tasks)
+        end
+        index += _fermi_dimension(remaining_n, M - 1)
+        _fermi_basis!(result, (1, postfix...), index, remaining_n - 1, Val(M - 1), n_tasks)
+    end
+end
+
+@inline function _fermi_basis!(
+    result::Vector{T}, postfix, index, remaining_n, ::Val{M}
+) where {M,T}
+    @assert M ≥ remaining_n
+    if remaining_n == 0
+        @inbounds result[index] = T((ntuple(Returns(0), Val(M))..., postfix...))
+    elseif remaining_n == M
+        @inbounds result[index] = T((ntuple(Returns(1), Val(M))..., postfix...))
+    elseif remaining_n == 1
+        _basis_basecase_N1!(result, postfix, index, Val(M))
+    else
+        _fermi_basis!(result, (0, postfix...), index, remaining_n, Val(M - 1))
+        index += _fermi_dimension(remaining_n, M - 1)
+        _fermi_basis!(result, (1, postfix...), index, remaining_n - 1, Val(M - 1))
+    end
+    return nothing
+end
+
+###
+### CompositeFS
+###
+function build_basis(::Type{C}) where {T,C<:CompositeFS{<:Any,<:Any,<:Any,T}}
+    sub_results = map(build_basis, reverse(Tuple(T.parameters)))
+    result = Vector{C}(undef, prod(length, sub_results))
+    Threads.@threads for i in eachindex(result)
+        @inbounds result[i] = C(_collect_addrs(sub_results, i))
+    end
+    return result
+end
+
+@inline _collect_addrs(::Tuple{}, _) = ()
+@inline function _collect_addrs((v, vs...), i)
+    rest, curr = fldmod1(i, length(v))
+    return (_collect_addrs(vs, rest)..., v[curr])
+end
+
+###
+### Base cases
+###
+# Base case for 1 particle in M modes.
+@inline function _basis_basecase_N1!(
+    result::Vector{T}, postfix, index, ::Val{M}
+) where {M,T}
+    rest = ntuple(Returns(0), Val(M))
+    for k in 1:M
+        @inbounds result[index + k - 1] = T((setindex(rest, 1, k)..., postfix...))
+    end
+end
+# Base case for bosons with `remaining_n` particles in 2 modes.
+@inline function _bose_basis_basecase_M2!(
+    result::Vector{T}, postfix, index, remaining_n
+) where {T}
+    for n1 in 0:remaining_n
+        @inbounds result[index + n1] = T((remaining_n - n1, n1, postfix...))
+    end
+end
+# Base case for bosons with `remaining_n` particles in 3 modes.
+@inline function _bose_basis_basecase_M3!(
+    result::Vector{T}, postfix, index, remaining_n
+) where {T}
+    k = 0
+    for n1 in 0:remaining_n, n2 in 0:(remaining_n - n1)
+        @inbounds result[index + k] = T((remaining_n - n1 - n2, n2, n1, postfix...))
+        k += 1
+    end
+end
diff --git a/src/ExactDiagonalization/basis_set_representation.jl b/src/ExactDiagonalization/basis_set_representation.jl
index 10577655a..8942ffd9d 100644
--- a/src/ExactDiagonalization/basis_set_representation.jl
+++ b/src/ExactDiagonalization/basis_set_representation.jl
@@ -1,173 +1,16 @@
-"""
-    build_sparse_matrix_from_LO(
-        ham, address=starting_address(ham);
-        cutoff, filter=nothing, nnzs, col_hint, sort=false, kwargs...
-    ) -> sparse_matrix, basis
-    build_sparse_matrix_from_LO(ham, addresses::AbstractVector; kwargs...)
-
-Create a sparse matrix `sparse_matrix` of all reachable matrix elements of a linear operator `ham`
-starting from `address`. Instead of a single address, a vector of `addresses` can be passed.
-The vector `basis` contains the addresses of basis configurations.
-
-Providing the number `nnzs` of expected calculated matrix elements and `col_hint` for the
-estimated number of nonzero off-diagonal matrix elements in each matrix column may improve
-performance.
-
-Providing an energy cutoff will skip the columns and rows with diagonal elements greater
-than `cutoff`. Alternatively, an arbitrary `filter` function can be used instead. These are
-not enabled by default. To generate the matrix truncated to the subspace spanned by the
-`addresses`, use `filter = Returns(false)`.
-
-Setting `sort` to `true` will sort the `basis` and order the matrix rows and columns
-accordingly. This is useful when the order of the columns matters, e.g. when comparing
-matrices. Any additional keyword arguments are passed on to `Base.sortperm`.
-
-See [`BasisSetRepresentation`](@ref).
-"""
-function build_sparse_matrix_from_LO(
-    ham, addr_or_vec=starting_address(ham);
-    cutoff=nothing,
-    filter=isnothing(cutoff) ? nothing : (a -> diagonal_element(ham, a) ≤ cutoff),
-    nnzs=0, col_hint=0, # sizehints are opt-in
-    sort=false, kwargs...
-)
-    # Set up `adds` as queue of addresses. Also returned as the basis.
-    adds = addr_or_vec isa Union{AbstractArray,Tuple} ? [addr_or_vec...] : [addr_or_vec]
-
-    T = eltype(ham)
-    dict = Dict(add => i for (i, add) in enumerate(adds)) # Map from addresses to indices
-    col = Dict{Int,T}()       # Temporary column storage
-    sizehint!(col, col_hint)
-
-    is = Int[] # row indices
-    js = Int[] # column indice
-    vs = T[]   # non-zero values
-
-    sizehint!(is, nnzs)
-    sizehint!(js, nnzs)
-    sizehint!(vs, nnzs)
-
-    i = 0
-    while i < length(adds)
-        i += 1
-        add = adds[i]
-        push!(is, i)
-        push!(js, i)
-        push!(vs, diagonal_element(ham, add))
-
-        empty!(col)
-        for (off, v) in offdiagonals(ham, add)
-            iszero(v) && continue
-            j = get(dict, off, nothing)
-            if isnothing(j)
-                # Energy cutoff: remember skipped addresses, but avoid adding them to `adds`
-                if !isnothing(filter) && !filter(off)
-                    dict[off] = 0
-                    j = 0
-                else
-                    push!(adds, off)
-                    j = length(adds)
-                    dict[off] = j
-                end
-            end
-            if !iszero(j)
-                col[j] = get(col, j, zero(T)) + v
-            end
-        end
-        # Copy the column into the sparse matrix vectors.
-        for (j, v) in col
-            iszero(v) && continue
-            push!(is, i)
-            push!(js, j)
-            push!(vs, v)
-        end
-    end
-
-    matrix = sparse(js, is, vs, length(adds), length(adds))
-    if sort
-        perm = sortperm(adds; kwargs...)
-        return permute!(matrix, perm, perm), permute!(adds, perm)
-    else
-        return matrix, adds
-    end
-end
-
-"""
-    build_basis(
-        ham, address=starting_address(ham);
-        cutoff, filter, sizelim, sort=false, minimum_size, kwargs...
-    ) -> basis
-    build_basis(ham, addresses::AbstractVector; kwargs...)
-
-Get all basis element of a linear operator `ham` that are reachable (via
-non-zero matrix elements) from the address `address`, returned as a vector.
-Instead of a single address, a vector of `addresses` can be passed.
-Does not return the matrix, for that purpose use [`BasisSetRepresentation`](@ref).
-
-Providing an energy cutoff will skip addresses with diagonal elements greater
-than `cutoff`. Alternatively, an arbitrary `filter` function can be used instead.
-Addresses passed as arguments are not filtered.
-A maximum basis size `sizelim` can be set which will throw an error if the expected
-dimension of `ham` is larger than `sizelim`. This may be useful when memory may be a
-concern.
-Setting a `minimum_size` will stop generating addresses once at least `minimum_size`
-addresses have been generated, rather than returning the full basis.
-These options are disabled by default.
-
-Setting `sort` to `true` will sort the basis. Any additional keyword arguments
-are passed on to `Base.sort!`.
-"""
-function build_basis(
-    ham, addr_or_vec=starting_address(ham);
-    cutoff=nothing,
-    filter=isnothing(cutoff) ? nothing : (a -> diagonal_element(ham, a) ≤ cutoff),
-    sort=false,
-    max_size=Inf, # retained for backwards compatibility; use sizelim instead
-    sizelim=max_size,
-    minimum_size=Inf,
-    kwargs...
-)
-    check_address_type(ham, addr_or_vec)
-    single_addr = addr_or_vec isa Union{AbstractArray,Tuple} ? addr_or_vec[1] : addr_or_vec
-    if minimum([dimension(ham, single_addr), minimum_size]) > sizelim
-        throw(ArgumentError("dimension larger than sizelim"))
-    end
-    # Set up `adds` as queue of addresses. Also returned as the basis.
-    adds = addr_or_vec isa Union{AbstractArray,Tuple} ? [addr_or_vec...] : [addr_or_vec]
-    known_basis = Set(adds)     # known addresses
-
-    i = 0
-    while i < length(adds) < minimum_size
-        i += 1
-        add = adds[i]
-
-        for (off, v) in offdiagonals(ham, add)
-            (iszero(v) || off in known_basis) && continue   # check if valid
-            push!(known_basis, off)
-            !isnothing(filter) && !filter(off) && continue  # check filter
-            push!(adds, off)
-        end
-    end
-
-    if sort
-        return sort!(adds, kwargs...)
-    else
-        return adds
-    end
-end
-
 """
     BasisSetRepresentation(
         hamiltonian::AbstractHamiltonian, addr=starting_address(hamiltonian);
-        sizelim=10^6, nnzs, cutoff, filter, sort=false, kwargs...
+        sizelim=10^7, cutoff, filter, max_depth, minimum_size, sort=false, kwargs...
     )
     BasisSetRepresentation(hamiltonian::AbstractHamiltonian, addresses::AbstractVector; kwargs...)
 
-Eagerly construct the basis set representation of the operator `hamiltonian` with all addresses
-reachable from `addr`. Instead of a single address, a vector of `addresses` can be passed.
+Eagerly construct the basis set representation of the operator `hamiltonian` with all
+addresses reachable from `addr`. Instead of a single address, a vector of `addresses` can be
+passed.
 
-An `ArgumentError` is thrown if `dimension(hamiltonian) > sizelim` in order to prevent memory
-overflow. Set `sizelim = Inf` in order to disable this behaviour.
+An `ArgumentError` is thrown if `dimension(hamiltonian) > sizelim` in order to prevent
+memory overflow. Set `sizelim = Inf` in order to disable this behaviour.
 
 Providing the number `nnzs` of expected calculated matrix elements and `col_hint` for the
 estimated number of nonzero off-diagonal matrix elements in each matrix column may improve
@@ -178,10 +21,19 @@ than `cutoff`. Alternatively, an arbitrary `filter` function can be used instead
 Addresses passed as arguments are not filtered. To generate the matrix truncated to the
 subspace spanned by the `addresses`, use `filter = Returns(false)`.
 
+Providing a `max_depth` will limit the size of the matrix and basis by only visiting
+addresses that are connected to the `starting_address` through `max_depth` hops through the
+Hamiltonian. Similarly, providing `minimum_size` will stop the bulding process after the
+basis reaches a length of at least `minimum_size`.
+
 Setting `sort` to `true` will sort the matrix rows and columns. This is useful when the
 order of the columns matters, e.g. when comparing matrices. Any additional keyword arguments
 are passed on to `Base.sortperm`.
 
+!!! warning
+        The order of the returned basis and matrix rows and columns is arbitrary and
+        non-deterministic. Use `sort=true` if the ordering matters.
+
 ## Fields
 * `sparse_matrix`: sparse matrix representing `hamiltonian` in the basis `basis`
 * `basis`: vector of addresses
@@ -192,15 +44,15 @@ are passed on to `Base.sortperm`.
 julia> hamiltonian = HubbardReal1D(BoseFS(1,0,0));
 
 julia> bsr = BasisSetRepresentation(hamiltonian)
-BasisSetRepresentation(HubbardReal1D(fs"|1 0 0⟩"; u=1.0, t=1.0)) with dimension 3 and 9 stored entries:3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 9 stored entries:
-  0.0  -1.0  -1.0
- -1.0   0.0  -1.0
- -1.0  -1.0   0.0
+BasisSetRepresentation(HubbardReal1D(fs"|1 0 0⟩"; u=1.0, t=1.0)) with dimension 3 and 6 stored entries:3×3 SparseArrays.SparseMatrixCSC{Float64, Int32} with 6 stored entries:
+   ⋅   -1.0  -1.0
+ -1.0    ⋅   -1.0
+ -1.0  -1.0    ⋅
 
 julia> BasisSetRepresentation(hamiltonian, bsr.basis[1:2]; filter = Returns(false)) # passing addresses and truncating
-BasisSetRepresentation(HubbardReal1D(fs"|1 0 0⟩"; u=1.0, t=1.0)) with dimension 2 and 4 stored entries:2×2 SparseArrays.SparseMatrixCSC{Float64, Int64} with 4 stored entries:
-  0.0  -1.0
- -1.0   0.0
+BasisSetRepresentation(HubbardReal1D(fs"|1 0 0⟩"; u=1.0, t=1.0)) with dimension 2 and 2 stored entries:2×2 SparseArrays.SparseMatrixCSC{Float64, Int32} with 2 stored entries:
+   ⋅   -1.0
+ -1.0    ⋅
 
 julia> using LinearAlgebra; round.(eigvals(Matrix(bsr)); digits = 4) # eigenvalues
 3-element Vector{Float64}:
@@ -262,13 +114,9 @@ end
 # default, does not enforce symmetries
 function _bsr_ensure_symmetry(
     ::LOStructure, hamiltonian::AbstractHamiltonian, addr_or_vec;
-    sizelim=10^6, test_approx_symmetry=true, kwargs...
+    test_approx_symmetry=true, kwargs...
 )
     single_addr = addr_or_vec isa Union{AbstractArray,Tuple} ? addr_or_vec[1] : addr_or_vec
-    d = dimension(hamiltonian, single_addr)
-    if d > sizelim
-        throw(ArgumentError("Dimension = $d larger than sizelim = $sizelim. Set a larger `sizelim` if this is safe."))
-    end
     check_address_type(hamiltonian, addr_or_vec)
     sparse_matrix, basis = build_sparse_matrix_from_LO(hamiltonian, addr_or_vec; kwargs...)
     return BasisSetRepresentation(sparse_matrix, basis, hamiltonian)
@@ -277,12 +125,9 @@ end
 # build the BasisSetRepresentation while enforcing hermitian symmetry
 function _bsr_ensure_symmetry(
     ::IsHermitian, hamiltonian::AbstractHamiltonian, addr_or_vec;
-    sizelim=10^6, test_approx_symmetry=true, kwargs...
+    test_approx_symmetry=true, kwargs...
 )
     single_addr = addr_or_vec isa Union{AbstractArray,Tuple} ? addr_or_vec[1] : addr_or_vec
-    if dimension(hamiltonian, single_addr) > sizelim
-        throw(ArgumentError("dimension larger than sizelim"))
-    end
     check_address_type(hamiltonian, addr_or_vec)
     sparse_matrix, basis = build_sparse_matrix_from_LO(hamiltonian, addr_or_vec; kwargs...)
     fix_approx_hermitian!(sparse_matrix; test_approx_symmetry) # enforce hermitian symmetry after building
@@ -291,6 +136,7 @@ end
 
 """
     fix_approx_hermitian!(A; test_approx_symmetry=true, kwargs...)
+
 Replaces the matrix `A` by `½(A + A')` in place. This will be successful and the result
 is guaranteed to pass the `ishermitian` test only if the matrix is square and already
 approximately hermitian.
@@ -325,6 +171,7 @@ end
 
 """
     isapprox_enforce_hermitian!(A::AbstractSparseMatrixCSC; kwargs...) -> Bool
+
 Checks whether the matrix `A` is approximately hermitian by checking each pair of transposed
 matrix elements with `isapprox`. Keyword arguments are passed on to `isapprox`.
 Returns boolean `true` is the test is passed and `false` if not.
diff --git a/src/ExactDiagonalization/exact_diagonalization_problem.jl b/src/ExactDiagonalization/exact_diagonalization_problem.jl
index 5814ffbd3..525923555 100644
--- a/src/ExactDiagonalization/exact_diagonalization_problem.jl
+++ b/src/ExactDiagonalization/exact_diagonalization_problem.jl
@@ -34,8 +34,10 @@ See [`BasisSetRepresentation`](@ref) for more information.
     sparse matrices and `10^5` for dense matrices.
 - `cutoff`: A cutoff value for the basis set representation.
 - `filter`: A filter function for the basis set representation.
-- `nnzs = 0`: The number of non-zero elements in the basis set representation. Setting a
-    non-zero value can speed up the computation.
+- `max_depth = Inf`: Limit the depth when building the matrix.
+- `minimum_size = Inf`: Stop building the matrix after this size is reached.
+- `nnzs = 0`: A hint for the number of non-zero elements in the basis set representation.
+  Setting a non-zero value can speed up the computation.
 - `col_hint = 0`: A hint for the number of columns in the basis set representation.
 - `sort = false`: Whether to sort the basis set representation.
 
diff --git a/src/ExactDiagonalization/init_and_solvers.jl b/src/ExactDiagonalization/init_and_solvers.jl
index a281bd1b5..62447dbf3 100644
--- a/src/ExactDiagonalization/init_and_solvers.jl
+++ b/src/ExactDiagonalization/init_and_solvers.jl
@@ -114,6 +114,8 @@ function CommonSolve.init(
     nnzs = get(kw, :nnzs, 0)
     col_hint = get(kw, :col_hint, 0)
     sort = get(kw, :sort, false)
+    max_depth = get(kw, :max_depth, Inf)
+    minimum_size = get(kw, :minimum_size, Inf)
 
     # determine the starting address or vector
     v0 = p.v0
@@ -136,11 +138,16 @@ function CommonSolve.init(
     @assert v0 isa Union{FrozenDVec{<:AbstractFockAddress},Nothing}
 
     # create the BasisSetRepresentation
-    bsr = BasisSetRepresentation(p.hamiltonian, addr_or_vec; sizelim, filter, nnzs, col_hint, sort)
+    bsr = BasisSetRepresentation(
+        p.hamiltonian, addr_or_vec;
+        sizelim, filter, nnzs, col_hint, sort, max_depth, minimum_size,
+    )
 
     # prepare kwargs for the solver
     kw = (; kw..., sizelim, cutoff, filter, nnzs, col_hint, sort)
-    kw_nt = delete(kw, (:sizelim, :cutoff, :filter, :nnzs, :col_hint, :sort))
+    kw_nt = delete(
+        kw, (:sizelim, :cutoff, :filter, :nnzs, :col_hint, :sort, :max_depth, :minimum_size)
+    )
 
     return MatrixEDSolver(algorithm, p, bsr, v0, kw_nt)
 end
diff --git a/src/Hamiltonians/MatrixHamiltonian.jl b/src/Hamiltonians/MatrixHamiltonian.jl
index ab2bb6883..4b5addd94 100644
--- a/src/Hamiltonians/MatrixHamiltonian.jl
+++ b/src/Hamiltonians/MatrixHamiltonian.jl
@@ -62,6 +62,7 @@ function offdiagonals(
     rows = rowvals(mh.m)[nzrange(mh.m, col)]
     vals = nonzeros(mh.m)[nzrange(mh.m, col)]
     drow = findprev(x->x==col, rows, min(col,length(rows))) # rows[drow]==col should be true
+    drow = isnothing(drow) ? 0 : drow
     return SparseMatrixOffdiagonals{Int,T,typeof(rows),typeof(vals)}(rows, vals, drow, col)
 end
 
@@ -72,10 +73,10 @@ struct SparseMatrixOffdiagonals{A,T,R,V} <: AbstractOffdiagonals{A,T}
     col::Int # colum of the matrix
 end
 function Base.getindex(smo::SparseMatrixOffdiagonals, chosen)
-    ind = ifelse(chosen < smo.drow, chosen, chosen + 1)
+    ind = ifelse(chosen < smo.drow, chosen, chosen + !iszero(smo.drow))
     return smo.rows[ind], smo.vals[ind]
 end
-Base.size(smo::SparseMatrixOffdiagonals) = (length(smo.rows) - 1,)
+Base.size(smo::SparseMatrixOffdiagonals) = (length(smo.rows) - !iszero(smo.drow),)
 
 function get_offdiagonal(
     mh::MatrixHamiltonian{<:Any, <:SparseArrays.AbstractSparseMatrixCSC},
diff --git a/test/BitStringAddresses.jl b/test/BitStringAddresses.jl
index 396f25032..1d7480e01 100644
--- a/test/BitStringAddresses.jl
+++ b/test/BitStringAddresses.jl
@@ -560,16 +560,20 @@ end
     @testset "OccupationNumberFS in hamiltonians" begin
         bfs = BoseFS(1, 2, 3)
         ofs = OccupationNumberFS(bfs)
-        @test sparse(HubbardMom1D(ofs)) == sparse(HubbardMom1D(bfs))
-        @test sparse(HubbardMom1DEP(ofs)) == sparse(HubbardMom1DEP(bfs))
-        @test sparse(HubbardReal1D(ofs)) == sparse(HubbardReal1D(bfs))
-        @test sparse(HubbardReal1DEP(ofs)) == sparse(HubbardReal1DEP(bfs))
-        @test sparse(HubbardRealSpace(ofs)) == sparse(HubbardRealSpace(bfs))
-        @test sparse(ExtendedHubbardReal1D(ofs)) == sparse(ExtendedHubbardReal1D(bfs))
+        for ham in (
+            HubbardMom1D,
+            HubbardMom1DEP,
+            HubbardReal1D,
+            HubbardReal1DEP,
+            HubbardRealSpace,
+            ExtendedHubbardReal1D,
+        )
+            @test sparse(ham(ofs); sort=true) == sparse(ham(bfs); sort=true)
+        end
         oham = HubbardReal1D(OccupationNumberFS(0, 2, 1))
         bham = HubbardReal1D(BoseFS(0, 2, 1))
-        @test sparse(ParitySymmetry(oham; odd=true)) ==
-            sparse(ParitySymmetry(bham; odd=true))
+        @test sparse(ParitySymmetry(oham; odd=true); sort=true) ==
+            sparse(ParitySymmetry(bham; odd=true); sort=true)
     end
 end
 
diff --git a/test/ExactDiagonalization.jl b/test/ExactDiagonalization.jl
index 3f011e885..5c8a2cff6 100644
--- a/test/ExactDiagonalization.jl
+++ b/test/ExactDiagonalization.jl
@@ -14,7 +14,7 @@ using Suppressor
         n = 10
         addr = near_uniform(BoseFS{n,m})
         ham = HubbardReal1D(addr)
-        bsr = BasisSetRepresentation(ham; nnzs=dimension(ham))
+        bsr = BasisSetRepresentation(ham)
         @test length(bsr.basis) == dimension(bsr) ≤ dimension(ham)
         @test_throws ArgumentError BasisSetRepresentation(ham, BoseFS((1, 2, 3))) # wrong address type
         @test Matrix(bsr) == Matrix(bsr.sparse_matrix) == Matrix(ham)
@@ -50,6 +50,25 @@ using Suppressor
         @test mat_cut == mat_cut_manual
     end
 
+    @testset "max_depth and minimum_size" begin
+        addr = BoseFS(5, 1 => 1)
+        ham = HubbardRealSpace(addr; geometry=CubicGrid((5,), (false,)))
+
+        basis_addr = build_basis(addr)
+
+        @test build_basis(ham; max_depth=0) == [addr]
+        @test build_basis(ham; max_depth=1) == basis_addr[1:2]
+        @test build_basis(ham; max_depth=2) == basis_addr[1:3]
+        @test build_basis(ham; max_depth=3) == basis_addr[1:4]
+        @test build_basis(ham; max_depth=4) == basis_addr
+
+        @test build_basis(ham; minimum_size=0) == [addr]
+        @test build_basis(ham; minimum_size=1) == basis_addr[1:2]
+        @test build_basis(ham; minimum_size=2) == basis_addr[1:3]
+        @test build_basis(ham; minimum_size=3) == basis_addr[1:4]
+        @test build_basis(ham; minimum_size=4) == basis_addr
+    end
+
     @testset "getindex" begin
         ham = HubbardReal1D(near_uniform(BoseFS{10,2}))
         bsr = BasisSetRepresentation(ham; sort=true)
@@ -128,6 +147,29 @@ using Suppressor
         basis = build_basis(ham, add; cutoff)
         @test basis == bsr.basis
     end
+
+    @testset "fock build basis" begin
+        for addr in (
+            BoseFS(10, 10),
+            FermiFS(1, 1, 0),
+            FermiFS(1, 0),
+            BoseFS(1, 1, 1, 1, 1, 2, 1, 1),
+            FermiFS(1, 1, 1, 0, 0, 0),
+            FermiFS2C((1, 1, 0, 0, 0, 0), (0, 1, 0, 1, 1, 0)),
+            CompositeFS(BoseFS(2, 0, 0), FermiFS(1, 0, 1), BoseFS(0, 2, 0), BoseFS(1, 0, 0)),
+        )
+            H = HubbardRealSpace(addr)
+            @test build_basis(addr) == build_basis(H; sort=true)
+        end
+        # These don't work with HubbardRealSpace
+        for addr in (
+            FermiFS(1),
+            BoseFS(1),
+        )
+            H = HubbardReal1D(addr)
+            @test build_basis(addr) == build_basis(H; sort=true)
+        end
+    end
 end
 
 Random.seed!(123) # for reproducibility, as some solvers start with random vectors