Shift constrained solves below the projection kernel by bernalde · Pull Request #133 · SECQUOIA/TenSolver.jl

bernalde · 2026-08-20T18:30:31Z

Fix the DMRG collapse on feasible constrained models.

Closes #132.

Mechanism

The projected Hamiltonian P'HP assigns energy zero to the infeasible subspace (the kernel of the projections). When every feasible objective value is positive, that kernel is the ground space of the projected operator, so the DMRG sweep is attracted into it; noise and SVD truncation leak amplitude there, the post-sweep re-projection annihilates the state, and the solve dies with constrained DMRG produced a state with zero feasible amplitude. The mechanism is independent of per-variable domains — uniform-domain fixtures on main crash 50–60% of runs (see Evidence).

Fix

Constrained solves now run DMRG on H − λ·Id with λ one more than a cheap box bound on the objective magnitude (max|domain|²·Σ|Q| + max|domain|·Σ|l|, and the analogous per-term bound for polynomial objectives). Every in-domain assignment then has negative energy, the kernel (still at zero) sits above the whole feasible spectrum, and the sweep itself penalizes leaked amplitude back out instead of being attracted to it. The shift is added back to every reported energy, so returned optima, per-iteration stats, logs, and callbacks are unchanged in scale. Unconstrained and zero-objective solves take the exact previous path.

Alternatives considered:

Retry groundstate from a fresh projected random state on collapse — does not remove the attractor; with the observed ~60–70% per-attempt failure rate, bounded retries still fail often.
Custom DMRG loop (the issue's suggestion) — a much larger redesign; the shift achieves the goal inside ITensorMPS.dmrg.
Disabling noise for constrained solves — truncation also leaks, and noise helps convergence.

Evidence

Red on unchanged main (4fa351a): the new regression fixtures crash 6/8 and ~5/10 runs with zero feasible amplitude.
Green with the fix: 0 crashes and 0 wrong optima in 30/30 runs of each fixture (well beyond the 8+4 iterations the committed test performs), so the test is not expected to flake.
Full Pkg.test passes locally on Julia 1.12.0.
The issue's original per-variable fixture: applying this same shift (with the per-variable bound M'·|Q|·M + |l|'·M) locally on top of Per-variable domains #131's head takes it from 7/10 crashes to 0/10, confirming the mechanism transfers. That patch is validation only — nothing was pushed to Per-variable domains #131.

Branch Hygiene

Base: main (default), branched from origin/main at 4fa351a.
Overlap with Per-variable domains #131 (per-variable domains): both PRs touch minimize_mpo in src/backends/dmrg.jl. Whichever lands second rebases; the adaptation for per-variable domains is confined to objective_box_bound (per-variable maxima instead of a scalar maximum(abs, domain)), verified as described above.
The regression test lives in test/constrained_solve.jl (topically where constrained-solve tests live, and untouched by any open PR) rather than a new file, so no test/runtests.jl include-line conflict is introduced.

Follow-up

After #131 lands, add the issue's per-variable fixture (dom = [[0,2],[-1,0],[1,3,4]], RelationConstraint(1, :(<=), 3)) as a second regression case — it is not expressible on main today.

bernalde

Maintainer review of c8f0f1f against main (4fa351a). The changed-file set matches GitHub's (2 files). I re-derived the mechanism and re-ran the evidence rather than trusting the PR narrative.

Mechanism check. The diagnosis is correct: project_hamiltonian gives the infeasible kernel energy 0, so any fixture whose feasible spectrum is entirely positive makes the kernel the projected operator's ground space, and the post-sweep re-projection then annihilates the collapsed state. The fix places the whole in-domain spectrum strictly below 0 (λ = box bound + 1), which removes the attractor rather than papering over the crash; reported energies, per-iteration stats, logs, and callbacks all restore the shift, and the unconstrained, zero-objective, and infeasible paths are bit-for-bit the old ones (shift = 0 short-circuits to the previous expressions).

Evidence re-run at this head.

Regression fixtures red on unchanged main: 6/8 and ~5/10 crash rates reproduced.
With the fix: 0 crashes / 0 wrong optima in 30 runs per fixture (the committed test does 8+4), so flake risk is low.
Full Pkg.test on Julia 1.12.0: all suites pass, Constrained solving 147/147.
Existing behavior preserved: constrained tests with negative spectra, degenerate-optima amplitude checks, maximize forwarding, and the n=1 deterministic fast path all pass; the n=1 enumeration and the shifted spectrum commute (uniform shift preserves the argmin).
Integer-typed objectives: minimize(::Matrix{Int}, ...) with constraints fails identically at base and head (pre-existing eps(::Type{Int64})/dispatch errors upstream of this change), so real(T)(objective_bound) narrows nothing reachable.
Issue #132's own per-variable fixture is not expressible on main; the PR's claim that the adapted shift fixes it on top of #131 (7/10 → 0/10) matches my independent run of the same experiment.

Linked issue. Closes #132 is appropriate: the reported defect — checkably feasible models erroring during constrained solves — is fixed with a red-proven regression test, and the follow-up (committing the per-variable fixture once #131 lands) is recorded in the body and below.

Coordination. #131 also rewrites minimize_mpo; the Branch Hygiene section documents the overlap and the one-function adaptation, which I verified is confined to objective_box_bound. Whichever PR lands second rebases.

Nonblocking: shift magnitude on large dense problems

objective_box_bound sums |Q| entrywise, so λ grows like O(n²·max|Q|·max|domain|²) while individual energy gaps stay O(1). Krylov eigensolvers converge on relative gaps, so a large λ compresses the spectrum relative to the operator norm and may slow local convergence on large dense instances. This is safe (correctness is unaffected; today's problem sizes are far from the regime) and the bound is deliberately cheap, so no change is requested now — but if constrained DMRG convergence degrades on big models, a tighter bound (e.g. an estimate of max q − min q instead of the entrywise sum) is the first knob to try. Worth keeping in mind rather than in code.

Review summary

Blocking (0). Nonblocking (1: 0 inline, 1 body) — the shift-magnitude scaling note above, explicitly not requesting changes. Questions (0).

PR-head verification: full Pkg.test pass on Julia 1.12.0 at c8f0f1f in a clean worktree; red-on-base and 30-run stability probes as listed; git diff --check clean.

Merge-result verification: the branch is level with main (4fa351a); the merge is clean.

Current-head checks: all 11 green (nine test-matrix jobs, build, Documentation; docs preview deployed).

This review is a COMMENT because the PR author and reviewer are the same account; the human review this PR is being prepared for is the real gate.

iagoleal

This is a good start but I do not consider it enough to close #132.
The theoretical problem is still there, because DMRG may still pass through the kernel (it's a variational method after all!) and break for a feasible model.

We can merge this but should leave #132 open or create a followup until the technical problem itself is solved. There is a structural problem with DMRG + constraints as projections that this methodology is not enough to address.

Also, this shift may add huge values to relatively moderate problems, so I worry about numerical stability. Do we have any checks on larger models? I remember that we used to do a similar shift on the initial versions of this solver to increase the spectral gap, but it ended up being more of a problem than an aid during the solves.

iagoleal · 2026-08-20T20:20:00Z

+# this bound so the projection kernel cannot undercut the feasible spectrum;
+# see the shift construction in `minimize_mpo` and issue #132.
+function objective_box_bound(Q::AbstractMatrix, l::AbstractVector, domain)
+  M = maximum(abs, domain)


I am almost sure this fails for a non-uniform domain. No?

You were right — after #131 merged, maximum(abs, domain) over a Domains throws outright (iteration yields per-variable vectors). Fixed in 82379e3: the bound takes each variable's own maximum(abs, d). The #132 model in the testset exercises exactly this path.

iagoleal · 2026-08-20T20:21:21Z

+function objective_box_bound(Q::AbstractMatrix, l::AbstractVector, domain)
+  M = maximum(abs, domain)
+  return M^2 * sum(abs, Q) + M * sum(abs, l)
+end
+
+function objective_box_bound(p::AbstractPolynomial{T}, domain) where T
+  M = maximum(abs, domain)
+  bound = zero(abs(one(T)) * M)
+  for t in terms(p)
+    isconstant(t) && continue
+    bound += abs(coefficient(t)) * M^sum(last, powers(t))
+  end
+  return bound
+end


Why are we using the overall maximum instead of the per-domain maximum?
This would make a huge difference if the domains for different variables have different magnitudes, e.g.,
domain = [[0, 1], [0, 1e9]].

Done in 82379e3 — per-variable maxima (M' |Q| M + |l|' M, and per-term products for polynomials). On a [[0, 1], [0, 1000]] model this gives λ ≈ 1.0e3 where the global-max bound gave ≈ 1.0e6, and the model solves correctly 10/10.

iagoleal · 2026-08-20T20:22:10Z

                     , cutoff      = 1e-8  #  a cutoff of 1E-5 gives sensible accuracy; a cutoff of 1E-8 is high accuracy; and a cutoff of 1E-12 is near exact accuracy. (https://itensor.org/docs.cgi?page=tutorials/dmrg_params)
                     , verbosity   = 1
                     , constraints = AbstractConstraint[]
+                     , objective_bound :: Union{Nothing, Real} = nothing


All callers provide a bound and we own all callers. Let's not add unused defaults.

Suggested change

, objective_bound :: Union{Nothing, Real} = nothing

, objective_bound :: T

Applied as suggested: objective_bound :: T, required, no default. Both callers pass it explicitly.

iagoleal · 2026-08-20T20:23:33Z

+  # The projected Hamiltonian P'HP assigns energy zero to the infeasible
+  # subspace (the kernel of the projections). When every feasible objective
+  # value is positive, that kernel is the ground space, so the DMRG sweep is
+  # attracted into it and the solve collapses with zero feasible amplitude
+  # (issue #132). Shifting the objective spectrum below zero by more than its
+  # magnitude bound makes the feasible minimum the true ground state again;
+  # the shift is added back to every reported energy.


State only why, not what you are doing. The code should speak for itself.
Also, it's better to add full github links instead of only issue numbers.

Suggested change

# The projected Hamiltonian P'HP assigns energy zero to the infeasible

# subspace (the kernel of the projections). When every feasible objective

# value is positive, that kernel is the ground space, so the DMRG sweep is

# attracted into it and the solve collapses with zero feasible amplitude

# (issue #132). Shifting the objective spectrum below zero by more than its

# magnitude bound makes the feasible minimum the true ground state again;

# the shift is added back to every reported energy.

# The projected Hamiltonian P'HP assigns energy zero to the infeasible

# subspace (the kernel of the projections). When every feasible objective

# value is positive, that kernel is the ground space, so the DMRG sweep is

# attracted into it and the solve collapses with zero feasible amplitude

# ([issue #132](https://github.com/SECQUOIA/TenSolver.jl/issues/132)). Shifting the objective spectrum below zero by more than its

# magnitude bound makes the feasible minimum the true ground state again;

Applied your suggested text verbatim (only closing the trailing semicolon as a period, since the clause after it was the part being dropped).

iagoleal · 2026-08-20T20:27:32Z

+  shift = if zero_objective || isempty(projections) || isnothing(objective_bound)
+    zero(real(T))
+  else
+    real(T)(objective_bound) + one(real(T))
+  end
+  H_solve = iszero(shift) ? H_obj : H_obj - shift * ITensorMPS.MPO(T, sites, "Id")


I am not sure this logic should be in here. Why don't we add this term to the Hamiltonian construction itself?

Also, the objective_bound should add the one itself... why is it adding one in here?

Both parts in 82379e3: objective_box_bound folds the one-unit margin in itself, and the shift + H_solve construction now live inside the Hamiltonian-construction block. I stopped short of pushing the term into tensorize: it is constraint-agnostic and shared with unconstrained solves, so the shift there would either apply where there is no kernel to avoid (extra bond dimension for nothing) or force tensorize to know about constraints. Happy to move it anyway if you prefer that layering.

iagoleal · 2026-08-20T20:29:06Z

                      , eigsolve_verbosity = 0
                 )
+    # The solve ran on the shifted spectrum; report objective-scale energies.
+    energy += shift


Add it together with the constant at report time. Else, we will have a lot of scattered changes to the energy value and it can get confusing. Prefer to centralize these things.

Done — one objective = energy + shift + c per iteration feeds the stats, the log line, and the callback; the raw DMRG energy is no longer mutated after groundstate.

iagoleal · 2026-08-20T20:30:00Z

+    # subspace. That subspace is the kernel of the projections. With the
+    # spectral shift applied by `minimize_mpo` the kernel sits above every
+    # feasible energy, so the sweep itself penalizes leaked amplitude back out
+    # and this re-projection is only a numerical cleanup.


How are you sure of this?

I wasn't — fair catch, and your review body makes the same point. The comment now says the shift makes the kernel energetically unfavorable, which suppresses leakage but cannot forbid it (local updates, deliberate noise), so the re-projection remains the feasibility guarantee. The structural question is tracked in #134.

iagoleal · 2026-08-20T20:30:53Z

+  @testset "Feasible model with all-positive spectrum does not collapse" begin
+    # Regression for #132: the projected Hamiltonian assigns the infeasible
+    # kernel energy zero, so when every feasible objective value is positive
+    # the kernel is the DMRG attractor and the solve used to die with
+    # "zero feasible amplitude" on most runs. The spectral shift in
+    # minimize_mpo keeps the feasible minimum below the kernel.
+    Q = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
+
+    # x1 >= x3 over {1, 2}: feasible minimum 3 at [1, 1, 1].
+    rel = AbstractConstraint[RelationConstraint(1, :(>=), 3)]
+    for _ in 1:8
+      E, psi = minimize(Q; domain = [1, 2], constraints = rel, verbosity = 0)
+      @test E ≈ 3.0
+      @test is_feasible(TenSolver.sample(psi), rel)
+    end
+
+    # Forbidding the unconstrained optimum [1, 1, 1] over {1, 2, 3}:
+    # feasible minimum 6 at permutations of [2, 1, 1].
+    noteq = AbstractConstraint[NotEqualsConstraint([1, 2, 3], [1, 1, 1])]
+    for _ in 1:4


The model from the issue is mussing here. That is the one we should be testing...

dom = [[0.0, 2.0], [-1.0, 0.0], [1.0, 3.0, 4.0]] Q = [1.0 0.3 0.0; 0.0 -1.0 0.2; 0.0 0.0 1.0] l = [0.1, -0.2, 0.3] minimize(Q, l; domain = dom, constraints = [RelationConstraint(1, :(<=), 3)], verbosity = 0)

Added in 82379e3 — expressible now that #131 is on main. Red without the fix at current main (6/10 crashes on that exact model), green with it (0 failures in 30 runs; the committed test performs 4).

bernalde · 2026-08-21T00:56:08Z

Addressed the review at 82379e3, which is a merge of current main (per-variable domains, #131) with the shift re-applied and revised per your comments.

Commits pushed: 82379e3 (merge + revisions; the branch was DIRTY against main after #131 landed, so the conflict resolution and the review revisions are one commit).

Main changes:

objective_box_bound now uses each variable's own domain magnitude (M' |Q| M + |l|' M and the per-term polynomial analogue) instead of the overall maximum, and folds the strictness margin in itself. For mixed-magnitude domains this is a large tightening: on a [[0, 1], [0, 1000]] model, λ ≈ 1.0e3 versus ≈ 1.0e6 under the old global-max bound, and the model solves correctly 10/10.
minimize_mpo takes objective_bound :: T as a required keyword (your suggestion applied); the shift and H_solve construction live inside the Hamiltonian-construction block; every reported energy is assembled once per iteration as energy + shift + c.
Your comment-text suggestion applied verbatim; the re-projection comment no longer overclaims (details in the thread replies).
The regression testset now includes the exact model from Infeasibility reached for feasible model on DMRG solve #132 — expressible now that Per-variable domains #131 is on main — alongside the two uniform-domain fixtures.

On "not enough to close #132": agreed that the structural problem remains — DMRG is variational and can still traverse the kernel; the shift removes the systematic attractor, not the possibility. Per your "leave #132 open or create a followup", I created #134 ("Structural robustness of constrained DMRG: sweeps can still traverse the projection kernel"), which records the mechanisms discussed (feasible-manifold DMRG loop, ITE via #93, in-sweep re-projection, graceful restart) and acceptance criteria. With #134 owning the structural work, this PR's Closes #132 closes the concrete collapse bug. If you would rather keep #132 itself open, I have not edited the PR body — flipping Closes to Refs is a one-word edit on your side or mine, just say which you prefer.

On numerical stability of the shift (your review-body question): measured rather than argued. On a random dense n = 15 QUBO with a cardinality constraint (λ ≈ 109, brute-force-verifiable), five shifted solves and five unshifted solves (same model, objective_bound = 0 internally) land on the same feasible local optimum with gap 0.525 — the shifted runs found the exact optimum once, the unshifted never did. So at this scale the shift is at parity or slightly better, and the optimality gap it does show is the pre-existing projected-DMRG local-minimum behavior (#19), reproduced identically without the shift. The per-variable bound above is the main guard against the "huge values on moderate problems" scenario, since λ now scales with each variable's own range. I did not benchmark beyond n = 15; if you remember which model family made the early-days spectral-gap shift harmful, that is worth adding to #134's validation matrix.

Tests: full Pkg.test green locally on Julia 1.12.0 at 82379e3 (Constrained solving 155/155, zero failures suite-wide). All three regression fixture families: 0 crashes / 0 wrong optima in 30 runs each; on current main without the fix they crash 5/10, 7/10, and 6/10 (the last being the #132 model).

Not addressed: nothing declined. The formal review decision remains CHANGES_REQUESTED until you re-review or dismiss.

bernalde

Re-review at 82379e3, following my review of c8f0f1f and @iagoleal's CHANGES_REQUESTED round. The head is the merge of current main (post-#131) with the shift re-applied; the merge-result diff against main is exactly the intended two files, and the changed-file set matches GitHub's.

Prior findings, verified at this head

All ten items from the CHANGES_REQUESTED round are addressed in the diff itself:

Per-variable bound: objective_box_bound takes each variable's own maximum(abs, d); the old overall-maximum form would in fact throw on the merged Domains type, so the non-uniform concern was even sharper than stated. Measured effect on [[0, 1], [0, 1000]]: λ ≈ 1.0e3 versus ≈ 1.0e6, solving correctly 10/10.
objective_bound :: T is a required keyword with no default; both owned callers pass it.
The strictness margin lives inside the bound helper; the shift and H_solve are built in the Hamiltonian-construction block; reported energies are assembled once per iteration as energy + shift + c (stats, log line, and callback all consume the same value).
The reviewer-prescribed comment text is applied verbatim, and the re-projection comment now states only what is defensible: the shift suppresses kernel leakage but cannot forbid it, so re-projection remains the feasibility guarantee. The structural limitation is tracked in #134, cross-linked from the PR discussion.
The exact #132 model is in the regression testset alongside the two uniform-domain fixtures.

My earlier nonblocking note (shift magnitude on large dense problems) is superseded by the per-variable bound plus the measured n = 15 parity result: shifted and unshifted constrained solves land on the same feasible local optimum (the shifted runs found the exact optimum once in five; the unshifted never did), so the residual optimality gap is the pre-existing projected-DMRG local-minimum behavior (#19), not shift-induced.

Findings

None — Blocking (0), Nonblocking (0), Questions (0).

PR-head verification: full Pkg.test on Julia 1.12.0 at exactly 82379e3: all suites pass, Constrained solving 155/155, zero failures. Regression fixtures: 0 crashes / 0 wrong optima in 30 runs per family, against 5/10, 7/10, and 6/10 crash rates on current main without the fix. git diff --check clean.

Merge-result verification: the head already contains the merge of main (5b49f2b); the diff against main is the two intended files and mergeStateStatus is CLEAN.

Current-head checks: all 11 green (nine test-matrix jobs, build, Documentation with deployed preview).

Coordination: #44 also touches src/backends/dmrg.jl (import lines, based on pre-#131 main); whichever lands second rebases, and #44 needs a rebase against current main regardless.

This PR carries @iagoleal's approval at exactly this head and no gate remains other than the merge decision itself. This review is a COMMENT because the PR author and reviewer are the same account.

The projected Hamiltonian P'HP assigns energy zero to the infeasible subspace, so a feasible model whose objective is positive over the whole feasible set makes that kernel the DMRG ground space: sweeps collapse into it and the solve dies with "zero feasible amplitude" (#132). Constrained solves now shift the objective spectrum below zero by more than a cheap box bound on its magnitude, making the feasible minimum the true ground state again; the shift is added back to every reported energy. Unconstrained and zero-objective solves are unchanged.

Re-applies the #132 spectral shift on top of the per-variable Domains API that landed with #131, revised per review: - objective_box_bound uses each variable's own domain magnitude instead of the overall maximum, tightening the shift by orders of magnitude for mixed-magnitude domains, and folds the strictness margin in itself. - minimize_mpo takes objective_bound :: T as a required keyword; no internal default. - The shift is built inside the Hamiltonian-construction block, and every reported energy is assembled once per iteration as energy + shift + c. - The re-projection comment no longer overclaims: the shift suppresses kernel leakage but re-projection remains the feasibility guarantee. - The regression testset gains the exact model from #132, expressible now that per-variable domains are on main.

bernalde marked this pull request as ready for review August 20, 2026 18:41

bernalde commented Aug 20, 2026

View reviewed changes

iagoleal requested changes Aug 20, 2026

View reviewed changes

bernalde mentioned this pull request Aug 21, 2026

Structural robustness of constrained DMRG: sweeps can still traverse the projection kernel #134

Open

2 tasks

bernalde requested a review from iagoleal August 21, 2026 00:57

iagoleal approved these changes Aug 21, 2026

View reviewed changes

iagoleal force-pushed the fix/issue-132-feasible-collapse branch 2 times, most recently from e82846b to 82379e3 Compare August 21, 2026 18:00

bernalde commented Aug 21, 2026

View reviewed changes

bernalde added 2 commits August 22, 2026 10:28

iagoleal force-pushed the fix/issue-132-feasible-collapse branch from 82379e3 to 531c6b3 Compare August 22, 2026 13:36

iagoleal merged commit af48e6b into main Aug 22, 2026
11 checks passed

iagoleal deleted the fix/issue-132-feasible-collapse branch August 22, 2026 13:48

	, objective_bound :: Union{Nothing, Real} = nothing
	, objective_bound :: T

Conversation

bernalde commented Aug 20, 2026

Mechanism

Fix

Evidence

Branch Hygiene

Follow-up

Uh oh!

bernalde left a comment

Choose a reason for hiding this comment

Nonblocking: shift magnitude on large dense problems

Review summary

Uh oh!

iagoleal left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

bernalde commented Aug 21, 2026

Uh oh!

bernalde left a comment

Choose a reason for hiding this comment

Prior findings, verified at this head

Findings

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

iagoleal left a comment •

edited

Loading