diff --git a/src/Core.jl b/src/Core.jl
index bc9210aed..2fc9a55cb 100644
--- a/src/Core.jl
+++ b/src/Core.jl
@@ -17,6 +17,9 @@ using .OptionsStructModule:
     AbstractOptions,
     Options,
     ComplexityMapping,
+    AbstractParetoOptions,
+    ParetoSingleOptions,
+    ParetoTopKOptions,
     specialized_options,
     operator_specialization
 using .OperatorsModule:
diff --git a/src/ExpressionBuilder.jl b/src/ExpressionBuilder.jl
index 00264be6a..e22bebd6d 100644
--- a/src/ExpressionBuilder.jl
+++ b/src/ExpressionBuilder.jl
@@ -10,7 +10,7 @@ using DynamicExpressions:
     AbstractExpressionNode, AbstractExpression, constructorof, with_metadata
 using StatsBase: StatsBase
 using ..CoreModule: AbstractOptions, Dataset
-using ..HallOfFameModule: HallOfFame
+using ..HallOfFameModule: HallOfFame, ParetoSingle, ParetoTopK
 using ..PopulationModule: Population
 using ..PopMemberModule: PopMember
 
@@ -124,20 +124,32 @@ end
         pop::Population, options::AbstractOptions, dataset::Dataset{T,L}
     ) where {T,L}
         return Population(
-            map(Fix{2}(Fix{3}(embed_metadata, dataset), options), pop.members)
+            map(member -> embed_metadata(member, options, dataset), pop.members)
+        )
+    end
+    function embed_metadata(
+        el::ParetoSingle, options::AbstractOptions, dataset::Dataset{T,L}
+    ) where {T,L}
+        return ParetoSingle(embed_metadata(el.member, options, dataset))
+    end
+    function embed_metadata(
+        el::ParetoTopK, options::AbstractOptions, dataset::Dataset{T,L}
+    ) where {T,L}
+        return ParetoTopK(
+            map(member -> embed_metadata(member, options, dataset), el.members), el.k
         )
     end
     function embed_metadata(
         hof::HallOfFame, options::AbstractOptions, dataset::Dataset{T,L}
     ) where {T,L}
         return HallOfFame(
-            map(Fix{2}(Fix{3}(embed_metadata, dataset), options), hof.members), hof.exists
+            map(el -> embed_metadata(el, options, dataset), hof.elements), hof.exists
         )
     end
     function embed_metadata(
-        vec::Vector{H}, options::AbstractOptions, dataset::Dataset{T,L}
+        sets::Vector{H}, options::AbstractOptions, dataset::Dataset{T,L}
     ) where {T,L,H<:Union{HallOfFame,Population,PopMember}}
-        return map(Fix{2}(Fix{3}(embed_metadata, dataset), options), vec)
+        return map(set -> embed_metadata(set, options, dataset), sets)
     end
 end
 
@@ -171,11 +183,23 @@ function strip_metadata(
 ) where {T,L}
     return Population(map(member -> strip_metadata(member, options, dataset), pop.members))
 end
+function strip_metadata(
+    el::ParetoSingle, options::AbstractOptions, dataset::Dataset{T,L}
+) where {T,L}
+    return ParetoSingle(strip_metadata(el.member, options, dataset))
+end
+function strip_metadata(
+    el::ParetoTopK, options::AbstractOptions, dataset::Dataset{T,L}
+) where {T,L}
+    return ParetoTopK(
+        map(member -> strip_metadata(member, options, dataset), el.members), el.k
+    )
+end
 function strip_metadata(
     hof::HallOfFame, options::AbstractOptions, dataset::Dataset{T,L}
 ) where {T,L}
     return HallOfFame(
-        map(member -> strip_metadata(member, options, dataset), hof.members), hof.exists
+        map(el -> strip_metadata(el, options, dataset), hof.elements), hof.exists
     )
 end
 
diff --git a/src/HallOfFame.jl b/src/HallOfFame.jl
index 7d870f78b..dd3f6f887 100644
--- a/src/HallOfFame.jl
+++ b/src/HallOfFame.jl
@@ -2,15 +2,75 @@ module HallOfFameModule
 
 using StyledStrings: @styled_str
 using DynamicExpressions: AbstractExpression, string_tree
+using Printf: @sprintf
 using ..UtilsModule: split_string, AnnotatedIOBuffer, dump_buffer
+using ..CoreModule: ParetoSingleOptions, ParetoTopKOptions
 using ..CoreModule: AbstractOptions, Dataset, DATA_TYPE, LOSS_TYPE, relu, create_expression
 using ..ComplexityModule: compute_complexity
 using ..PopMemberModule: PopMember
 using ..InterfaceDynamicExpressionsModule: format_dimensions, WILDCARD_UNIT_STRING
-using Printf: @sprintf
+using ..PopulationModule: Population
+
+"""
+    AbstractParetoElement{P<:PopMember}
 
+Abstract type for storing elements on the Pareto frontier.
+
+# Subtypes
+- `ParetoSingle`: Stores a single member at each complexity level
+- `ParetoTopK`: Stores multiple members at each complexity level in a fixed-size bucket
 """
-    HallOfFame{T<:DATA_TYPE,L<:LOSS_TYPE}
+abstract type AbstractParetoElement{P<:PopMember} end
+
+pop_member_type(::Type{<:AbstractParetoElement{P}}) where {P} = P
+
+struct ParetoSingle{T,L,N,P<:PopMember{T,L,N}} <: AbstractParetoElement{P}
+    member::P
+end
+struct ParetoTopK{T,L,N,P<:PopMember{T,L,N}} <: AbstractParetoElement{P}
+    members::Vector{P}
+    k::Int
+end
+
+Base.copy(el::ParetoSingle) = ParetoSingle(copy(el.member))
+Base.copy(el::ParetoTopK) = ParetoTopK(sizehint!(map(copy, el.members), el.k + 1), el.k)
+
+Base.first(el::ParetoSingle) = el.member
+Base.first(el::ParetoTopK) = first(el.members)
+
+Base.iterate(el::ParetoSingle) = (el.member, nothing)
+Base.iterate(::ParetoSingle, ::Nothing) = nothing
+Base.iterate(el::ParetoTopK) = iterate(el.members)
+Base.iterate(el::ParetoTopK, state) = iterate(el.members, state)
+
+function Base.show(io::IO, mime::MIME"text/plain", el::ParetoSingle)
+    print(io, "ParetoSingle(")
+    show(io, mime, el.member)
+    print(io, ")")
+    return nothing
+end
+
+function _depwarn_pareto_single(funcsym::Symbol)
+    Base.depwarn(
+        "Hall of fame `.members` is now `.elements` which is a vector of `AbstractParetoElement` objects. ",
+        funcsym,
+    )
+    return nothing
+end
+
+@inline function Base.getproperty(s::ParetoSingle, name::Symbol)
+    name == :member && return getfield(s, :member)
+    _depwarn_pareto_single(:getproperty)
+    return getproperty(s.member, name)
+end
+@inline function Base.setproperty!(s::ParetoSingle, name::Symbol, value)
+    name == :member && return setfield!(s, :member, value)
+    _depwarn_pareto_single(:setproperty!)
+    return setproperty!(s.member, name, value)
+end
+
+"""
+    HallOfFame{T<:DATA_TYPE,L<:LOSS_TYPE,N<:AbstractExpression{T}}
 
 List of the best members seen all time in `.members`, with `.members[c]` being
 the best member seen at complexity c. Including only the members which actually
@@ -22,23 +82,39 @@ have been set, you can run `.members[exists]`.
     These are ordered by complexity, with `.members[1]` the member with complexity 1.
 - `exists::Array{Bool,1}`: Whether the member at the given complexity has been set.
 """
-struct HallOfFame{T<:DATA_TYPE,L<:LOSS_TYPE,N<:AbstractExpression{T}}
-    members::Array{PopMember{T,L,N},1}
-    exists::Array{Bool,1} #Whether it has been set
+struct HallOfFame{
+    T<:DATA_TYPE,
+    L<:LOSS_TYPE,
+    N<:AbstractExpression{T},
+    H<:AbstractParetoElement{<:PopMember{T,L,N}},
+}
+    elements::Vector{H}
+    exists::Vector{Bool}
+end
+pop_member_type(::Type{<:HallOfFame{T,L,N,H}}) where {T,L,N,H} = pop_member_type(H)
+@inline function Base.getproperty(hof::HallOfFame, name::Symbol)
+    if name == :members
+        Base.depwarn(
+            "HallOfFame.members is deprecated. Use HallOfFame.elements instead.",
+            :getproperty,
+        )
+        return getfield(hof, :elements)
+    end
+    return getfield(hof, name)
 end
 function Base.show(io::IO, mime::MIME"text/plain", hof::HallOfFame{T,L,N}) where {T,L,N}
     println(io, "HallOfFame{...}:")
-    for i in eachindex(hof.members, hof.exists)
-        s_member, s_exists = if hof.exists[i]
-            sprint((io, m) -> show(io, mime, m), hof.members[i]), "true"
+    for i in eachindex(hof.elements, hof.exists)
+        s_element, s_exists = if hof.exists[i]
+            sprint((io, m) -> show(io, mime, m), hof.elements[i]), "true"
         else
             "undef", "false"
         end
         println(io, " "^4 * ".exists[$i] = $s_exists")
-        print(io, " "^4 * ".members[$i] =")
-        splitted = split(strip(s_member), '\n')
+        print(io, " "^4 * ".elements[$i] =")
+        splitted = split(strip(s_element), '\n')
         if length(splitted) == 1
-            println(io, " " * s_member)
+            println(io, " " * s_element)
         else
             println(io)
             foreach(line -> println(io, " "^8 * line), splitted)
@@ -61,58 +137,153 @@ Arguments:
 - `dataset`: Dataset containing the input data.
 """
 function HallOfFame(
-    options::AbstractOptions, dataset::Dataset{T,L}
+    options::AbstractOptions, dataset::Dataset{T,L};
 ) where {T<:DATA_TYPE,L<:LOSS_TYPE}
     base_tree = create_expression(zero(T), options, dataset)
+    N = typeof(base_tree)
+    member = PopMember(
+        base_tree, L(0), L(Inf), options; parent=-1, deterministic=options.deterministic
+    )
 
-    return HallOfFame{T,L,typeof(base_tree)}(
+    return HallOfFame(
         [
-            PopMember(
-                copy(base_tree),
-                L(0),
-                L(Inf),
-                options;
-                parent=-1,
-                deterministic=options.deterministic,
-            ) for i in 1:(options.maxsize)
+            init_pareto_element(options.pareto_element_options, member) for
+            i in 1:(options.maxsize)
         ],
         [false for i in 1:(options.maxsize)],
     )
 end
+Base.copy(hof::HallOfFame) = HallOfFame(map(copy, hof.elements), copy(hof.exists))
 
-function Base.copy(hof::HallOfFame)
-    return HallOfFame(
-        [copy(member) for member in hof.members], [exists for exists in hof.exists]
-    )
+function init_pareto_element(::Union{ParetoSingleOptions,ParetoSingle}, member::PopMember)
+    return ParetoSingle(copy(member))
+end
+function init_pareto_element(opt::Union{ParetoTopKOptions,ParetoTopK}, member::PopMember)
+    members = sizehint!(typeof(member)[], opt.k + 1)
+    push!(members, copy(member))
+    return ParetoTopK(members, opt.k)
+end
+
+function Base.push!(hof::HallOfFame, (size, member)::Pair{<:Integer,<:PopMember})
+    maxsize = length(hof.elements)
+    if 0 < size <= maxsize
+        if !hof.exists[size]
+            hof.elements[size] = init_pareto_element(hof.elements[size], member)
+            hof.exists[size] = true
+        else
+            hof.elements[size] = push!(hof.elements[size], member.score => member)
+        end
+    end
+    return hof
+end
+
+function Base.push!(el::ParetoSingle, (score, member)::Pair{<:LOSS_TYPE,<:PopMember})
+    return el.member.score > score ? ParetoSingle(copy(member)) : el
+end
+function Base.push!(el::ParetoTopK, (score, member)::Pair{<:LOSS_TYPE,<:PopMember})
+    if isempty(el.members)
+        push!(el.members, copy(member))
+        return el
+    elseif el.members[end].score <= score
+        # No update needed
+        return el
+    elseif el.members[1].score > score
+        pushfirst!(el.members, copy(member))
+    else
+        # Find the first member with worse score
+        i = findfirst(m -> m.score > score, el.members)::Int
+        # member assumes that position, and pushes the array forward
+        insert!(el.members, i, copy(member))
+    end
+
+    if length(el.members) > el.k
+        pop!(el.members)
+    end
+
+    return el
+end
+
+function Base.append!(hof::HallOfFame, pop::Population; options::AbstractOptions)
+    for member in pop.members
+        size = compute_complexity(member, options)
+        push!(hof, size => member)
+    end
+    return hof
+end
+
+function Base.merge!(hof1::HallOfFame, hof2::HallOfFame)
+    for i in eachindex(hof1.elements, hof1.exists, hof2.elements, hof2.exists)
+        if hof1.exists[i] && hof2.exists[i]
+            hof1.elements[i] = merge(hof1.elements[i], hof2.elements[i])
+        elseif !hof1.exists[i] && hof2.exists[i]
+            hof1.elements[i] = copy(hof2.elements[i])
+            hof1.exists[i] = true
+        else
+            # do nothing, as !hof2.exists[i]
+        end
+    end
+    return hof1
+end
+function Base.merge(el1::ParetoSingle, el2::ParetoSingle)
+    # Remember: we want the MIN score (bad API choice, but we're stuck with it for now)
+    return el1.member.score <= el2.member.score ? el1 : copy(el2)
+end
+function Base.merge(el1::ParetoTopK, el2::ParetoTopK)
+    P = pop_member_type(typeof(el1))
+    new_neighborhood = sizehint!(P[], el1.k + 1)
+    i1 = firstindex(el1.members)
+    n1 = length(el1.members)
+    i2 = firstindex(el2.members)
+    n2 = length(el2.members)
+    i = 1
+    while i1 <= n1 && i2 <= n2 && i <= el1.k
+        m1 = el1.members[i1]
+        m2 = el2.members[i2]
+        if m1.score <= m2.score
+            # TODO: Is it safe that we don't copy here? I think so; since we are merging
+            #       onto el1 (see `Base.merge!`), but perhaps someone could misuse this.
+            push!(new_neighborhood, m1)
+            i1 += 1
+        else
+            push!(new_neighborhood, copy(m2))
+            i2 += 1
+        end
+        i += 1
+    end
+    return ParetoTopK(new_neighborhood, el1.k)
 end
 
 """
-    calculate_pareto_frontier(hallOfFame::HallOfFame{T,L,P}) where {T<:DATA_TYPE,L<:LOSS_TYPE}
+    calculate_pareto_frontier(hof::HallOfFame)
+
+Compute the dominating pareto curve - each returned member must be better than all simpler equations.
 """
-function calculate_pareto_frontier(hallOfFame::HallOfFame{T,L,N}) where {T,L,N}
-    # TODO - remove dataset from args.
-    P = PopMember{T,L,N}
-    # Dominating pareto curve - must be better than all simpler equations
+function calculate_pareto_frontier(hof::HallOfFame)
+    P = pop_member_type(typeof(hof))
     dominating = P[]
-    for size in eachindex(hallOfFame.members)
-        if !hallOfFame.exists[size]
+    for i in eachindex(hof.elements)
+        if !hof.exists[i]
             continue
         end
-        member = hallOfFame.members[size]
-        # We check if this member is better than all members which are smaller than it and
-        # also exist.
-        betterThanAllSmaller = true
-        for i in 1:(size - 1)
-            if !hallOfFame.exists[i]
+        element = hof.elements[i]
+        member = first(element)
+        # We check if this member is better than all
+        # elements which are smaller than it and also exist.
+        is_dominating = true
+        for j in 1:(i - 1)
+            if !hof.exists[j]
                 continue
             end
-            simpler_member = hallOfFame.members[i]
-            if member.loss >= simpler_member.loss
-                betterThanAllSmaller = false
+            smaller_element = hof.elements[j]
+            smaller_member = first(smaller_element)
+            if member.loss >= smaller_member.loss
+                is_dominating = false
                 break
             end
+            # TODO: Why are we using loss and not score? In other words,
+            #       why are we _pushing_ based on score and not loss?
         end
-        if betterThanAllSmaller
+        if is_dominating
             push!(dominating, copy(member))
         end
     end
diff --git a/src/MLJInterface.jl b/src/MLJInterface.jl
index 2c37c4583..33def62bc 100644
--- a/src/MLJInterface.jl
+++ b/src/MLJInterface.jl
@@ -27,7 +27,14 @@ using LossFunctions: SupervisedLoss
 using ..InterfaceDynamicQuantitiesModule: get_dimensions_type
 using ..InterfaceDynamicExpressionsModule: InterfaceDynamicExpressionsModule as IDE
 using ..CoreModule:
-    Options, Dataset, AbstractMutationWeights, MutationWeights, LOSS_TYPE, ComplexityMapping
+    Options,
+    Dataset,
+    AbstractMutationWeights,
+    MutationWeights,
+    LOSS_TYPE,
+    ComplexityMapping,
+    AbstractParetoOptions,
+    ParetoSingleOptions
 using ..CoreModule.OptionsModule: DEFAULT_OPTIONS, OPTION_DESCRIPTIONS
 using ..ComplexityModule: compute_complexity
 using ..HallOfFameModule: HallOfFame, format_hall_of_fame
diff --git a/src/Options.jl b/src/Options.jl
index b6f1dd847..afaa256a6 100644
--- a/src/Options.jl
+++ b/src/Options.jl
@@ -29,7 +29,8 @@ using ..OperatorsModule:
     safe_atanh
 using ..MutationWeightsModule: AbstractMutationWeights, MutationWeights, mutations
 import ..OptionsStructModule: Options
-using ..OptionsStructModule: ComplexityMapping, operator_specialization
+using ..OptionsStructModule:
+    ComplexityMapping, operator_specialization, AbstractParetoOptions, ParetoSingleOptions
 using ..UtilsModule: @save_kwargs, @ignore
 
 """Build constraints on operator-level complexity from a user-passed dict."""
@@ -523,6 +524,7 @@ $(OPTION_DESCRIPTIONS)
     ###           hof_migration
     ###           fraction_replaced
     ###           fraction_replaced_hof
+    ###           pareto_element_options
     ###           topn
     ## 9. Data Preprocessing:
     ###           [none]
@@ -581,6 +583,7 @@ $(OPTION_DESCRIPTIONS)
     hof_migration::Bool=true,
     fraction_replaced::Union{Real,Nothing}=nothing,
     fraction_replaced_hof::Union{Real,Nothing}=nothing,
+    pareto_element_options::AbstractParetoOptions=ParetoSingleOptions(),
     topn::Union{Nothing,Integer}=nothing,
     ## 9. Data Preprocessing:
     ## 10. Stopping Criteria:
@@ -873,6 +876,7 @@ $(OPTION_DESCRIPTIONS)
         expression_type,
         typeof(expression_options),
         typeof(set_mutation_weights),
+        typeof(pareto_element_options),
         turbo,
         bumper,
         deprecated_return_state::Union{Bool,Nothing},
@@ -913,6 +917,7 @@ $(OPTION_DESCRIPTIONS)
         ncycles_per_iteration,
         fraction_replaced,
         fraction_replaced_hof,
+        pareto_element_options,
         topn,
         verbosity,
         Val(print_precision),
diff --git a/src/OptionsStruct.jl b/src/OptionsStruct.jl
index 169b30574..f9ff86313 100644
--- a/src/OptionsStruct.jl
+++ b/src/OptionsStruct.jl
@@ -8,6 +8,22 @@ using LossFunctions: SupervisedLoss
 
 import ..MutationWeightsModule: AbstractMutationWeights
 
+"""
+    AbstractParetoOptions
+
+Abstract type for different Pareto front storage strategies.
+
+# Subtypes
+- `ParetoSingleOptions`: Store only the single best member at each complexity
+- `ParetoTopKOptions`: Store multiple members at each complexity in a fixed-size bucket
+"""
+abstract type AbstractParetoOptions end
+
+struct ParetoSingleOptions <: AbstractParetoOptions end
+Base.@kwdef struct ParetoTopKOptions <: AbstractParetoOptions
+    k::Int
+end
+
 """
 This struct defines how complexity is calculated.
 
@@ -184,6 +200,7 @@ struct Options{
     E<:AbstractExpression,
     EO<:NamedTuple,
     MW<:AbstractMutationWeights,
+    PO<:AbstractParetoOptions,
     _turbo,
     _bumper,
     _return_state,
@@ -224,6 +241,7 @@ struct Options{
     ncycles_per_iteration::Int
     fraction_replaced::Float32
     fraction_replaced_hof::Float32
+    pareto_element_options::PO
     topn::Int
     verbosity::Union{Int,Nothing}
     v_print_precision::Val{print_precision}
diff --git a/src/SearchUtils.jl b/src/SearchUtils.jl
index 89ca03ee2..e5a1cdd25 100644
--- a/src/SearchUtils.jl
+++ b/src/SearchUtils.jl
@@ -13,11 +13,12 @@ using Logging: AbstractLogger
 
 using DynamicExpressions: AbstractExpression, string_tree
 using ..UtilsModule: subscriptify
-using ..CoreModule: Dataset, AbstractOptions, Options, RecordType, max_features
+using ..CoreModule:
+    Dataset, AbstractOptions, AbstractParetoOptions, Options, RecordType, max_features
 using ..ComplexityModule: compute_complexity
 using ..PopulationModule: Population
 using ..PopMemberModule: PopMember
-using ..HallOfFameModule: HallOfFame, string_dominating_pareto_curve
+using ..HallOfFameModule: HallOfFame, string_dominating_pareto_curve, init_pareto_element
 using ..ProgressBarsModule: WrappedProgressBar, manually_iterate!, barlen
 using ..AdaptiveParsimonyModule: RunningSearchStatistics
 
@@ -231,6 +232,17 @@ end
 
 const DefaultWorkerOutputType{P,H} = Tuple{P,H,RecordType,Float64}
 
+function get_hall_of_fame_type(
+    ::Type{T},
+    ::Type{L},
+    example_ex::AbstractExpression,
+    pareto_element_options::AbstractParetoOptions,
+) where {T,L}
+    example_member = PopMember(example_ex, zero(L), zero(L); deterministic=false)
+    example_pareto_element = init_pareto_element(pareto_element_options, example_member)
+    ParetoElementType = typeof(example_pareto_element)
+    return HallOfFame{T,L,typeof(example_ex),ParetoElementType}
+end
 function get_worker_output_type(
     ::Val{PARALLELISM}, ::Type{PopType}, ::Type{HallOfFameType}
 ) where {PARALLELISM,PopType,HallOfFameType}
@@ -356,8 +368,10 @@ function _check_for_loss_threshold(_, ::Nothing, ::AbstractOptions)
 end
 function _check_for_loss_threshold(halls_of_fame, f::F, options::AbstractOptions) where {F}
     return all(halls_of_fame) do hof
-        any(hof.members[hof.exists]) do member
-            f(member.loss, compute_complexity(member, options))::Bool
+        any(hof.elements[hof.exists]) do element
+            any(element) do member
+                f(member.loss, compute_complexity(member, options))::Bool
+            end
         end
     end
 end
@@ -518,7 +532,7 @@ end
 load_saved_population(::Nothing; kws...) = nothing
 
 """
-    AbstractSearchState{T,L,N}
+    AbstractSearchState{T,L,N,H}
 
 An abstract type encapsulating the internal state of the search process during symbolic regression.
 
@@ -535,17 +549,18 @@ Look through the source of `equation_search` to see how this is used.
 - [`AbstractOptions`](@ref SymbolicRegression.CoreModule.OptionsStruct.AbstractOptions): See how to extend abstract types for customizing options.
 
 """
-abstract type AbstractSearchState{T,L,N<:AbstractExpression{T}} end
+abstract type AbstractSearchState{T,L,N<:AbstractExpression{T},H<:HallOfFame{T,L,N}} end
 
 """
-    SearchState{T,L,N,WorkerOutputType,ChannelType} <: AbstractSearchState{T,L,N}
+    SearchState{T,L,N,H,WorkerOutputType,ChannelType} <: AbstractSearchState{T,L,N,H}
 
 The state of the search, including the populations, worker outputs, tasks, and
 channels. This is used to manage the search and keep track of runtime variables
 in a single struct.
 """
-Base.@kwdef struct SearchState{T,L,N<:AbstractExpression{T},WorkerOutputType,ChannelType} <:
-                   AbstractSearchState{T,L,N}
+Base.@kwdef struct SearchState{
+    T,L,N<:AbstractExpression{T},H<:HallOfFame{T,L,N},WorkerOutputType,ChannelType
+} <: AbstractSearchState{T,L,N,H}
     procs::Vector{Int}
     we_created_procs::Bool
     worker_output::Vector{Vector{WorkerOutputType}}
@@ -553,7 +568,7 @@ Base.@kwdef struct SearchState{T,L,N<:AbstractExpression{T},WorkerOutputType,Cha
     channels::Vector{Vector{ChannelType}}
     worker_assignment::WorkerAssignments
     task_order::Vector{Tuple{Int,Int}}
-    halls_of_fame::Vector{HallOfFame{T,L,N}}
+    halls_of_fame::Vector{H}
     last_pops::Vector{Vector{Population{T,L,N}}}
     best_sub_pops::Vector{Vector{Population{T,L,N}}}
     all_running_search_statistics::Vector{RunningSearchStatistics}
@@ -676,22 +691,4 @@ function construct_datasets(
     ]
 end
 
-function update_hall_of_fame!(
-    hall_of_fame::HallOfFame, members::Vector{PM}, options::AbstractOptions
-) where {PM<:PopMember}
-    for member in members
-        size = compute_complexity(member, options)
-        valid_size = 0 < size <= options.maxsize
-        if !valid_size
-            continue
-        end
-        not_filled = !hall_of_fame.exists[size]
-        better_than_current = member.score < hall_of_fame.members[size].score
-        if not_filled || better_than_current
-            hall_of_fame.members[size] = copy(member)
-            hall_of_fame.exists[size] = true
-        end
-    end
-end
-
 end
diff --git a/src/SingleIteration.jl b/src/SingleIteration.jl
index 2d36e6c87..a7c2bd9fc 100644
--- a/src/SingleIteration.jl
+++ b/src/SingleIteration.jl
@@ -56,9 +56,9 @@ function s_r_cycle(
         num_evals += tmp_num_evals
         for (i, member) in enumerate(pop.members)
             size = compute_complexity(member, options)
-            score = if options.batching
+            if options.batching
                 oid = member.tree
-                if loss_cache[i].oid != oid || first_loop
+                same_batch_score = if loss_cache[i].oid != oid || first_loop
                     # Evaluate on fixed batch so that we can more accurately
                     # compare expressions with a batched loss (though the batch
                     # changes each iteration, and we evaluate on full-batch outside,
@@ -73,8 +73,14 @@ function s_r_cycle(
                     # the cached score
                     loss_cache[i].score
                 end
+
+                stored_member = copy(member)
+                stored_member.score = same_batch_score
+                # ^I think we don't want to _store_ this score in this member
+                #  since it would fundamentally mutate the population
+                push!(best_examples_seen, size => stored_member)
             else
-                member.score
+                push!(best_examples_seen, size => member)
             end
             # TODO: Note that this per-population hall of fame only uses the batched
             #       loss, and is therefore inaccurate. Therefore, some expressions
@@ -83,13 +89,6 @@ function s_r_cycle(
             # - Could just recompute losses here (expensive)
             # - Average over a few batches
             # - Store multiple expressions in hall of fame
-            if 0 < size <= options.maxsize && (
-                !best_examples_seen.exists[size] ||
-                score < best_examples_seen.members[size].score
-            )
-                best_examples_seen.exists[size] = true
-                best_examples_seen.members[size] = copy(member)
-            end
         end
         first_loop = false
     end
diff --git a/src/SymbolicRegression.jl b/src/SymbolicRegression.jl
index fd95b0851..0845efcca 100644
--- a/src/SymbolicRegression.jl
+++ b/src/SymbolicRegression.jl
@@ -234,6 +234,9 @@ using .CoreModule:
     ComplexityMapping,
     AbstractMutationWeights,
     MutationWeights,
+    AbstractParetoOptions,
+    ParetoSingleOptions,
+    ParetoTopKOptions,
     get_safe_op,
     max_features,
     is_weighted,
@@ -282,7 +285,12 @@ using .LossFunctionsModule: eval_loss, score_func, update_baseline_loss!
 using .PopMemberModule: PopMember, reset_birth!
 using .PopulationModule: Population, best_sub_pop, record_population, best_of_sample
 using .HallOfFameModule:
-    HallOfFame, calculate_pareto_frontier, string_dominating_pareto_curve
+    HallOfFame,
+    ParetoSingle,
+    ParetoTopK,
+    calculate_pareto_frontier,
+    string_dominating_pareto_curve,
+    init_pareto_element
 using .MutateModule: mutate!, condition_mutation_weights!, MutationResult
 using .SingleIterationModule: s_r_cycle, optimize_and_simplify_population
 using .ProgressBarsModule: WrappedProgressBar
@@ -296,6 +304,7 @@ using .SearchUtilsModule:
     WorkerAssignments,
     DefaultWorkerOutputType,
     assign_next_worker!,
+    get_hall_of_fame_type,
     get_worker_output_type,
     extract_from_worker,
     @sr_spawner,
@@ -317,7 +326,6 @@ using .SearchUtilsModule:
     construct_datasets,
     save_to_file,
     get_cur_maxsize,
-    update_hall_of_fame!,
     logging_callback!
 using .LoggingModule: AbstractSRLogger, SRLogger, get_logger
 using .TemplateExpressionModule: TemplateExpression, TemplateStructure
@@ -424,10 +432,10 @@ which is useful for debugging and profiling.
 
 # Returns
 - `hallOfFame::HallOfFame`: The best equations seen during the search.
-    hallOfFame.members gives an array of `PopMember` objects, which
-    have their tree (equation) stored in `.tree`. Their score (loss)
-    is given in `.score`. The array of `PopMember` objects
-    is enumerated by size from `1` to `options.maxsize`.
+    map(first, hallOfFame.elements) gives an array of `PopMember` objects
+    These have their tree (equation) stored in `.tree`. Their score (loss)
+    is given in `.score`. The array is enumerated by size from `1`
+    to `options.maxsize`.
 """
 function equation_search(
     X::AbstractMatrix{T},
@@ -591,9 +599,9 @@ end
     nout = length(datasets)
     example_dataset = first(datasets)
     example_ex = create_expression(zero(T), options, example_dataset)
-    NT = typeof(example_ex)
-    PopType = Population{T,L,NT}
-    HallOfFameType = HallOfFame{T,L,NT}
+    ExpressionType = typeof(example_ex)
+    PopType = Population{T,L,ExpressionType}
+    HallOfFameType = get_hall_of_fame_type(T, L, example_ex, options.pareto_element_options)
     WorkerOutputType = get_worker_output_type(
         Val(ropt.parallelism), PopType, HallOfFameType
     )
@@ -650,7 +658,7 @@ end
         j in 1:nout
     ]
 
-    return SearchState{T,L,typeof(example_ex),WorkerOutputType,ChannelType}(;
+    return SearchState{T,L,typeof(example_ex),HallOfFameType,WorkerOutputType,ChannelType}(;
         procs=procs,
         we_created_procs=we_created_procs,
         worker_output=worker_output,
@@ -670,12 +678,12 @@ end
     )
 end
 function _initialize_search!(
-    state::AbstractSearchState{T,L,N},
+    state::AbstractSearchState{T,L,N,H},
     datasets,
     ropt::AbstractRuntimeOptions,
     options::AbstractOptions,
     saved_state,
-) where {T,L,N}
+) where {T,L,N,H}
     nout = length(datasets)
 
     init_hall_of_fame = load_saved_hall_of_fame(saved_state)
@@ -688,12 +696,23 @@ function _initialize_search!(
         # case the dataset changed:
         for j in eachindex(init_hall_of_fame, datasets, state.halls_of_fame)
             hof = strip_metadata(init_hall_of_fame[j], options, datasets[j])
-            for member in hof.members[hof.exists]
-                score, result_loss = score_func(datasets[j], member, options)
+            new_hof = HallOfFame(options, datasets[j])
+            for el in hof.elements[hof.exists], member in el
+                size = compute_complexity(member, options)
+                score, result_loss = score_func(
+                    datasets[j], member, options; complexity=size
+                )
                 member.score = score
                 member.loss = result_loss
+
+                # In case the new loss changes the elements stored in the Pareto
+                # frontier, we push them individually:
+                push!(new_hof, size => member)
             end
-            state.halls_of_fame[j] = hof
+            # TODO: Would be nice if there was a way to mark `init_hall_of_fame[j]`
+            # as being dead now. We want to make sure we never use it at this point
+            # in the code.
+            state.halls_of_fame[j] = new_hof
         end
     end
 
@@ -748,11 +767,11 @@ function _initialize_search!(
     return nothing
 end
 function _warmup_search!(
-    state::AbstractSearchState{T,L,N},
+    state::AbstractSearchState{T,L,N,H},
     datasets,
     ropt::AbstractRuntimeOptions,
     options::AbstractOptions,
-) where {T,L,N}
+) where {T,L,N,H}
     nout = length(datasets)
     for j in 1:nout, i in 1:(options.populations)
         dataset = datasets[j]
@@ -769,9 +788,7 @@ function _warmup_search!(
         last_pop = state.worker_output[j][i]
         updated_pop = @sr_spawner(
             begin
-                in_pop = first(
-                    extract_from_worker(last_pop, Population{T,L,N}, HallOfFame{T,L,N})
-                )
+                in_pop = first(extract_from_worker(last_pop, Population{T,L,N}, H))
                 _dispatch_s_r_cycle(
                     in_pop,
                     dataset,
@@ -883,9 +900,10 @@ function _main_search_loop!(
                 update_frequencies!(state.all_running_search_statistics[j]; size)
             end
             #! format: off
-            update_hall_of_fame!(state.halls_of_fame[j], cur_pop.members, options)
-            update_hall_of_fame!(state.halls_of_fame[j], best_seen.members[best_seen.exists], options)
+            append!(state.halls_of_fame[j], cur_pop; options)
+            merge!(state.halls_of_fame[j], best_seen)
             #! format: on
+            # TODO: Confirm that `best_seen.members` have full-batch scores
 
             # Dominating pareto curve - must be better than all simpler equations
             dominating = calculate_pareto_frontier(state.halls_of_fame[j])
@@ -904,7 +922,12 @@ function _main_search_loop!(
                 )
             end
             if options.hof_migration && length(dominating) > 0
-                migrate!(dominating => cur_pop, options; frac=options.fraction_replaced_hof)
+                all_hof_members = [
+                    member for el in state.halls_of_fame[j].elements for member in el
+                ]
+                migrate!(
+                    all_hof_members => cur_pop, options; frac=options.fraction_replaced_hof
+                )
             end
             ###################################################################
 
@@ -1101,15 +1124,23 @@ end
         dataset, out_pop, options, cur_maxsize, record
     )
     num_evals += evals_from_optimize
-    if options.batching
-        for i_member in 1:(options.maxsize)
-            score, result_loss = score_func(dataset, best_seen.members[i_member], options)
-            best_seen.members[i_member].score = score
-            best_seen.members[i_member].loss = result_loss
+    return_hof = if options.batching
+        # Compute full-dataset scores for all members of the Pareto front.
+        new_hof = HallOfFame(options, dataset)::typeof(best_seen)
+        for el in best_seen.elements[best_seen.exists], member in el
+            size = compute_complexity(member, options)
+            score, result_loss = score_func(dataset, member, options)
+            member.score = score
+            member.loss = result_loss
             num_evals += 1
+
+            push!(new_hof, size => member)
         end
+        new_hof
+    else
+        best_seen
     end
-    return (out_pop, best_seen, record, num_evals)
+    return (out_pop, return_hof, record, num_evals)
 end
 function _info_dump(
     state::AbstractSearchState,
diff --git a/test/test_custom_operators_multiprocessing.jl b/test/test_custom_operators_multiprocessing.jl
index 22c978771..9b3eb3726 100644
--- a/test/test_custom_operators_multiprocessing.jl
+++ b/test/test_custom_operators_multiprocessing.jl
@@ -48,6 +48,6 @@ hof = equation_search(
 )
 
 @test any(
-    early_stop(member.loss, my_complexity(member.tree)) for
-    member in hof.members[hof.exists]
+    early_stop(member.loss, my_complexity(member.tree)) for el in hof.elements[hof.exists],
+    member in el
 )
diff --git a/test/test_early_stop.jl b/test/test_early_stop.jl
index 3ba36e555..f203402df 100644
--- a/test/test_early_stop.jl
+++ b/test/test_early_stop.jl
@@ -15,5 +15,6 @@ options = SymbolicRegression.Options(;
 hof = equation_search(X, y; options=options, niterations=1_000_000_000)
 
 @test any(
-    early_stop(member.loss, count_nodes(member.tree)) for member in hof.members[hof.exists]
+    early_stop(member.loss, count_nodes(member.tree)) for el in hof.elements[hof.exists],
+    member in el
 )
diff --git a/test/test_pretty_printing.jl b/test/test_pretty_printing.jl
index dcbc3f59c..86a3b260f 100644
--- a/test/test_pretty_printing.jl
+++ b/test/test_pretty_printing.jl
@@ -25,7 +25,7 @@ end
 
 @testitem "pretty print hall of fame" tags = [:part1] begin
     using SymbolicRegression
-    using SymbolicRegression: embed_metadata
+    using SymbolicRegression: embed_metadata, ParetoSingle
     using SymbolicRegression.CoreModule: safe_pow
 
     options = Options(; binary_operators=[+, safe_pow], maxsize=7)
@@ -46,24 +46,24 @@ end
 
     hof = HallOfFame(options, dataset)
     hof = embed_metadata(hof, options, dataset)
-    hof.members[5] = member
+    hof.elements[5] = ParetoSingle(member)
     hof.exists[5] = true
     s_hof = strip(shower(hof))
     true_s = "HallOfFame{...}:
     .exists[1] = false
-    .members[1] = undef
+    .elements[1] = undef
     .exists[2] = false
-    .members[2] = undef
+    .elements[2] = undef
     .exists[3] = false
-    .members[3] = undef
+    .elements[3] = undef
     .exists[4] = false
-    .members[4] = undef
+    .elements[4] = undef
     .exists[5] = true
-    .members[5] = PopMember(tree = ((x ^ 2.0) + 1.5), loss = 16.25, score = 1.0)
+    .elements[5] = ParetoSingle(PopMember(tree = ((x ^ 2.0) + 1.5), loss = 16.25, score = 1.0))
     .exists[6] = false
-    .members[6] = undef
+    .elements[6] = undef
     .exists[7] = false
-    .members[7] = undef"
+    .elements[7] = undef"
 
     @test s_hof == true_s
 end