From d075f600271f9aaaa0d00edc4c2fedfa32cba7f5 Mon Sep 17 00:00:00 2001
From: Miles <miles.stoudenmire@gmail.com>
Date: Sun, 14 Dec 2025 21:54:11 -0500
Subject: [PATCH] Simplify and fix issue with densitymatrix subspace expansion
 by having it take truncation parameters from factorize_kwargs.

---
 Project.toml                          | 2 +-
 src/solvers/eigsolve.jl               | 1 -
 src/solvers/subspace/densitymatrix.jl | 3 +--
 test/solvers/test_eigsolve.jl         | 8 ++++++--
 4 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/Project.toml b/Project.toml
index 2ce216a7..fad02f9f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ITensorNetworks"
 uuid = "2919e153-833c-4bdc-8836-1ea460a35fc7"
 authors = ["Matthew Fishman <mfishman@flatironinstitute.org>, Joseph Tindall <jtindall@flatironinstitute.org> and contributors"]
-version = "0.15.5"
+version = "0.15.6"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl
index a922d5a8..23a42d25 100644
--- a/src/solvers/eigsolve.jl
+++ b/src/solvers/eigsolve.jl
@@ -67,7 +67,6 @@ function eigsolve(
         nsweeps;
         nsites,
         factorize_kwargs,
-        subspace_expand!_kwargs = (; eigen_kwargs = factorize_kwargs),
         sweep_kwargs...,
     )
     prob = problem(sweep_solve!(sweep_iter))
diff --git a/src/solvers/subspace/densitymatrix.jl b/src/solvers/subspace/densitymatrix.jl
index ba1b39df..50a899c0 100644
--- a/src/solvers/subspace/densitymatrix.jl
+++ b/src/solvers/subspace/densitymatrix.jl
@@ -8,7 +8,6 @@ function subspace_expand!(
         expansion_factor = 1.5,
         maxexpand = typemax(Int),
         north_pass = 1,
-        eigen_kwargs = (;),
     )
     prob = problem(region_iter)
 
@@ -29,7 +28,7 @@ function subspace_expand!(
     isnothing(a) && return region_iter, local_state
     basis_size = prod(dim.(uniqueinds(A, C)))
 
-    trunc_kwargs = truncation_parameters(region_iter.which_sweep; eigen_kwargs...)
+    trunc_kwargs = truncation_parameters(region_iter.which_sweep; region_kwargs(factorize, region_iter)...)
     expanded_maxdim = compute_expansion(
         dim(a), basis_size; expansion_factor, maxexpand, trunc_kwargs.maxdim
     )
diff --git a/test/solvers/test_eigsolve.jl b/test/solvers/test_eigsolve.jl
index df7c47de..f69ca8f6 100644
--- a/test/solvers/test_eigsolve.jl
+++ b/test/solvers/test_eigsolve.jl
@@ -1,6 +1,6 @@
 using Test: @test, @testset
 using ITensors
-using ITensorNetworks: siteinds, ttn, dmrg
+using ITensorNetworks: dmrg, maxlinkdim, siteinds, ttn
 using Graphs: dst, edges, src, vertices
 using ITensorMPS: OpSum
 using TensorOperations: TensorOperations #For contraction order finding
@@ -58,12 +58,16 @@ include("utilities/tree_graphs.jl")
     #
     nsites = 1
     nsweeps = 5
-    factorize_kwargs = (; cutoff = [1.0e-5, 1.0e-6], maxdim = [8, 16, 32])
+    maxdim = [8, 16, 32]
+    factorize_kwargs = (; cutoff = [1.0e-5, 1.0e-6], maxdim)
     extract!_kwargs = (; subspace_algorithm = "densitymatrix")
     E, psi = dmrg(H, psi0; extract!_kwargs, factorize_kwargs, nsites, nsweeps, outputlevel)
     (outputlevel >= 1) && println("1-site+subspace DMRG energy = ", E)
     @test E ≈ Ex atol = 1.0e-5
 
+    # Regression test that subspace expansion feature obeys maxdim limit
+    @test maxlinkdim(psi) <= last(maxdim)
+
     #
     # Test passing cutoff and maxdim as a vector of values
     #