diff --git a/benchmarks/Project.toml b/benchmarks/Project.toml
new file mode 100644
index 0000000..84f0a8d
--- /dev/null
+++ b/benchmarks/Project.toml
@@ -0,0 +1,11 @@
+[deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+TenSolver = "cfd60579-01bf-4b47-9b3b-7c8dc91a23e2"
+
+[compat]
+BenchmarkTools = "1"
+TenSolver = "0.2"
+julia = "1.10"
diff --git a/benchmarks/knapsack/KnapsackBenchmark.jl b/benchmarks/knapsack/KnapsackBenchmark.jl
new file mode 100644
index 0000000..ce7d0a8
--- /dev/null
+++ b/benchmarks/knapsack/KnapsackBenchmark.jl
@@ -0,0 +1,497 @@
+module KnapsackBenchmark
+
+using DelimitedFiles: writedlm
+using LinearAlgebra: dot
+using Random: AbstractRNG, MersenneTwister, Random, rand
+
+import BenchmarkTools
+import TenSolver
+
+export projection_scaling_rows, run_benchmarks, write_csv
+
+const RNG_SEED = 66
+const INSTANCE_RNG = MersenneTwister(RNG_SEED)
+Random.seed!(RNG_SEED)
+
+# -----------------------------------------------------------------------------
+# Types and infrastructure
+# -----------------------------------------------------------------------------
+
+struct KnapsackInstance
+  name::String
+  weights::Vector{Int}
+  values::Vector{Float64}
+  capacity::Int
+
+  function KnapsackInstance(name, weights, values, capacity)
+    length(weights) == length(values) ||
+      throw(DimensionMismatch("weights and values must have the same length"))
+    isempty(weights) &&
+      throw(ArgumentError("a knapsack instance must contain at least one item"))
+    all(>(0), weights) || throw(ArgumentError("item weights must be positive integers"))
+    all(>(0), values) || throw(ArgumentError("item values must be positive"))
+    capacity >= 0 || throw(ArgumentError("capacity must be nonnegative"))
+    return new(String(name), Int.(weights), Float64.(values), Int(capacity))
+  end
+end
+
+struct BenchmarkSettings
+  iterations::Int
+  reads::Int
+  cutoff::Float64
+  time_limit::Float64
+  timing_samples::Int
+end
+
+struct KnapsackMethod
+  name::String
+  penalty_factor::Union{Missing, Float64}
+  formulate::Function
+  solve::Function
+end
+
+struct SuiteEntry
+  instance_key::String
+  method_key::String
+  outputs::Vector{Any}
+end
+
+struct BenchmarkReport
+  metadata::Vector{Pair{String, String}}
+  rows::Vector{NamedTuple}
+end
+
+function validate(settings::BenchmarkSettings, penalty_factors)
+  settings.iterations > 0 || throw(ArgumentError("iterations must be positive"))
+  settings.reads > 0 || throw(ArgumentError("reads must be positive"))
+  settings.cutoff > 0 || throw(ArgumentError("cutoff must be positive"))
+  settings.time_limit > 0 || throw(ArgumentError("time_limit must be positive"))
+  settings.timing_samples > 0 || throw(ArgumentError("timing_samples must be positive"))
+  isodd(settings.timing_samples) ||
+    throw(ArgumentError("timing_samples must be odd"))
+  all(>(0), penalty_factors) ||
+    throw(ArgumentError("penalty factors must be positive"))
+  return nothing
+end
+
+function runtime_metadata(settings::BenchmarkSettings, benchmark_id)
+  return [
+    "benchmark_id" => string(benchmark_id),
+    "julia_version" => string(VERSION),
+    "julia_threads" => string(Threads.nthreads()),
+    "system" => string(Sys.KERNEL),
+    "architecture" => string(Sys.ARCH),
+    "tensolver_version" => string(Base.pkgversion(TenSolver)),
+    "benchmarktools_version" => string(Base.pkgversion(BenchmarkTools)),
+    "rng_seed" => string(RNG_SEED),
+    "iterations" => string(settings.iterations),
+    "reads" => string(settings.reads),
+    "cutoff" => string(settings.cutoff),
+    "time_limit_seconds" => string(settings.time_limit),
+    "timing_samples" => string(settings.timing_samples),
+  ]
+end
+
+function write_csv(io::IO, rows)
+  isempty(rows) && return nothing
+  columns = propertynames(first(rows))
+  table = Matrix{Any}(undef, length(rows) + 1, length(columns))
+  table[1, :] .= string.(columns)
+  for (i, row) in enumerate(rows), (j, column) in enumerate(columns)
+    value = getproperty(row, column)
+    table[i + 1, j] = ismissing(value) ? "" : value
+  end
+  writedlm(io, table, ',')
+  return nothing
+end
+
+function write_csv(io::IO, report::BenchmarkReport)
+  for (key, value) in report.metadata
+    println(io, "# $(key)=$(value)")
+  end
+  return write_csv(io, report.rows)
+end
+
+# -----------------------------------------------------------------------------
+# Instance definitions and building helpers
+# -----------------------------------------------------------------------------
+
+reference_instance() =
+  KnapsackInstance("reference_4", [4, 3, 2, 3], [8, 4, 5, 3], 6)
+
+function pisinger_instance(
+  rng::AbstractRNG,
+  kind,
+  n;
+  coefficient_range = 10,
+)
+  n > 0 || throw(ArgumentError("instance size must be positive"))
+  coefficient_range >= 10 ||
+    throw(ArgumentError("coefficient range must be at least 10"))
+
+  weights = rand(rng, 1:coefficient_range, n)
+  correlation_range = div(coefficient_range, 10)
+  values = if kind == :uncorrelated
+    rand(rng, 1:coefficient_range, n)
+  elseif kind == :weakly_correlated
+    max.(1, weights .+ rand(rng, (-correlation_range):correlation_range, n))
+  elseif kind == :strongly_correlated
+    weights .+ correlation_range
+  elseif kind == :subset_sum
+    copy(weights)
+  else
+    throw(ArgumentError("unsupported Pisinger instance class: $(repr(kind))"))
+  end
+
+  capacity = max(maximum(weights), div(sum(weights), 2))
+  return KnapsackInstance("pisinger_$(kind)_n$(n)", weights, values, capacity)
+end
+
+function default_instances(rng::AbstractRNG = copy(INSTANCE_RNG))
+  specifications = (
+    (:uncorrelated, 8),
+    (:weakly_correlated, 12),
+    (:strongly_correlated, 16),
+    (:subset_sum, 16),
+  )
+  generated = map(specifications) do (kind, n)
+    return pisinger_instance(rng, kind, n)
+  end
+  return [reference_instance(), generated...]
+end
+
+item_weight(instance::KnapsackInstance, items) = dot(instance.weights, items)
+item_value(instance::KnapsackInstance, items) = dot(instance.values, items)
+is_capacity_feasible(instance::KnapsackInstance, items) =
+  item_weight(instance, items) <= instance.capacity
+
+function brute_force_optimum(instance::KnapsackInstance)
+  best = (value = -Inf, weight = typemax(Int), items = Int[])
+  for assignment in Iterators.product(fill(0:1, length(instance.weights))...)
+    items = collect(assignment)
+    weight = item_weight(instance, items)
+    value = item_value(instance, items)
+    if weight <= instance.capacity &&
+       (value > best.value || (value == best.value && weight < best.weight))
+      best = (; value, weight, items)
+    end
+  end
+  return best
+end
+
+function slack_weights(capacity::Integer)
+  capacity >= 0 || throw(ArgumentError("capacity must be nonnegative"))
+  encoded = Int[]
+  remaining = Int(capacity)
+  power = 1
+  while remaining > 0
+    weight = min(power, remaining)
+    push!(encoded, weight)
+    remaining -= weight
+    power *= 2
+  end
+  return encoded
+end
+
+function penalty_qubo(instance::KnapsackInstance, penalty::Real)
+  penalty > 0 || throw(ArgumentError("penalty must be positive"))
+  slack = slack_weights(instance.capacity)
+  coefficients = Float64.([instance.weights; slack])
+  values = [instance.values; zeros(length(slack))]
+  lambda = Float64(penalty)
+  Q = lambda .* (coefficients * coefficients')
+  l = -values .- (2lambda * instance.capacity) .* coefficients
+  constant = lambda * instance.capacity^2
+  return (; Q, l, constant, nitems = length(instance.weights))
+end
+
+penalty_value(model, assignment) =
+  dot(assignment, model.Q, assignment) +
+  dot(model.l, assignment) +
+  model.constant
+
+item_bits(sample, nitems) = round.(Int, sample[1:nitems])
+
+function best_sample(samples, decode, rank)
+  best = decode(first(samples))
+  best_rank = rank(best)
+  for sample in Iterators.drop(samples, 1)
+    candidate = decode(sample)
+    candidate_rank = rank(candidate)
+    if candidate_rank < best_rank
+      best = candidate
+      best_rank = candidate_rank
+    end
+  end
+  return best
+end
+
+function projection_method()
+  formulate = function (instance)
+    return TenSolver.SumConstraint(
+      collect(eachindex(instance.weights)),
+      instance.weights,
+      instance.capacity;
+      relation = :(<=),
+    )
+  end
+  solve = function (instance, constraint, settings)
+    reported_objective, solution = TenSolver.maximize(
+      instance.values;
+      constraints = [constraint],
+      solver_options(settings)...,
+    )
+    decode = sample -> item_bits(sample, length(instance.weights))
+    rank = items -> (
+      !is_capacity_feasible(instance, items),
+      -item_value(instance, items),
+      item_weight(instance, items),
+    )
+    items = best_sample(TenSolver.sample(solution, settings.reads), decode, rank)
+    return (;
+      reported_objective,
+      solution,
+      items,
+      nvariables = length(instance.weights),
+      penalty = missing,
+      penalized_objective = missing,
+    )
+  end
+  return KnapsackMethod("projection", missing, formulate, solve)
+end
+
+function penalty_method(penalty_factor)
+  factor = Float64(penalty_factor)
+  formulate = function (instance)
+    return penalty_qubo(instance, factor * sum(instance.values))
+  end
+  solve = function (instance, model, settings)
+    reported_objective, solution = TenSolver.minimize(
+      model.Q,
+      model.l,
+      model.constant;
+      solver_options(settings)...,
+    )
+    decode = function (sample)
+      assignment = round.(Int, sample)
+      return (; assignment, items = item_bits(assignment, model.nitems))
+    end
+    rank = candidate -> (
+      penalty_value(model, candidate.assignment),
+      -item_value(instance, candidate.items),
+      item_weight(instance, candidate.items),
+    )
+    best = best_sample(TenSolver.sample(solution, settings.reads), decode, rank)
+    return (;
+      reported_objective,
+      solution,
+      items = best.items,
+      nvariables = length(best.assignment),
+      penalty = factor * sum(instance.values),
+      penalized_objective = penalty_value(model, best.assignment),
+    )
+  end
+  return KnapsackMethod("penalty_$(factor)", factor, formulate, solve)
+end
+
+function benchmark_methods(penalty_factors)
+  return (projection_method(), (penalty_method(factor) for factor in penalty_factors)...)
+end
+
+function projection_scaling_instances()
+  probes = NamedTuple[]
+  for capacity in (1, 2, 4, 8)
+    nitems = 16
+    instance =
+      KnapsackInstance("capacity_$(capacity)", ones(Int, nitems), ones(Int, nitems), capacity)
+    push!(probes, (sweep = "capacity", instance))
+  end
+  for nitems in (8, 16, 32)
+    instance =
+      KnapsackInstance("items_$(nitems)", ones(Int, nitems), ones(Int, nitems), 3)
+    push!(probes, (sweep = "item_count", instance))
+  end
+  for scale in (4, 8, 32, 128)
+    weights = vcat([1, 2, 3], fill(scale, 5))
+    instance =
+      KnapsackInstance("weight_scale_$(scale)", weights, ones(Int, length(weights)), 3)
+    push!(probes, (sweep = "weight_magnitude", instance))
+  end
+  return probes
+end
+
+# -----------------------------------------------------------------------------
+# Runners
+# -----------------------------------------------------------------------------
+
+function solver_options(settings::BenchmarkSettings)
+  return (
+    iterations = settings.iterations,
+    time_limit = settings.time_limit,
+    cutoff = settings.cutoff,
+    inidim = 8,
+    maxdim = [10, 20, 40, 80, 120, 200],
+    noise = [1e-6, 1e-8, 0.0],
+    check_variance_every_iteration = typemax(Int),
+    vtol = -Inf,
+    verbosity = 0,
+  )
+end
+
+function result_row(
+  instance::KnapsackInstance,
+  exact,
+  method::KnapsackMethod,
+  output,
+  settings::BenchmarkSettings,
+)
+  stats = output.solution.stats
+  bonds = stats.max_bonds
+  elapsed = isempty(stats.elapsed_times) ? 0.0 : last(stats.elapsed_times)
+  feasible = is_capacity_feasible(instance, output.items)
+  value = item_value(instance, output.items)
+  return (
+    instance = instance.name,
+    method = method.name,
+    nitems = length(instance.weights),
+    nvariables = output.nvariables,
+    capacity = instance.capacity,
+    penalty_factor = method.penalty_factor,
+    penalty = output.penalty,
+    exact_value = exact.value,
+    original_value = value,
+    feasible,
+    optimality_gap = feasible ? exact.value - value : missing,
+    penalized_objective = output.penalized_objective,
+    solver_reported_objective = output.reported_objective,
+    sweeps = length(stats.energies),
+    time_limit_reached =
+      length(stats.energies) < settings.iterations && elapsed >= settings.time_limit,
+    solver_elapsed_seconds = elapsed,
+    solution_max_bond = isempty(stats.bond_dims) ? 0 : maximum(stats.bond_dims),
+    objective_mpo_bond = bonds.objective,
+    projection_mpo_bond =
+      isempty(bonds.projections) ? missing : maximum(bonds.projections),
+    effective_hamiltonian_bond = bonds.hamiltonian,
+  )
+end
+
+function execute_case(
+  instance::KnapsackInstance,
+  exact,
+  method::KnapsackMethod,
+  settings::BenchmarkSettings,
+)
+  model = method.formulate(instance)
+  output = method.solve(instance, model, settings)
+  return result_row(instance, exact, method, output, settings)
+end
+
+function benchmark_suite(instances, methods, settings::BenchmarkSettings)
+  suite = BenchmarkTools.BenchmarkGroup()
+  entries = SuiteEntry[]
+  samples = settings.timing_samples
+  for instance in instances
+    suite[instance.name] = BenchmarkTools.BenchmarkGroup()
+    exact = brute_force_optimum(instance)
+    for method in methods
+      outputs = Any[]
+      runner = function ()
+        push!(outputs, execute_case(instance, exact, method, settings))
+        return nothing
+      end
+      suite[instance.name][method.name] =
+        BenchmarkTools.@benchmarkable $runner() samples=samples evals=1 seconds=3600
+      push!(entries, SuiteEntry(instance.name, method.name, outputs))
+    end
+  end
+  return suite, entries
+end
+
+function measured_rows(trials, entries)
+  rows = NamedTuple[]
+  for entry in entries
+    trial = trials[entry.instance_key][entry.method_key]
+    sample_count = length(trial.times)
+    length(entry.outputs) >= sample_count ||
+      error("BenchmarkTools did not capture every measured output")
+    outputs = last(entry.outputs, sample_count)
+    median_index = sortperm(trial.times)[cld(sample_count, 2)]
+    row = merge(
+      outputs[median_index],
+      (end_to_end_wall_seconds = trial.times[median_index] / 1e9,),
+    )
+    push!(rows, row)
+  end
+  return rows
+end
+
+"""
+Run the penalty-versus-projection benchmark as a declarative BenchmarkTools
+suite. Pass keywords directly when calling this function from the REPL.
+"""
+function run_benchmarks(
+  instances = default_instances();
+  penalty_factors = (0.001, 0.01, 0.1, 1.1),
+  iterations = 6,
+  reads = 64,
+  cutoff = 1e-10,
+  time_limit = 120.0,
+  timing_samples = 3,
+  benchmark_id = "unversioned",
+  verbose = true,
+)
+  settings = BenchmarkSettings(
+    Int(iterations),
+    Int(reads),
+    Float64(cutoff),
+    Float64(time_limit),
+    Int(timing_samples),
+  )
+  validate(settings, penalty_factors)
+  suite, entries = benchmark_suite(
+    instances,
+    benchmark_methods(penalty_factors),
+    settings,
+  )
+  trials = BenchmarkTools.run(suite; verbose)
+  return BenchmarkReport(
+    runtime_metadata(settings, benchmark_id),
+    measured_rows(trials, entries),
+  )
+end
+
+"""
+Build the projection resource table through TenSolver's public solver API.
+"""
+function projection_scaling_rows(; cutoff = 1e-10, time_limit = 60.0)
+  settings = BenchmarkSettings(1, 1, Float64(cutoff), Float64(time_limit), 1)
+  return map(projection_scaling_instances()) do probe
+    instance = probe.instance
+    constraint = TenSolver.SumConstraint(
+      collect(eachindex(instance.weights)),
+      instance.weights,
+      instance.capacity;
+      relation = :(<=),
+    )
+    _, solution = TenSolver.maximize(
+      instance.values;
+      constraints = [constraint],
+      solver_options(settings)...,
+    )
+    bonds = solution.stats.max_bonds
+    return (
+      sweep = probe.sweep,
+      instance = instance.name,
+      nitems = length(instance.weights),
+      capacity = instance.capacity,
+      max_weight = maximum(instance.weights),
+      capacity_state_bound = instance.capacity + 2,
+      objective_mpo_bond = bonds.objective,
+      projection_mpo_bond = maximum(bonds.projections),
+      effective_hamiltonian_bond = bonds.hamiltonian,
+    )
+  end
+end
+
+end
diff --git a/benchmarks/knapsack/README.md b/benchmarks/knapsack/README.md
new file mode 100644
index 0000000..406ffdd
--- /dev/null
+++ b/benchmarks/knapsack/README.md
@@ -0,0 +1,99 @@
+# Knapsack penalty vs projection benchmark
+
+This benchmark compares two encodings of the same binary knapsack instances:
+
+- a conventional unconstrained QUBO with bounded-binary slack variables and a
+  squared capacity penalty;
+- TenSolver's native `SumConstraint`, lowered to a hard projection MPO.
+
+The comparison follows CoTenN's direct-constraint tensor-network framing
+([Sharma et al., PLDI 2026](https://doi.org/10.1145/3808272)).
+
+## Run
+
+From the repository root, instantiate the environment once:
+
+```bash
+julia --project=benchmarks -e 'using Pkg; Pkg.instantiate()'
+```
+
+Run the default benchmark:
+
+```bash
+julia --project=benchmarks benchmarks/knapsack/run.jl
+```
+
+The runner writes a CSV report to standard output. To change the workload or
+save the report, call `run_benchmarks` from the REPL:
+
+```julia
+include("benchmarks/knapsack/KnapsackBenchmark.jl")
+using .KnapsackBenchmark
+
+report = run_benchmarks(
+  [KnapsackBenchmark.reference_instance()];
+  penalty_factors = (0.1, 1.1),
+  iterations = 3,
+  reads = 8,
+  timing_samples = 1,
+  time_limit = 30,
+  benchmark_id = "pilot",
+)
+
+open("results.csv", "w") do io
+  write_csv(io, report)
+end
+```
+
+`run_benchmarks` builds a `BenchmarkTools.BenchmarkGroup` with one entry for
+each instance and formulation. Each entry measures the complete formulation,
+solve, and sampling path with one evaluation per sample.
+
+## Instances and methods
+
+The workload contains a four-item hand-checkable instance plus uncorrelated,
+weakly correlated, strongly correlated, and subset-sum classes from Martello,
+Pisinger, and Toth's standard 0-1 knapsack generator
+([paper](https://doi.org/10.1287/mnsc.45.3.414),
+[generator archive](https://hjemmesider.diku.dk/~pisinger/codes.html)).
+
+The generated instances use one shared RNG and contain 8, 12, or 16 items.
+Integer weights and capacity make the instances compatible with the projection
+formulation, and their small size permits exact brute-force reference values.
+
+For each instance, the suite benchmarks the projection formulation and penalty
+factors `0.001`, `0.01`, `0.1`, and `1.1` times the sum of item values. Every
+result is scored against the original knapsack objective and exact feasible
+optimum.
+
+## Resource scaling
+
+Generate the projection resource table from the REPL:
+
+```julia
+include("benchmarks/knapsack/KnapsackBenchmark.jl")
+using .KnapsackBenchmark
+
+write_csv(stdout, projection_scaling_rows())
+```
+
+It controls capacity, item count, and weight magnitude in separate sweeps. The
+table reports the projected Hamiltonian bond separately because it is the
+network used by constrained DMRG. Use the main solver CSV's
+`solution_max_bond` column to compare this controlled projection scaling with
+the observed penalty-QUBO DMRG bond growth.
+
+## Output
+
+The report header records run-wide provenance and settings once: runtime and
+package versions, thread count, system and architecture, RNG seed, sweep/read
+limits, tensor cutoff, timing samples, and benchmark identifier.
+
+Each result row records only case-specific information:
+
+- instance, formulation, capacity, item count, and encoded variable count;
+- penalty factor and coefficient for penalty-QUBO rows;
+- original value, feasibility, and exact optimality gap;
+- BenchmarkTools end-to-end time and solver-reported elapsed time;
+- completed sweeps and soft-time-limit status;
+- solution, objective, projection, and projected-Hamiltonian bond dimensions.
diff --git a/benchmarks/knapsack/run.jl b/benchmarks/knapsack/run.jl
new file mode 100644
index 0000000..7522663
--- /dev/null
+++ b/benchmarks/knapsack/run.jl
@@ -0,0 +1,5 @@
+include("KnapsackBenchmark.jl")
+
+using .KnapsackBenchmark
+
+write_csv(stdout, run_benchmarks())
diff --git a/docs/src/examples.md b/docs/src/examples.md
index 735ea26..77a338a 100644
--- a/docs/src/examples.md
+++ b/docs/src/examples.md
@@ -245,8 +245,7 @@ See [`TenSolver.SolverStatistics`](@ref) for more details.
 `psi.stats.energies`, `psi.stats.bond_dims`, and `psi.stats.elapsed_times` contain
 one value per completed iteration. `psi.stats.variances` has the same length and
 contains either the checked variance or `nothing` when that iteration did not
-perform the configured variance check. The former top-level properties
-`psi.energies`, `psi.bond_dims`, and `psi.elapsed_times` are deprecated aliases.
+perform the configured variance check.
 
 For per-iteration sampling, pass an `on_iteration` callback.
 The callback receives the MPS for that iteration alongside metadata as keyword arguments.
diff --git a/src/backends/dmrg.jl b/src/backends/dmrg.jl
index 2d94317..0c4039f 100644
--- a/src/backends/dmrg.jl
+++ b/src/backends/dmrg.jl
@@ -237,7 +237,7 @@ function minimize_mpo( H_obj :: MPO
                      , iterations :: Union{Nothing, Int} = nothing
                      , time_limit = +Inf
                      , vtol       = cutoff
-                     , check_variance_every_iteration = 10
+                     , check_variance_every_iteration :: Int = 10
                      # DMRG keywords
                      , inidim     = 40
                      , maxdim     = [10, 10, 10, 20, 50, 100, 100, 200, 300, 300, 400, 400, 800, 900, 1000]
@@ -252,6 +252,9 @@ function minimize_mpo( H_obj :: MPO
                      , permutation :: Vector{Int} = collect(1:length(H_obj))
                      ) where {T}
   callback_every >= 1 || throw(ArgumentError("`callback_every` must be >= 1, got $callback_every"))
+  check_variance_every_iteration >= 1 || throw(ArgumentError(
+    "`check_variance_every_iteration` must be >= 1, got $check_variance_every_iteration",
+  ))
   initial_time = time()
 
   # Quantization
diff --git a/src/solution.jl b/src/solution.jl
index f429cd4..835d3c6 100644
--- a/src/solution.jl
+++ b/src/solution.jl
@@ -84,24 +84,6 @@ struct Solution{T <: Real}
   end
 end
 
-function Base.getproperty(solution::Solution, name::Symbol)
-  if name === :energies || name === :bond_dims || name === :elapsed_times
-    Base.depwarn(
-      "`solution.$name` is deprecated; use `solution.stats.$name` instead.",
-      name,
-    )
-    return getproperty(getfield(solution, :stats), name)
-  end
-
-  return getfield(solution, name)
-end
-
-function Base.propertynames(::Solution, _private::Bool=false)
-  fields = fieldnames(Solution)
-  aliases = (:energies, :bond_dims, :elapsed_times)
-  return (fields..., aliases...)
-end
-
 function infeasible_solution(::Type{T}, domain, stats) where {T <: Real}
   return Solution{T}(nothing, domain, Int[], stats)
 end
diff --git a/src/solver.jl b/src/solver.jl
index d822fca..dc12db3 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -106,6 +106,8 @@ Keyword arguments:
 - `preprocess :: Bool` - Defaults to `false`. If `true`, permute QUBO variables before constructing the MPS Hamiltonian
   so coupled variables are closer in the one-dimensional tensor order. Samples are returned in the
   caller's original variable order. This is an experimental feature and may be subject to changes.
+- `check_variance_every_iteration :: Int` - Calculate and record the Hamiltonian
+  variance every N iterations. Must be >= 1. Defaults to `10`.
 - `on_iteration :: Function` - Called after each recorded iteration as
   `f(psi::MPS; iteration, objective, bond_dim, elapsed_time)`.
   `objective` is the expected objective function ⟨ψ|H|ψ⟩ at this iteration.
@@ -124,8 +126,6 @@ Keyword arguments:
   may have limited support depending on the backend.
 
 The returned `Solution` carries per-iteration convergence data in `solution.stats`.
-The former top-level fields `solution.energies`, `solution.bond_dims`, and
-`solution.elapsed_times` remain available as deprecated aliases.
 
 Provably infeasible constrained models are reported as a status:
 `minimize` logs a warning and returns `+Inf` (the minimum over an empty feasible set)
diff --git a/test/constrained_solve.jl b/test/constrained_solve.jl
index bda3cac..8923fe8 100644
--- a/test/constrained_solve.jl
+++ b/test/constrained_solve.jl
@@ -187,11 +187,15 @@ end
   @testset "Infeasible constraints report status, not exception" begin
     impossible = AbstractConstraint[SumConstraint([1, 2], [1, 1], 3; relation=:(==))]
 
-    E, psi = @test_logs (:warn, r"empty feasible subspace") minimize(
-      zeros(2, 2);
-      constraints=impossible,
-      verbosity=0,
-    )
+    io = IOBuffer()
+    E, psi = with_logger(ConsoleLogger(io, Logging.Debug)) do
+      minimize(
+        zeros(2, 2);
+        constraints=impossible,
+        verbosity=0,
+      )
+    end
+    debug_log = String(take!(io))
 
     @test E == Inf
     @test !is_feasible(psi)
@@ -199,25 +203,8 @@ end
     @test isempty(psi.stats.bond_dims)
     @test isempty(psi.stats.elapsed_times)
     @test isempty(psi.stats.variances)
-    @test length(psi.stats.max_bonds.projections) == 1
-    @test all(>(0), psi.stats.max_bonds.projections)
-    @test psi.stats.max_bonds.objective > 0
     @test psi.stats.max_bonds.initial_state == 0
-    @test psi.stats.max_bonds.hamiltonian > 0
     @test [0, 0] ∉ psi
-    @test_throws DomainError TenSolver.sample(psi)
-
-    io = IOBuffer()
-    E_debug, psi_debug = with_logger(ConsoleLogger(io, Logging.Debug)) do
-      minimize(
-        zeros(2, 2);
-        constraints=impossible,
-        verbosity=0,
-      )
-    end
-    debug_log = String(take!(io))
-    @test E_debug == Inf
-    @test !is_feasible(psi_debug)
     @test occursin("empty feasible subspace", debug_log)
     @test !occursin("Exception while generating log record", debug_log)
 
diff --git a/test/qubo.jl b/test/qubo.jl
index 9a6b362..815b31f 100644
--- a/test/qubo.jl
+++ b/test/qubo.jl
@@ -127,7 +127,13 @@
 
   @testset "Iteration stats tracking" begin
     @testset "Solution carries stats" begin
-      E, psi = minimize([1.0 0; 0 -1.0]; iterations=5, verbosity = 0)
+      E, psi = minimize(
+        [1.0 0; 0 -1.0];
+        iterations=5,
+        check_variance_every_iteration=2,
+        vtol=-Inf,
+        verbosity=0,
+      )
       @test psi isa TenSolver.Solution
       @test length(psi.stats.energies)      == 5
       @test length(psi.stats.bond_dims)     == 5
@@ -136,39 +142,13 @@
       @test all(isfinite, psi.stats.energies)
       @test issorted(psi.stats.elapsed_times)
       @test all(>(0), psi.stats.bond_dims)
-      @test all(isnothing, psi.stats.variances)
+      @test map(isnothing, psi.stats.variances) == [true, false, true, false, true]
+      @test all(isfinite, filter(!isnothing, psi.stats.variances))
       @test isfinite(last(psi.stats.energies))
       @test last(psi.stats.energies) ≈ E
 
       @test isempty(psi.stats.max_bonds.projections)
-      @test psi.stats.max_bonds.objective > 0
-      @test psi.stats.max_bonds.initial_state > 0
-      @test psi.stats.max_bonds.hamiltonian > 0
-    end
-
-    @testset "Solution preserves deprecated stats aliases" begin
-      _, psi = minimize([1.0 0; 0 -1.0]; iterations=2, verbosity = 0)
-
-      @test :energies in propertynames(psi)
-      @test :bond_dims in propertynames(psi)
-      @test :elapsed_times in propertynames(psi)
-      @test_deprecated psi.energies === psi.stats.energies
-      @test_deprecated psi.bond_dims === psi.stats.bond_dims
-      @test_deprecated psi.elapsed_times === psi.stats.elapsed_times
-    end
-
-    @testset "Variance stats distinguish unchecked iterations" begin
-      _, psi = minimize(
-        [1.0 0; 0 -1.0];
-        iterations=3,
-        check_variance_every_iteration=1,
-        vtol=-Inf,
-        verbosity=0,
-      )
-
-      @test length(psi.stats.variances) == 3
-      @test all(!isnothing, psi.stats.variances)
-      @test all(isfinite, psi.stats.variances)
+      @test psi.stats.max_bonds.objective == psi.stats.max_bonds.hamiltonian
     end
 
     @testset "on_iteration callback preserves its documented signature" begin
@@ -190,9 +170,11 @@
       @test calls == [3, 6, 9]
     end
 
-    @testset "callback_every < 1 throws" begin
+    @testset "iteration cadence < 1 throws" begin
       @test_throws ArgumentError minimize([1.0 0; 0 -1.0]; verbosity = 0, iterations=1, callback_every=0)
       @test_throws ArgumentError minimize([1.0 0; 0 -1.0]; verbosity = 0, iterations=1, callback_every=-1)
+      @test_throws ArgumentError minimize([1.0 0; 0 -1.0]; verbosity = 0, iterations=1, check_variance_every_iteration=0)
+      @test_throws ArgumentError minimize([1.0 0; 0 -1.0]; verbosity = 0, iterations=1, check_variance_every_iteration=-1)
     end
 
     @testset "callback receives a fresh MPS each iteration" begin