API Reference

This page documents the public API of TenSolver.jl.

Optimization Functions

TenSolver.minimizeFunction
minimize([Q::Matrix], [l::Vector], [c::Number] ; domain, kwargs...)
minimize(p::AbstractPolynomial ; domain, kwargs...)

Solve a polynomial discrete optimization problem

min  p(x)
s.t. x_i in domain
     constraints

In the matrix version, the objective is limited to quadratic forms x -> x'Qx + l'x + c. Missing arguments (quadratic, linear or constant term) are allowed and taken to be zero.

Return the optimal value E and a probability distribution ψ over optimal solutions. You can use sample to get an actual solution vector from ψ.

There are multiple backends available, selected through the keyword backend. By default, it uses DMRG to calculate the optimal solution.

Keyword arguments:

  • constraints :: AbstractVector{<:AbstractConstraint} - Experimental native Julia hard constraints. Defaults to AbstractConstraint[]. In constrained DMRG solves, TenSolver lowers each constraint to a projection MPO, solves the projected Hamiltonian, and returns a feasible sampled assignment. For polynomial objectives, constraints are expressed in the same order as their effective_variables. If the constraints admit no solution at all, the solve does not error: it logs a warning and returns +Inf together with an infeasible Solution (see is_feasible).

  • domain - Allowed values per variable. Can be either a vector of per-variable domains or a single uniform domain valid for all variables. Defaults to [0, 1]. Unconstrained DMRG optimization accepts any finite collection of real values; individual constraint types can impose narrower requirements. Use [-1, 1] for Ising spins. Domains are sorted and deduplicated before solving.

  • iterations :: Int - Maximum iterations the solver should run. Defaults to 10.

  • cutoff :: Float64 - Any absolute value below this threshold is considered zero. Defaults to 1e-8. You can use this keyword to control the solver's accuracy vs resources trade-off.

  • time_limit :: Float64 - If specified, determines the maximum running time in seconds. It only determines whether a new iteration should start or not, thus the solver may run for longer if the threshold happens during an iteration.

  • device = cpu - Accelerator device used during computation. See the section below for how to run on GPUs.

  • 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. Use to collect statistics or serialize intermediate states. psi is the MPS for that iteration. Default: nothing (no callback).

  • callback_every :: Int - Invoke the callback every N iterations. Must be >= 1. Default: 1.

  • backend - Solver backend. Defaults to the current DMRG implementation. Use backend = :dmrg or backend = DMRGBackend() to select it explicitly. Other backends are reserved for optional extensions.

    Other keywords might be available depending on the chosen backend. See the documentation for each backend for comprehensive lists.

    Some keywords, such as constraints and domain, may have limited support depending on the backend.

The returned Solution carries per-iteration convergence data in solution.stats.

Provably infeasible constrained models are reported as a status: minimize logs a warning and returns +Inf (the minimum over an empty feasible set) together with an infeasible Solution, which cannot be sampled. Check it with is_feasible.

Running on GPU:

The optional keyword device controls whether the solver should run on CPU or GPU. For using a GPU, you can import the respective package, e.g. CUDA.jl, and pass its accelerator as argument.

import CUDA
minimize(Q; device = CUDA.cu)

import Metal
minimize(Q; device = Metal.mtl)

See also maximize.

source
TenSolver.maximizeFunction
maximize([Q::Matrix], [l::Vector], [c::Number] ; domain, kwargs...)
maximize(p::AbstractPolynomial ; domain, kwargs...)

Solve a polynomial discrete optimization problem

max  p(x)
s.t. x_i in domain
     constraints

In the matrix version, the objective is limited to quadratic forms x -> x'Qx + l'x + c. Missing arguments (quadratic, linear or constant term) are allowed and taken to be zero.

All keywords accepted by minimize can also be used for maximization problems. Provably infeasible constrained models return -Inf (the supremum over an empty feasible set) together with an infeasible Solution.

See also minimize.

source

Solver Backends

TenSolver.AbstractTenSolverBackendType
AbstractTenSolverBackend

Abstract solver backend marker for TenSolver implementations.

Backends must provide backend-specific minimize methods for the normalized optimization inputs they support. Matrix backends implement minimize(::MyBackend, Q::AbstractMatrix, l, c; kwargs...); polynomial backends implement minimize(::MyBackend, p::AbstractPolynomial; kwargs...). Extensions that support symbolic selection must also define normalize_backend(::Val{:my_backend}) = MyBackend(...).

See also

DMRGBackend, normalize_backend

source
TenSolver.normalize_backendFunction
normalize_backend(backend)

Normalize a user-facing backend selector into a backend object.

Backends can support backend = :my_backend by defining normalize_backend(::Val{:my_backend}) = MyBackend(...).

source

Solution

TenSolver.SolutionType
Solution{T}

The result of running minimize or maximize: an MPS wave function over the optimal solution space, together with per-iteration convergence stats.

Use sample to draw vectors from it.

Fields

  • tensor: the underlying MPS, or nothing when the model is infeasible.
  • domain: possible variable values.
  • permutation: original variable index represented by each tensor site.
  • stats: per-iteration convergence stats. See SolverStatistics.

Provably infeasible models produce a Solution with no MPS and empty stats vectors; check with is_feasible before sampling.

source
TenSolver.SolverStatisticsType
SolverStatistics{T}

Fields

  • energies: per-iteration calculated energy/objective value;
  • bond_dims: per-iteration solution maximum bond dimension;
  • elapsed_times: per-iteration total until the iteration completed;
  • variances: per-iteration Hamiltonian variance, or nothing on iterations where the variance was not checked;
  • max_bonds: structure containing bond dimensions for multiple tensors used throughout the iteration
    • initial_state: bond for the initial MPS guess;
    • objective: bond for the MPO representing the objective function H;
    • projections: bonds for the MPOs representing each constraint;
    • hamiltonian: bond for the MPO representing the actual Hamiltonian P'HP representing both objective and constraints;
source

Sampling Functions

TenSolver.sampleFunction
sample(psi)

Sample a vector from a (quantum) probability distribution.

Throw a DomainError when psi is infeasible (see is_feasible), since there is no solution to query.

source

Boolean/Spin Conversions

TenSolver.bool_to_spinFunction
bool_to_spin(x)

Convert a Boolean bit vector x_i in {0, 1} to Ising spins s_i in {-1, +1} using QUBOTools' BoolDomain => SpinDomain cast.

source
TenSolver.spin_to_boolFunction
spin_to_bool(s)

Convert an Ising spin vector s_i in {-1, +1} to Boolean bits x_i in {0, 1} using QUBOTools' SpinDomain => BoolDomain cast.

source
TenSolver.qubo_to_isingFunction
qubo_to_ising(Q[, l[, c]]; convention = :spin)
qubo_to_ising(form; convention = :spin)

Convert TenSolver's Boolean QUBO objective

dot(x, Q, x) + dot(l, x) + c

with x_i in {0, 1} into a sparse QUBOTools form in SpinDomain. For non-symmetric matrices, the QUBOTools form constructor preserves TenSolver's dot(x, Q, x) convention by storing the effective pair coefficient Q[i, j] + Q[j, i] once in the upper triangle.

The Boolean-to-spin conversion introduces halves and quarters. Integer coefficient inputs therefore return floating-point forms, while rational inputs preserve exact rational arithmetic.

The returned form includes the constant offset, so

dot(x, Q, x) + dot(l, x) + c == QUBOTools.value(bool_to_spin(x), qubo_to_ising(Q, l, c))

for every Boolean vector x.

source
TenSolver.ising_to_quboFunction
ising_to_qubo(form)
ising_to_qubo(J, h[, offset])

Convert a QUBOTools spin-domain Ising form back to a sparse Boolean-domain QUBOTools form. The matrix/vector method first builds a QUBOTools spin form from

dot(s, J, s) + dot(h, s) + offset

with s_i in {-1, +1}. QUBOTools folds diagonal quadratic spin terms into the constant offset and stores each off-diagonal unordered pair once in the upper triangle.

As with qubo_to_ising, integer coefficient inputs return floating-point forms when the conversion introduces fractional coefficients, while rational inputs preserve exact rational arithmetic.

source

Constraints

Hard constraints are enforced by lowering each one to an exact projection MPO, following CoTenN (Sharma, Peng, Dangwal, and Achour, "CoTenN: Constrained Optimization with Tensor Networks," PLDI 2026). See Constrained Optimization for a worked example.

TenSolver.AbstractConstraintType
AbstractConstraint

Supertype for conditions over a vector x addressed by 1-based site indices.

API

These are an interface for feasibility constraints. Any concrete subtype is expected to implement

Constraint types are experimental. They currently provide TenSolver's Julia lowering target for projection-MPO constrained solves; future JuMP/MOI integration may change which constraint abstraction is considered stable public API.

See also SumConstraint, SumModConstraint, NotEqualsConstraint, AssignmentConstraint, and RelationConstraint.

source
TenSolver.SumConstraintType
SumConstraint{T} <: AbstractConstraint
SumConstraint(sites, weights, relation, rhs)
SumConstraint(sites, weights, rhs; relation)

Weighted-sum constraint over a vector x:

sum(weights[i] * x[sites[i]] for i in eachindex(sites)) relation rhs.

sites must be unique positive integers, weights must be the same length as sites and only contain nonnegative integers, and relation must be one of :(==), :(!=), :(<=), or :(>=).

Warning: The == and != relations use exact arithmetic comparison.

source
TenSolver.SumModConstraintType
SumModConstraint{T} <: AbstractConstraint
SumModConstraint(sites, weights, rhs; mod)

Modular weighted-sum constraint over a vector x:

sum(weights[i] * x[sites[i]] for i in eachindex(sites)) ≡ rhs (mod m).

sites must be unique positive integers, weights must be the same length as sites, weights and rhs must be integer-valued, and mod must be a positive integer.

Weights and the rhs are stored as their least nonnegative residues modulo mod.

source
TenSolver.NotEqualsConstraintType
NotEqualsConstraint <: AbstractConstraint
NotEqualsConstraint(sites, values)

Excludes a single assignment over a vector x: at least one component of x[sites] must differ from values. Equivalently, the partial assignment x[sites] == values is forbidden.

sites must be unique positive integers, and values must have the same length as sites.

source
TenSolver.AssignmentConstraintType
AssignmentConstraint{T} <: AbstractConstraint
AssignmentConstraint(sites, values, relation, rhs)

Restrict how many sites satisfy x[site] in values, i.e.,

count(x[site] in values for site in sites) relation rhs.

rhs must be a nonnegative integer. The count and rhs are stored independently of the numeric element type of values.

A common application is to restrict exactly one variable to be a certain value,

ExactlyOne(sites, value) = AssignmentConstraint(sites, [value], :(==), 1)
source
TenSolver.RelationConstraintType
RelationConstraint <: AbstractConstraint
RelationConstraint(left_site, relation, right_site)

Pairwise constraint over a vector x: x[left_site] relation x[right_site].

left_site and right_site must be distinct positive integers, and relation must be one of :(==), :(!=), :(<=), or :(>=).

source
TenSolver.is_feasibleFunction
is_feasible(x, constraint::AbstractConstraint)

Test whether the vector x satisfies a single constraint.

source

Utility Functions

Base.inMethod
in(xs, psi::Solution [; cutoff)

Whether the vector xs has a positive probability of being sampleable from psi. When setting cutoff, it will be used as the minimum probability considered positive. Always false for infeasible solutions.

source
TenSolver.permuteFunction
permute(c::AbstractContraint, p)

Reorder the constraint sites of a constraint according to a permutation p. This effectively converts the constraint to one with the same semantics but applied to an optimization model with reordered variables.

source

Internal Functions

These functions are part of the internal implementation and are not exported. They are documented here for advanced users who may need to understand the internals. Notice: As unexported method and types, they are subject to change without warning.

Objective Construction

TenSolver.tensorizeFunction
tensorize(p)

Turn a polynomial function action on bitstrings into an equivalent MPO Hamiltonian acting on Qudit sites. The conversion consists of exchanging each integer variable x_i for a matrix P_i whose eigenvalues represent its feasible set K_i.

∑ Q_ij x_i x_j + ∑ l_i x_i --> H = Σ Q_ij D_i D_j + ∑ l_i D_i
source
TenSolver.qmatrix_permutationFunction
qmatrix_permutation(Q; cutoff)

Return a deterministic permutation that places coupled QUBO variables closer together in the one-dimensional MPS ordering.

cutoff controls which quadratic couplings are included in the ordering graph: an undirected edge between variables i and j is used when abs(Q[i, j] + Q[j, i]) > cutoff. The default cutoff = 0 preserves every nonzero coupling for callers that only want the ordering.

The returned permutation maps each tensor site to its original variable index: entry k is the original QUBO variable represented by tensor site k.

Examples

Q = [0.0 0.0 1.0;
     0.0 0.0 0.0;
     1.0 0.0 0.0]

permutation = qmatrix_permutation(Q; cutoff=0)
Q[permutation, permutation]

# output

3×3 Matrix{Float64}:
 0.0  1.0  0.0
 1.0  0.0  0.0
 0.0  0.0  0.0
source
TenSolver.preprocess_modelFunction
preprocess_model(Q, l, c; domain, constraints, cutoff)

Permute variables before Hamiltonian construction so coupled variables are closer in the one-dimensional tensor order.

source

MPO Construction

TenSolver.DFAType
DFA{S, A}

Deterministic finite automaton with step-dependent and partial transitions.

Fields:

  • states: DFA states, used to define the MPO bond dimension.
  • alphabet: Per-stage DFA alphabet.
  • initial: start state.
  • accepting: set of accepting states.
  • transitions: one transition table per step; each table maps (state, symbol) to the next state. Missing entries are rejected.
source
TenSolver.mapreduce_dfaFunction
mapreduce_dfa(f, op, constraint, alphabets; initial, predicate, states)

Build a DFA by mapping each constrained site symbol through f and combining the result in a state accumulator with op.

The function f is assumed to take the states to a set where op acts as a monoid operation, i.e., its associative and initial is the identity element. The predicate must be a Boolean-valued function deciding whether a state is accepting or not.

This is an internal method encapsulating a common pattern for constraint representation.

source
TenSolver.dfa_to_mpoFunction
dfa_to_mpo([T], dfa, sites)

Build an exact diagonal projection MPO from a step-dependent DFA.

The MPO bond dimension is at most the number of states.

source
TenSolver.projection_mpoFunction
projection_mpo([T], constraint, sites; domain)

Build a projection MPO representing a constraint applicable to any MPS over sites. Constraint site numbers must use the same 1-based register indexing as sites.

Known constraints

  • SumConstraint uses a exact integer partial-sum automaton. For a constraint with rhs k, its maximum bond dimension is k+2.
  • SumModConstraint uses a modular partial-sum automaton. Its m residue states give it bond dimension m.
  • NotEqualsConstraint uses a MPO with bond dimension 2, independently of the rhs.
  • AssignmentConstraint uses a membership counting automaton. For rhs k, the maximum bond dimension is k+2.
  • RelationConstraint uses a MPO with bond dimension equal to the first variable's domain size.
source
TenSolver.projection_mposFunction
projection_mpos([T], constraints, sites; domain)

Build a list of projection MPOs representing constraints applicable to any MPS over sites.

This is a convenience wrapper around projection_mpo. T controls the numeric element type of the assembled MPO tensors.

source

Projected Hamiltonian Construction

TenSolver.project_hamiltonianFunction
project_hamiltonian(H, projections; formulation=:commuting, cutoff, kwargs...)

Project a Hamiltonian MPO with one or more projection MPOs.

If Q = P₁ * ⋯ * Pₙ is the combined projector, the effective Hamiltonian has the semantics Q' * H * Q.

With the default formulation=:commuting, H and all Pᵢ must be mutually commuting, while each Pᵢ must an orthogonal projection (Hermitian and idempotent). The construction then simplifies to H * Q, with bond dimension bounded by the product of H's links and each projection link. TenSolver's objective and constraint MPOs satisfy these assumptions because they are diagonal.

Use formulation=:sandwich for general, potentially noncommuting MPOs. It constructs Q' * H * Q directly, so each projection link contributes twice to the bond-dimension bound.

source
TenSolver.project_stateFunction
project_state(psi, projections; kwargs...)

Apply one or more diagonal projection MPOs to an MPS.

The result has zero amplitude on basis states rejected by any projection, while keeping the original unprimed site indices so it can be used as a DMRG input state.

source

Variable Domains

TenSolver.DomainsType
Domains{T}

Represent a finite domain for each variable in an optimization problem.

Indexing or Iterating over it yields the local variable domains, while "wholesale" operations such as in or rand treat it as a collection of possible values for a variable.

source
TenSolver.domain_residueFunction
domain_residue(p, domains)
domain_residue(Q, l, c, domains)

Simplify a (polynomial) function representation given that its variables are restricted to finite Domains.

The function name stems from calculating the residual of a polynomial modulo the ideal generated by the domain.

Theory

Each finite variable domain x_i in U_i = {u1, ..., ud} is equivalent to the root set of a single variable polynomial q_i(x) = (x_i - u1)...(x_i - ud).

By dividing p // q_i, we get

p(x) = m(x)q_i(x) + r(x).

Notice that for any a in U, q(a) = 0, and

p(a) = m(a)*0 + r(a) = r(a).

Thus, the transformation p -> r acts as degree reduction procedure. Performing it for all variables finds the residue.

source

Index