Conversation
The `MPO_AC2_Hamiltonian{<:MPSBondTensor, Nothing, Nothing, <:MPSBondTensor}`
action (the operator-free / overlap two-site derivative) was missing a
semicolon in the `rightenv` index, declaring a `(2, 0)` index pattern for a
`(1, 1)` bond tensor. This threw a `TensorOperations.IndexError` for the planar
bond-tensor case (e.g. `approximate(x, b, DMRG2())` on planar states); the
ComplexSpace path happened to avoid the broken method.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Use the fully-qualified `KrylovKit.linsolve` at the internal transfer-matrix / quasiparticle / fidelity call sites, so that the top-level `linsolve` name can be owned by MPSKit for the upcoming MPS linear-solver category without shadowing these tensor-level solves. Behavior-preserving. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Add a `truncation_factor` field to `DynamicTol` and let `adapt_solver` floor the retuned tolerance at `truncation_factor * eps_trunc`, i.e. `tol = clamp(max(truncation_factor*eps_trunc, tol_factor*g_global/sqrt(iter)), tol_min, tol_max)`. The default `truncation_factor = 0` keeps existing behavior (DMRG's gauge / eigensolver uses are unchanged); a positive value is used by truncating two-site sweeps to avoid solving a bond far below what the SVD truncation will discard. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Add a first-class `linsolve`/`linsolve!` for finite MPS, solving `(a0 + a1*A)*x = b` with `A` an MPO/MPOHamiltonian and `x`, `b` MPS. The signature is state-first (`linsolve(x0, A, b, [alg], [envs]; a0=0, a1=1, ...)`) and returns `(x, envs, eps)`, matching find_groundstate/approximate; the shift follows KrylovKit's convention. Algorithms are DMRG-style sweeps `DMRGSolve` (single-site) and `DMRGSolve2` (two-site, rank-adaptive), each parameterized by a local `LinsolveFormulation`: `Galerkin` (solve the effective system directly; CG/GMRES/BiCGStab) and `LeastSquares` (normal equations via squared environments). Convergence uses the local Galerkin residual `||(a0 + a1*A_eff)*AC - b_eff|| / ||b||`, formed as a local tensor (machine-precision floor, no catastrophic cancellation). The two-site stop is truncation-aware. Local solves default to adaptive tolerances via `Defaults.alg_linsolve()` (a structure-selected CG/MINRES/GMRES wrapped in `DynamicTol`, with a truncation floor for two-site); pass a plain solver to disable adaptation. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Cover the resolvent `(z - H)x = V|gs>` against the analytic `1/(z - E0)` and the existing `propagator` (both Galerkin and LeastSquares, single- and two-site), a positive-definite solve via the keyword interface (CG path), adaptive tolerances including BiCGStab, and a general RHS checked against a dense solve. Residuals are verified densely (MPS-level +/- does not reliably cancel, cf. #473). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Three fixes to the local solves: - `Defaults.alg_linsolve(; ishermitian = true)` returned a `MINRES` solver, but KrylovKit exports the type without implementing `linsolve` for it (see its own "TODO: implement MINRES for symmetric but not posdef; for now use GMRES"), so the hermitian-indefinite path threw a `MethodError`. Fall back to GMRES too. - `@doc str a, b` attaches only to the last binding, so exported `linsolve` had no docstring at all. Share it explicitly with `@doc (@doc linsolve) linsolve!`. - A local solve no longer warns about non-convergence when it undershot its own adaptive tolerance but still reached `tol * norm(b)`, the accuracy the outer sweep needs. Adaptive tolerances routinely dip below the round-off floor of the local problem, which made the warning fire on nearly every tightly converged sweep. `LeastSquares` now forms the `A†b` term of the normal equations from the mixed sandwich <x|A|b> instead of materializing `A*b` as an MPS, avoiding a state of bond dimension chi_A*chi_b together with its environments. The JordanMPO effective operator assumes `below === above` and branching on that inside `AC_hamiltonian` would make the return type of that hot function a union, so the mixed case goes through the generic sparse MPO derivative in `_mixed_projection` -- one application per bond, so its cost is irrelevant. This also makes the formulation usable on `WindowMPS`, which the `A*b` route never could be: there is no `*(::WindowMPOHamiltonian, ::WindowMPS)`. Also adds the missing in-place `linsolve!(x0, A, b, ::UnionAlg, envs...)`, the counterpart of the existing copying method. Tests cover two-site `LeastSquares` and a `WindowMPS` resolvent for both formulations. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`corvector.jl` carried its own DMRG sweep, which had drifted from the one in the `linsolve` category on three counts: it converged on the per-site update norm ||AC' - AC|| (which cannot tell converged from stalled) rather than on the residual of the linear system, it used a fixed-tolerance GMRES rather than the structure-selected adaptive solver, and it was single-site only. `propagator` is now a single method for both flavours whose body is a `linsolve!` call. Both flavours keep their variational problem unchanged. `NaiveInvert` is a `Galerkin` solve of `(z - H)x = |psi0>`. `Jeckelmann` still minimizes functional (14) of Jeckelmann2002 and reconstructs G(z) through equation (11): its operator is assembled as `LinearCombination((H^2, H), (1, -2w))` over a hand-built `LazyLincoCache` -- hand-built so that `squaredenvs`' boundary conditions, which matter for `WindowMPS`, survive -- and its stationarity condition `((w - H)^2 + eta^2) psi = -eta |psi0>` is scaled by `-1/eta` so that the right-hand side is `|psi0>` itself. No scaled copy of the right-hand side is needed and `x` is the identical vector; the shift comes out real, so a real initial guess stays real. Consequent changes to the defaults of `DynamicalDMRG`: - `tol` defaults to `Defaults.tol` rather than `10 * Defaults.tol`, and now measures the relative residual of the linear system. - `solver` defaults to `Defaults.alg_linsolve()` (adaptive tolerances), chosen per flavour through an explicit keyword constructor. - a new `trscheme` field prepends a truncated two-site sweep, so the bond dimension no longer has to be fixed by the initial guess. - iteration logs are labelled `linsolve` instead of `DDMRG`. - the `Jeckelmann` flavour throws for real `z` instead of returning `NaN`; the equation-(11) reconstruction divides by `imag(z)`. The docstrings claimed that the second return value is `1/(z - H)|psi0>` for both flavours. That is only true for `NaiveInvert`; for `Jeckelmann` it is the minimizer of the functional, i.e. the imaginary part of the propagator. Say so. Adds a `linsolve` section to the algorithms manual: the category was exported but appeared nowhere in the docs. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Member
|
It might be good (or not good because to broad) to consider the names of all the high level methods, and try to make them more consistent and streamlined. This would of course be massively breaking, but on the other hand, renaming a method is probably the easiest type of breaking change in terms of deprecations, and the old name can stay around for a very very long time. So what do we currently have:
and maybe
The current naming is certainly all over the place 😄 . |
Member
|
That being said, I do think the current situation is already much nicer than some of the other tensor network packages, which often use the specific algorithm as a method name, and certainly also have a wild variety of names around. |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This PR adds dedicated routines for
linsolve, i.e. solvingAx = bwith MPO/MPS.This was already partially implemented in the
propagatorand dynamical DMRG stuff, but the goal is to make it a bit more of a first-class citizen and not hardcode the specific linear problems to solve.Can also be used for time-stepping purposes, doing things like trapezoid rules.
There definitely are some design choices that I'm happy to have discussions about, considering names etc, but one particular choice I made here is to not overload KrylovKit's version.
The reasoning behind this is that the
envswere a bit awkward to get in there anyways, and KrylovKit does not have an in-place version to begin with. It does have the annoying side-effect of nameclashes, since I definitely want to export this, so maybe a different name could be useful?