diff --git a/Project.toml b/Project.toml
index c82008c..2be8129 100644
--- a/Project.toml
+++ b/Project.toml
@@ -14,6 +14,14 @@ QUBODrivers = "a3f166f7-2cd3-47b6-9e1e-6fbfe0449eb0"
 QUBOTools = "60eb5b62-0a39-4ddc-84c5-97d2adff9319"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
+[weakdeps]
+SpinGlassEngine = "0563570f-ea1b-4080-8a64-041ac6565a4e"
+SpinGlassNetworks = "b7f6bd3e-55dc-4da6-96a9-ef9dbec6ac19"
+SpinGlassTensors = "7584fc6a-5a23-4eeb-8277-827aab0146ea"
+
+[extensions]
+TenSolverSpinGlassPEPSExt = ["SpinGlassEngine", "SpinGlassNetworks", "SpinGlassTensors"]
+
 [compat]
 Combinatorics = "1"
 ITensorMPS = "0.3"
@@ -22,6 +30,9 @@ LinearAlgebra = "1.10"
 MultivariatePolynomials = "0.5.9"
 Printf = "1"
 SparseArrays = "1.10"
+SpinGlassEngine = "1.6"
+SpinGlassNetworks = "1.4"
+SpinGlassTensors = "1.3"
 QUBODrivers = "0.6.1"
 QUBOTools = "0.13"
 julia = "1.10"
diff --git a/docs/src/api.md b/docs/src/api.md
index ae9daa8..c4053e1 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -7,6 +7,7 @@ This page documents the public API of TenSolver.jl.
 ```@docs
 TenSolver.minimize
 TenSolver.maximize
+TenSolver.solve_ising
 ```
 
 ## Solver Backends
diff --git a/docs/src/spinglasspeps_integration.md b/docs/src/spinglasspeps_integration.md
index 0a1f901..a8169d1 100644
--- a/docs/src/spinglasspeps_integration.md
+++ b/docs/src/spinglasspeps_integration.md
@@ -181,8 +181,31 @@ learn the SpinGlassPEPS API. The current behavior is:
 - `backend = :dmrg` selects the current path explicitly;
 - `backend = DMRGBackend()` selects the current path explicitly through the
   backend-object interface; and
-- `backend = :peps` errors clearly until a later optional bridge package or
-  extension provides the structured backend.
+- unavailable backend symbols error clearly without changing default DMRG
+  behavior.
+
+This stack step keeps the direct PEPS path as non-public scaffolding. The core
+package contains internal backend, topology, and result boundaries, while the
+exported `solve_ising` function is the public Ising boundary that optional
+structured backends may implement. `TenSolverSpinGlassPEPSExt` owns the
+SpinGlass component imports and calls. This keeps ordinary TenSolver installs
+on the existing dependency footprint and avoids documenting an activation path
+that cannot be tested from registered packages.
+
+The extension remains gated while the upstream dependency stack settles. In
+local checks against SpinGlassNetworks 1.4, SpinGlassEngine 1.6, and
+SpinGlassTensors 1.3, the current registered component compat bounds do not
+resolve with TenSolver's ITensors/QUBOTools environment. The source bridge and
+gated tests are kept in this stack step so the TenSolver boundary is concrete,
+but the PEPS backend types are not exported or listed in the public API until
+CI can exercise the SpinGlass component stack.
+
+The initial internal structured topology scaffolding covers one-spin-per-site
+and multi-spin-per-site square/king grids. QUBO inputs are converted through
+[`qubo_to_ising`](@ref) before the PEPS extension builds a SpinGlassNetworks
+Ising graph, clusters it with `super_square_lattice`, constructs the Potts
+Hamiltonian, runs `MpsContractor` plus `low_energy_spectrum`, and decodes
+retained states back to TenSolver Boolean vectors.
 
 Later PRs should add QUBODrivers/JuMP raw optimizer attributes for backend and
 PEPS parameters.
diff --git a/ext/TenSolverSpinGlassPEPSExt.jl b/ext/TenSolverSpinGlassPEPSExt.jl
new file mode 100644
index 0000000..053f900
--- /dev/null
+++ b/ext/TenSolverSpinGlassPEPSExt.jl
@@ -0,0 +1,263 @@
+module TenSolverSpinGlassPEPSExt
+
+import SparseArrays: findnz
+
+import TenSolver
+import TenSolver:
+  IsingModel,
+  KingGrid,
+  PEPSBackend,
+  PEPSSolution,
+  SquareGrid,
+  ising_energy,
+  spin_to_bool
+
+import SpinGlassEngine
+import SpinGlassNetworks
+import SpinGlassTensors
+
+function _peps_type(::IsingModel{T}) where {T}
+  return typeof(float(one(T)))
+end
+
+function _ising_instance(model::IsingModel{T}) where {T}
+  S = _peps_type(model)
+  instance = Dict{Tuple{Int, Int}, S}()
+
+  for i in eachindex(model.h)
+    instance[(i, i)] = S(model.h[i])
+  end
+
+  rows, cols, vals = findnz(model.J)
+  for k in eachindex(vals)
+    i = rows[k]
+    j = cols[k]
+    if i < j && !iszero(vals[k])
+      instance[(i, j)] = get(instance, (i, j), zero(S)) + S(vals[k])
+    end
+  end
+
+  return instance
+end
+
+function _check_topology_size(topology, model::IsingModel)
+  expected = TenSolver._topology_size(topology)
+  actual = length(model.h)
+  actual == expected ||
+    throw(DimensionMismatch("PEPS topology $(repr(topology)) expects $expected spins, but the Ising model has $actual spins."))
+end
+
+_edge_supported(::SquareGrid, a::Tuple, b::Tuple) = abs(a[1] - b[1]) + abs(a[2] - b[2]) == 1
+_edge_supported(::KingGrid, a::Tuple, b::Tuple) = maximum(abs.(a .- b)) == 1
+
+function _check_layout_edges(topology, model::IsingModel, lattice)
+  rows, cols, vals = findnz(model.J)
+  for k in eachindex(vals)
+    i = rows[k]
+    j = cols[k]
+    if i < j && !iszero(vals[k])
+      ci = lattice[i]
+      cj = lattice[j]
+      if ci != cj && !_edge_supported(topology, ci, cj)
+        throw(ArgumentError(
+          "Ising coupling ($i, $j) is not compatible with $(repr(topology)). " *
+          "Use a compatible structured topology or backend = :dmrg."
+        ))
+      end
+    end
+  end
+end
+
+function _potts_hamiltonian(backend::PEPSBackend, ig, lattice)
+  if isnothing(backend.local_dimension)
+    return SpinGlassNetworks.potts_hamiltonian(
+      ig;
+      spectrum = SpinGlassNetworks.full_spectrum,
+      cluster_assignment_rule = lattice,
+    )
+  end
+
+  return SpinGlassNetworks.potts_hamiltonian(
+    ig,
+    backend.local_dimension;
+    spectrum = SpinGlassNetworks.full_spectrum,
+    cluster_assignment_rule = lattice,
+  )
+end
+
+function _transformations(transformations)
+  transformations === :all && return SpinGlassEngine.all_lattice_transformations
+  transformations === :identity && return (SpinGlassEngine.rotation(0),)
+  if transformations isa Symbol
+    throw(ArgumentError("Unsupported PEPS transformations $(repr(transformations)). Use :all, :identity, a transformation, or a collection of transformations."))
+  end
+  transformations isa Tuple && return transformations
+  transformations isa AbstractVector && return Tuple(transformations)
+  return (transformations,)
+end
+
+function _strategy(backend::PEPSBackend)
+  backend.contraction in (:auto, :svd, :svd_truncate) && return SpinGlassEngine.SVDTruncate
+  backend.contraction === :zipper && return SpinGlassEngine.Zipper
+  throw(ArgumentError("Unsupported PEPS contraction $(repr(backend.contraction))."))
+end
+
+function _network(topology::SquareGrid, potts_h, transform, ::Type{T}) where {T}
+  return SpinGlassEngine.PEPSNetwork{
+    SpinGlassEngine.SquareSingleNode{SpinGlassEngine.GaugesEnergy},
+    SpinGlassEngine.Dense,
+    T,
+  }(topology.m, topology.n, potts_h, transform)
+end
+
+function _network(topology::KingGrid, potts_h, transform, ::Type{T}) where {T}
+  return SpinGlassEngine.PEPSNetwork{
+    SpinGlassEngine.KingSingleNode{SpinGlassEngine.GaugesEnergy},
+    SpinGlassEngine.Dense,
+    T,
+  }(topology.m, topology.n, potts_h, transform)
+end
+
+function _decoded_records(model::IsingModel, potts_h, sol, transform)
+  records = NamedTuple[]
+  for i in eachindex(sol.states)
+    decoded = SpinGlassNetworks.decode_potts_hamiltonian_state(potts_h, sol.states[i])
+    spins = [decoded[j] for j in eachindex(model.h)]
+    state = spin_to_bool(spins)
+    push!(records, (;
+      state,
+      spins,
+      energy = ising_energy(model, spins),
+      probability = sol.probabilities[i],
+      transformation = transform,
+      raw_energy = sol.energies[i],
+    ))
+  end
+  return records
+end
+
+function _deduplicated_records(records)
+  sort!(records; by = r -> (r.energy, -r.probability))
+
+  deduped = NamedTuple[]
+  positions = Dict{Any, Int}()
+  for record in records
+    key = Tuple(record.state)
+    index = get(positions, key, nothing)
+    if isnothing(index)
+      push!(deduped, record)
+      positions[key] = lastindex(deduped)
+    else
+      existing = deduped[index]
+      deduped[index] = (; existing..., probability = existing.probability + record.probability)
+    end
+  end
+
+  return deduped
+end
+
+function _metadata(backend::PEPSBackend, records, raw_results, failures)
+  best = first(records)
+  raw = raw_results[best.transformation]
+  return Dict{String, Any}(
+    "backend" => "SpinGlassPEPS",
+    "topology" => TenSolver._topology_name(backend.topology),
+    "topology_size" => TenSolver._topology_tuple(backend.topology),
+    "beta" => backend.beta,
+    "bond_dim" => backend.bond_dim,
+    "max_states" => backend.max_states,
+    "cutoff_prob" => backend.cutoff_prob,
+    "onGPU" => backend.onGPU,
+    "contraction" => String(backend.contraction),
+    "num_sweeps" => backend.num_sweeps,
+    "graduate_truncation" => backend.graduate_truncation,
+    "local_dimension" => backend.local_dimension,
+    "transformations_tried" => collect(string.(keys(raw_results))),
+    "transformations_failed" => [string(failure.transformation) for failure in failures],
+    "selected_transformation" => string(best.transformation),
+    "spin_glass_energies" => collect(raw.solution.energies),
+    "spin_glass_probabilities" => collect(raw.solution.probabilities),
+    "largest_discarded_probability" => raw.solution.largest_discarded_probability,
+  )
+end
+
+function TenSolver._solve_ising(backend::PEPSBackend, model::IsingModel; cutoff = nothing, verbosity = 1, kwargs...)
+  if !isempty(kwargs)
+    names = join(string.(keys(kwargs)), ", ")
+    throw(ArgumentError("Unsupported PEPS backend keyword(s): $names. Configure PEPSBackend instead."))
+  end
+
+  _check_topology_size(backend.topology, model)
+
+  S = _peps_type(model)
+  instance = _ising_instance(model)
+  ig = SpinGlassNetworks.ising_graph(S, instance)
+  lattice = SpinGlassNetworks.super_square_lattice(TenSolver._topology_tuple(backend.topology))
+  _check_layout_edges(backend.topology, model, lattice)
+  potts_h = _potts_hamiltonian(backend, ig, lattice)
+  params = SpinGlassEngine.MpsParameters{S}(;
+    bond_dim = backend.bond_dim,
+    num_sweeps = backend.num_sweeps,
+  )
+  search_params = SpinGlassEngine.SearchParameters(;
+    max_states = backend.max_states,
+    cutoff_prob = backend.cutoff_prob,
+  )
+  strategy = _strategy(backend)
+
+  records = NamedTuple[]
+  raw_results = Dict{Any, Any}()
+  failures = NamedTuple[]
+  for transform in _transformations(backend.transformations)
+    try
+      net = _network(backend.topology, potts_h, transform, S)
+      ctr = SpinGlassEngine.MpsContractor(
+        strategy,
+        net,
+        params;
+        onGPU = backend.onGPU,
+        beta = S(backend.beta),
+        graduate_truncation = backend.graduate_truncation,
+      )
+      merge_strategy = SpinGlassEngine.merge_branches(ctr; merge_prob = :none)
+      sol, info = SpinGlassEngine.low_energy_spectrum(
+        ctr,
+        search_params,
+        merge_strategy;
+        no_cache = backend.no_cache,
+      )
+
+      raw_results[transform] = (; solution = sol, info)
+      append!(records, _decoded_records(model, potts_h, sol, transform))
+    catch err
+      push!(failures, (;
+        transformation = transform,
+        error = sprint(showerror, err),
+      ))
+      verbosity > 0 && @warn "SpinGlassPEPS transformation failed" transformation = transform exception = (err, catch_backtrace())
+    finally
+      SpinGlassEngine.clear_memoize_cache()
+    end
+  end
+
+  if isempty(records)
+    if isempty(failures)
+      throw(ArgumentError("SpinGlassPEPS did not return any states."))
+    end
+
+    failure_summary = join(("$(failure.transformation): $(failure.error)" for failure in failures), "; ")
+    throw(ArgumentError("SpinGlassPEPS did not return any states. Failed transformations: $failure_summary"))
+  end
+  records = _deduplicated_records(records)
+  states = [record.state for record in records]
+  energies = S[record.energy for record in records]
+  probabilities = S[record.probability for record in records]
+  metadata = _metadata(backend, records, raw_results, failures)
+  raw = (; results = raw_results, failures)
+
+  verbosity > 0 && @info "SpinGlassPEPS backend finished" energy = first(energies) states = length(states)
+
+  return first(energies), PEPSSolution{S}(states, energies, probabilities, metadata, raw)
+end
+
+end # module
diff --git a/src/TenSolver.jl b/src/TenSolver.jl
index 6cf8339..28f6da8 100644
--- a/src/TenSolver.jl
+++ b/src/TenSolver.jl
@@ -14,7 +14,7 @@ include("ising.jl")
 export bool_to_spin, spin_to_bool, qubo_to_ising, ising_to_qubo
 
 include("solver.jl")
-export minimize, maximize
+export minimize, maximize, solve_ising
 export AbstractTenSolverBackend, DMRGBackend
 
 # Convergence logging
diff --git a/src/ising.jl b/src/ising.jl
index 7387769..c8e19e1 100644
--- a/src/ising.jl
+++ b/src/ising.jl
@@ -1,4 +1,52 @@
-import SparseArrays: dropzeros!, sparse
+import SparseArrays: SparseMatrixCSC, dropzeros!, findnz, sparse
+
+struct IsingModel{T <: Real}
+    J      :: SparseMatrixCSC{T, Int}
+    h      :: Vector{T}
+    offset :: T
+end
+
+function IsingModel(J::AbstractMatrix{<:Real}, h::AbstractVector{<:Real}, offset::Real=0)
+  issquare(J) || throw(DimensionMismatch("The Ising coupling matrix must be square. Encountered dimensions $(size(J))."))
+  size(J, 1) == length(h) || throw(DimensionMismatch("The Ising field vector length must match the coupling matrix size. Encountered dimensions $(size(J)) and length $(length(h))."))
+
+  T = promote_type(eltype(J), eltype(h), typeof(offset))
+  couplings, diagonal_offset = _canonical_ising_couplings(J, T)
+  return IsingModel{T}(couplings, collect(T, h), T(offset) + diagonal_offset)
+end
+
+function _canonical_ising_couplings(J::AbstractMatrix, ::Type{T}) where {T}
+  couplings = Dict{Tuple{Int, Int}, T}()
+  diagonal_offset = zero(T)
+  rows, cols, vals = findnz(sparse(T.(J)))
+
+  for k in eachindex(vals)
+    i = rows[k]
+    j = cols[k]
+    if i == j
+      diagonal_offset += vals[k]
+      continue
+    end
+
+    a, b = minmax(i, j)
+    key = (a, b)
+    couplings[key] = get(couplings, key, zero(T)) + vals[k]
+  end
+
+  out_rows = Int[]
+  out_cols = Int[]
+  out_vals = T[]
+  for (i, j) in sort!(collect(keys(couplings)))
+    coupling = couplings[(i, j)]
+    if !iszero(coupling)
+      push!(out_rows, i)
+      push!(out_cols, j)
+      push!(out_vals, coupling)
+    end
+  end
+
+  return sparse(out_rows, out_cols, out_vals, size(J, 1), size(J, 2)), diagonal_offset
+end
 
 function _conversion_type(Q::AbstractMatrix, l, c)
   T = promote_type(eltype(Q), isnothing(l) ? eltype(Q) : eltype(l), typeof(c))
@@ -53,6 +101,18 @@ function _drop_form_zeros!(form::QUBOTools.AbstractForm)
   return form
 end
 
+function _scaled_form_parts(form::QUBOTools.AbstractForm)
+  _, l, Q, scale, offset, _, _ = form
+  T = promote_type(eltype(Q), eltype(l), typeof(scale), typeof(offset))
+  return T(scale) .* sparse(T.(Q)), T(scale) .* collect(T, l), T(scale) * T(offset)
+end
+
+function IsingModel(form::QUBOTools.AbstractForm)
+  _check_form_domain(form, QUBOTools.SpinDomain, "Ising model")
+  J, h, offset = _scaled_form_parts(form)
+  return IsingModel(J, h, offset)
+end
+
 function _qubo_form(Q::AbstractMatrix, l::Union{Nothing, AbstractVector}, c::Real)
   _check_qubo_dimensions(Q, l)
 
@@ -138,6 +198,23 @@ function qubo_to_ising(form::QUBOTools.AbstractForm; convention::Symbol=:spin)
   return _drop_form_zeros!(QUBOTools.cast(QUBOTools.SpinDomain, form))
 end
 
+function ising_energy(model::IsingModel{T}, s::AbstractVector) where {T}
+  length(s) == length(model.h) || throw(DimensionMismatch("Spin vector length must match the Ising model size. Encountered length $(length(s)) and model size $(length(model.h))."))
+  spin_to_bool(s)
+
+  energy = model.offset + sum(model.h[i] * T(s[i]) for i in eachindex(model.h))
+  rows, cols, vals = findnz(model.J)
+  for k in eachindex(vals)
+    i = rows[k]
+    j = cols[k]
+    if i < j
+      energy += vals[k] * T(s[i]) * T(s[j])
+    end
+  end
+
+  return energy
+end
+
 """
     ising_to_qubo(form)
     ising_to_qubo(J, h[, offset])
@@ -163,6 +240,8 @@ function ising_to_qubo(form::QUBOTools.AbstractForm)
   return _drop_form_zeros!(QUBOTools.cast(QUBOTools.BoolDomain, form))
 end
 
+ising_to_qubo(model::IsingModel) = ising_to_qubo(_ising_form(model.J, model.h, model.offset))
+
 function ising_to_qubo(J::AbstractMatrix, h::AbstractVector, offset::Real=0)
   return ising_to_qubo(_ising_form(J, h, offset))
 end
diff --git a/src/solution.jl b/src/solution.jl
index 67bd63f..7af19ce 100644
--- a/src/solution.jl
+++ b/src/solution.jl
@@ -26,6 +26,21 @@ struct Solution{T <: Real}
     elapsed_times :: Vector{Float64}
 end
 
+# Internal result scaffold for the optional SpinGlassPEPS extension.
+#
+# This type is not exported while the optional SpinGlassPEPS component stack is
+# not covered by CI. The states are decoded to TenSolver Boolean vectors, and
+# `energies` are objective values in the original TenSolver convention.
+# Backend-specific diagnostics live in `metadata` and the raw extension result
+# is stored in `raw`.
+struct PEPSSolution{T <: Real}
+    states        :: Vector{Vector{Int}}
+    energies      :: Vector{T}
+    probabilities :: Vector{T}
+    metadata      :: Dict{String, Any}
+    raw           :: Any
+end
+
 # Sample from |ψ> in the {0, 1} world instead of 1-based Julia index world.
 """
     sample(psi)
@@ -36,13 +51,45 @@ sample(psi::Solution) = ITensorMPS.sample!(psi.tensor) .- 1
 
 sample(psi::Solution, n :: Integer) = [sample(psi) for _ in 1:n]
 
+function sample(psi::PEPSSolution)
+  isempty(psi.states) && throw(ArgumentError("Cannot sample an empty PEPS solution."))
+  if isempty(psi.probabilities)
+    return copy(first(psi.states))
+  end
+  length(psi.probabilities) == length(psi.states) ||
+    throw(ArgumentError("PEPS solution probabilities must match the number of retained states."))
+  any(probability -> probability < 0, psi.probabilities) &&
+    throw(ArgumentError("PEPS solution probabilities must be nonnegative."))
+
+  total = sum(psi.probabilities)
+  total > 0 || throw(ArgumentError("PEPS solution probabilities must have positive total weight."))
+
+  threshold = rand() * total
+  cumulative = zero(total)
+  for (state, probability) in zip(psi.states, psi.probabilities)
+    cumulative += probability
+    if threshold <= cumulative
+      return copy(state)
+    end
+  end
+
+  return copy(last(psi.states))
+end
+
+sample(psi::PEPSSolution, n :: Integer) = [sample(psi) for _ in 1:n]
+
 """
     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.
 """
-function Base.in(bs, psi::Solution; cutoff = 1e-8)
+function Base.in(bs::AbstractVector, psi::Solution; cutoff = 1e-8)
+  return prob(psi, bs) > cutoff
+end
+
+# Whether `xs` is one of the decoded Boolean states retained by the PEPS backend.
+function Base.in(bs::AbstractVector, psi::PEPSSolution; cutoff = 1e-8)
   return prob(psi, bs) > cutoff
 end
 
@@ -50,6 +97,17 @@ function prob(psi::Solution, bs)
   return abs2(coeff(psi, bs))
 end
 
+function prob(psi::PEPSSolution{T}, bs) where {T}
+  target = collect(Int, bs)
+  p = zero(T)
+  for (state, probability) in zip(psi.states, psi.probabilities)
+    if state == target
+      p += probability
+    end
+  end
+  return p
+end
+
 function coeff(psi::Solution, bs)
   tn    = psi.tensor
   sites = siteinds(tn)
diff --git a/src/solver.jl b/src/solver.jl
index 81d38ee..79de155 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -54,6 +54,110 @@ Select TenSolver's default ITensorMPS DMRG backend.
 """
 struct DMRGBackend <: AbstractTenSolverBackend end
 
+abstract type AbstractStructuredTopology end
+
+# Structured square-grid topology for optional PEPS solves.
+#
+# Variables are assumed to be ordered according to SpinGlassNetworks'
+# `super_square_lattice((m, n, spins_per_site))` convention.
+struct SquareGrid <: AbstractStructuredTopology
+  m :: Int
+  n :: Int
+  spins_per_site :: Int
+
+  function SquareGrid(m::Integer, n::Integer, spins_per_site::Integer=1)
+    m > 0 || throw(ArgumentError("SquareGrid requires m > 0. Got $m."))
+    n > 0 || throw(ArgumentError("SquareGrid requires n > 0. Got $n."))
+    spins_per_site > 0 || throw(ArgumentError("SquareGrid requires spins_per_site > 0. Got $spins_per_site."))
+    return new(Int(m), Int(n), Int(spins_per_site))
+  end
+end
+
+# Structured king-grid topology for optional PEPS solves. It uses the same
+# variable ordering as `SquareGrid`, but the PEPS compatibility graph also
+# allows diagonal interactions between neighboring grid cells.
+struct KingGrid <: AbstractStructuredTopology
+  m :: Int
+  n :: Int
+  spins_per_site :: Int
+
+  function KingGrid(m::Integer, n::Integer, spins_per_site::Integer=1)
+    m > 0 || throw(ArgumentError("KingGrid requires m > 0. Got $m."))
+    n > 0 || throw(ArgumentError("KingGrid requires n > 0. Got $n."))
+    spins_per_site > 0 || throw(ArgumentError("KingGrid requires spins_per_site > 0. Got $spins_per_site."))
+    return new(Int(m), Int(n), Int(spins_per_site))
+  end
+end
+
+_topology_size(topology::AbstractStructuredTopology) = topology.m * topology.n * topology.spins_per_site
+_topology_tuple(topology::AbstractStructuredTopology) = (topology.m, topology.n, topology.spins_per_site)
+_topology_name(::SquareGrid) = "square"
+_topology_name(::KingGrid) = "king"
+
+# Internal scaffold for the optional SpinGlassPEPS structured backend.
+#
+# The backend is implemented by the `TenSolverSpinGlassPEPSExt` package
+# extension, which loads only when `SpinGlassNetworks`, `SpinGlassEngine`, and
+# `SpinGlassTensors` are available. Without those packages this backend errors
+# clearly and the default DMRG backend remains unchanged.
+#
+# This constructor is intentionally not exported while the current registered
+# SpinGlass component dependency stack does not resolve in the same environment
+# as TenSolver's ITensors/QUBOTools stack and CI cannot exercise the extension.
+struct PEPSBackend{T <: AbstractStructuredTopology, S} <: AbstractTenSolverBackend
+  topology :: T
+  beta :: Float64
+  bond_dim :: Int
+  max_states :: Int
+  cutoff_prob :: Float64
+  onGPU :: Bool
+  contraction :: Symbol
+  num_sweeps :: Int
+  graduate_truncation :: Bool
+  transformations :: S
+  local_dimension :: Union{Nothing, Int}
+  no_cache :: Bool
+end
+
+function PEPSBackend(topology::AbstractStructuredTopology;
+                     beta::Real = 2.0,
+                     bond_dim::Integer = 16,
+                     max_states::Integer = 2^8,
+                     cutoff_prob::Real = 1e-4,
+                     onGPU::Bool = false,
+                     contraction::Symbol = :auto,
+                     num_sweeps::Integer = 1,
+                     graduate_truncation::Bool = true,
+                     transformations = :all,
+                     local_dimension::Union{Nothing, Integer} = nothing,
+                     no_cache::Bool = false)
+  beta > 0 && isfinite(beta) || throw(ArgumentError("PEPSBackend requires finite beta > 0. Got $beta."))
+  bond_dim >= 1 || throw(ArgumentError("PEPSBackend requires bond_dim >= 1. Got $bond_dim."))
+  max_states >= 1 || throw(ArgumentError("PEPSBackend requires max_states >= 1. Got $max_states."))
+  cutoff_prob >= 0 || throw(ArgumentError("PEPSBackend requires cutoff_prob >= 0. Got $cutoff_prob."))
+  num_sweeps >= 1 || throw(ArgumentError("PEPSBackend requires num_sweeps >= 1. Got $num_sweeps."))
+  contraction in (:auto, :svd, :svd_truncate, :zipper) ||
+    throw(ArgumentError("Unsupported PEPS contraction $(repr(contraction)). Use :auto, :svd, :svd_truncate, or :zipper."))
+  if !isnothing(local_dimension)
+    local_dimension >= 1 || throw(ArgumentError("PEPSBackend requires local_dimension >= 1 when provided. Got $local_dimension."))
+  end
+
+  return PEPSBackend{typeof(topology), typeof(transformations)}(
+    topology,
+    Float64(beta),
+    Int(bond_dim),
+    Int(max_states),
+    Float64(cutoff_prob),
+    onGPU,
+    contraction,
+    Int(num_sweeps),
+    graduate_truncation,
+    transformations,
+    isnothing(local_dimension) ? nothing : Int(local_dimension),
+    no_cache,
+  )
+end
+
 const default_backend = DMRGBackend()
 
 function _backend_error(backend)
@@ -64,7 +168,10 @@ function _backend_error(backend)
   return ArgumentError("No `_minimize` method is available for backend $(repr(backend)). Use backend = :dmrg or provide a backend-specific `_minimize` method.")
 end
 
+_backend_error(::PEPSBackend) = ArgumentError("PEPSBackend is not available. Install/load SpinGlassNetworks, SpinGlassEngine, and SpinGlassTensors to activate the PEPS extension, or use backend = :dmrg.")
+
 _normalize_backend(backend::DMRGBackend) = backend
+_normalize_backend(backend::PEPSBackend) = backend
 _normalize_backend(backend::AbstractTenSolverBackend) = backend
 _normalize_backend(backend::Symbol) = _normalize_backend(Val(backend))
 _normalize_backend(::Val{:dmrg}) = default_backend
@@ -227,6 +334,44 @@ function _minimize(backend::AbstractTenSolverBackend, args...; kwargs...)
   throw(_backend_error(backend))
 end
 
+function _minimize(backend::PEPSBackend, Q::AbstractMatrix{T}, l::Union{AbstractVector{T}, Nothing}=nothing, c::T=zero(T); kwargs...) where T
+  return _solve_ising(backend, IsingModel(qubo_to_ising(Q, l, c)); kwargs...)
+end
+
+function _minimize(backend::PEPSBackend, p::AbstractPolynomial; kwargs...)
+  throw(ArgumentError("PEPSBackend does not support polynomial inputs directly. Convert to a structured QUBO or call solve_ising with a supported topology."))
+end
+
+"""
+    solve_ising(model; backend = DMRGBackend(), kwargs...)
+    solve_ising(J, h[, offset]; backend = DMRGBackend(), kwargs...)
+
+Solve an Ising model with spins `s_i in {-1, +1}`.
+
+The returned solution still samples TenSolver Boolean vectors using
+`x_i = (s_i + 1) / 2`. The default DMRG path converts the Ising model back to a
+QUBO and calls [`minimize`](@ref). Optional structured backends can implement
+this boundary directly.
+"""
+function solve_ising end
+
+function solve_ising(model::IsingModel; backend=default_backend, kwargs...)
+  return _solve_ising(_normalize_backend(backend), model; kwargs...)
+end
+
+function solve_ising(J::AbstractMatrix, h::AbstractVector, offset::Real=0; backend=default_backend, kwargs...)
+  return solve_ising(IsingModel(J, h, offset); backend, kwargs...)
+end
+
+function _solve_ising(backend::AbstractTenSolverBackend, model::IsingModel; kwargs...)
+  throw(_backend_error(backend))
+end
+
+function _solve_ising(::DMRGBackend, model::IsingModel; kwargs...)
+  Q, l, c = _scaled_form_parts(ising_to_qubo(model))
+  return _minimize(default_backend, Q, l, c; kwargs...)
+end
+
 function _minimize(::DMRGBackend, Q::AbstractMatrix{T}, l::Union{AbstractVector{T}, Nothing}=nothing, c::T=zero(T); cutoff=1e-8, kwargs...) where T
   H      = tensorize(Q, isnothing(l) ? diag(Q) : diag(Q) + l; cutoff)
   obj(x) = dot(x, Q, x) + c + maybe(l -> dot(l,x), l; default=zero(T))
diff --git a/test/backend.jl b/test/backend.jl
index d20b532..be777bf 100644
--- a/test/backend.jl
+++ b/test/backend.jl
@@ -30,6 +30,30 @@ end
     @test TenSolver.sample(psi) == [1]
   end
 
+  @testset "solve_ising preserves pair couplings and offsets" begin
+    J = [0.0 0.5; 1.0 0.0]
+    h = [-0.25, 0.75]
+    offset = 2.0
+    model = TenSolver.IsingModel(J, h, offset)
+    spin_states = [[s1, s2] for s1 in (-1, 1) for s2 in (-1, 1)]
+    energies = [TenSolver.ising_energy(model, spin) for spin in spin_states]
+    expected_energy, expected_index = findmin(energies)
+    expected_spin = spin_states[expected_index]
+
+    for solve_call in (
+      () -> solve_ising(model; backend=:dmrg, verbosity=0),
+      () -> solve_ising(J, h, offset; backend=:dmrg, verbosity=0),
+    )
+      energy, solution = solve_call()
+      sample_bits = TenSolver.sample(solution)
+      sample_spin = TenSolver.bool_to_spin(sample_bits)
+
+      @test energy ≈ expected_energy
+      @test sample_spin == expected_spin
+      @test TenSolver.ising_energy(model, sample_spin) ≈ expected_energy
+    end
+  end
+
   @testset "Symbol backends can be provided by extensions" begin
     Q = reshape([1.0], 1, 1)
     E, payload = minimize(Q; backend=:test_symbol_backend, verbosity=0)
diff --git a/test/peps_backend.jl b/test/peps_backend.jl
new file mode 100644
index 0000000..4fe9930
--- /dev/null
+++ b/test/peps_backend.jl
@@ -0,0 +1,145 @@
+@testset "PEPS backend core" begin
+  @test !(:Solution in names(TenSolver))
+  @test !(:PEPSBackend in names(TenSolver))
+  @test !(:SquareGrid in names(TenSolver))
+  @test !(:KingGrid in names(TenSolver))
+  @test !(:PEPSSolution in names(TenSolver))
+
+  @test TenSolver.SquareGrid(2, 3).m == 2
+  @test TenSolver.SquareGrid(2, 3).n == 3
+  @test TenSolver.SquareGrid(2, 3).spins_per_site == 1
+  @test TenSolver.KingGrid(2, 3, 2).spins_per_site == 2
+  @test_throws ArgumentError TenSolver.SquareGrid(0, 1)
+  @test_throws ArgumentError TenSolver.KingGrid(1, 0)
+
+  backend = TenSolver.PEPSBackend(
+    TenSolver.SquareGrid(1, 1);
+    beta = 1.5,
+    bond_dim = 4,
+    max_states = 2,
+    cutoff_prob = 0.0,
+    onGPU = false,
+    contraction = :svd,
+    transformations = :identity,
+  )
+
+  @test backend.topology == TenSolver.SquareGrid(1, 1)
+  @test backend.beta == 1.5
+  @test backend.bond_dim == 4
+  @test backend.max_states == 2
+  @test backend.cutoff_prob == 0.0
+  @test backend.contraction == :svd
+  @test_throws ArgumentError TenSolver.PEPSBackend(TenSolver.SquareGrid(1, 1); beta = 0)
+  @test_throws ArgumentError TenSolver.PEPSBackend(TenSolver.SquareGrid(1, 1); bond_dim = 0)
+  @test_throws ArgumentError TenSolver.PEPSBackend(TenSolver.SquareGrid(1, 1); max_states = 0)
+  @test_throws ArgumentError TenSolver.PEPSBackend(TenSolver.SquareGrid(1, 1); contraction = :unknown)
+
+  Q = reshape([-1.0], 1, 1)
+  peps_error = try
+    minimize(Q; backend, verbosity = 0)
+  catch err
+    err
+  end
+  @test peps_error isa ArgumentError
+  @test occursin("PEPSBackend is not available", sprint(showerror, peps_error))
+  @test occursin("SpinGlassNetworks", sprint(showerror, peps_error))
+
+  model = TenSolver.IsingModel(zeros(1, 1), [-1.0])
+  energy, solution = solve_ising(model; backend = :dmrg, verbosity = 0)
+  @test energy ≈ -1.0
+  @test sample(solution) == [1]
+
+  peps_solution = TenSolver.PEPSSolution{Float64}(
+    [[1, 0], [0, 1]],
+    [-2.0, -1.0],
+    [0.0, 1.0],
+    Dict{String, Any}("backend" => "SpinGlassPEPS"),
+    nothing,
+  )
+  @test sample(peps_solution) == [0, 1]
+  @test sample(peps_solution, 2) == [[0, 1], [0, 1]]
+  @test [0, 1] in peps_solution
+  @test !([1, 0] in peps_solution)
+  @test !([0, 0] in peps_solution)
+  @test TenSolver.prob(peps_solution, [0, 1]) ≈ 1.0
+
+  duplicate_state_solution = TenSolver.PEPSSolution{Float64}(
+    [[1, 0], [0, 1], [1, 0]],
+    [-2.0, -1.0, -2.0],
+    [0.2, 0.3, 0.4],
+    Dict{String, Any}("backend" => "SpinGlassPEPS"),
+    nothing,
+  )
+  @test TenSolver.prob(duplicate_state_solution, [1, 0]) ≈ 0.6
+  @test TenSolver.prob(duplicate_state_solution, [0, 1]) ≈ 0.3
+
+  @test_throws ArgumentError sample(TenSolver.PEPSSolution{Float64}(
+    [[1, 0], [0, 1]],
+    [-2.0, -1.0],
+    [1.0],
+    Dict{String, Any}(),
+    nothing,
+  ))
+  @test_throws ArgumentError sample(TenSolver.PEPSSolution{Float64}(
+    [[1, 0], [0, 1]],
+    [-2.0, -1.0],
+    [1.0, -0.5],
+    Dict{String, Any}(),
+    nothing,
+  ))
+  @test_throws ArgumentError sample(TenSolver.PEPSSolution{Float64}(
+    [[1, 0], [0, 1]],
+    [-2.0, -1.0],
+    [0.0, 0.0],
+    Dict{String, Any}(),
+    nothing,
+  ))
+end
+
+@testset "Optional SpinGlassPEPS extension" begin
+  has_spinglasspeps_components = all(pkg -> !isnothing(Base.find_package(pkg)), (
+    "SpinGlassNetworks",
+    "SpinGlassEngine",
+    "SpinGlassTensors",
+  ))
+
+  if !has_spinglasspeps_components
+    @test_skip "SpinGlassPEPS component packages are not available in this environment."
+  else
+    import SpinGlassNetworks
+    import SpinGlassEngine
+    import SpinGlassTensors
+
+    backend = TenSolver.PEPSBackend(
+      TenSolver.SquareGrid(2, 2);
+      beta = 2.0,
+      bond_dim = 4,
+      max_states = 4,
+      cutoff_prob = 0.0,
+      onGPU = false,
+      contraction = :svd,
+      transformations = :identity,
+    )
+
+    Q = [
+      -1.0 0.5 0.0 0.0
+       0.0 -0.5 0.0 0.0
+       0.0 0.0 -0.25 0.25
+       0.0 0.0 0.0 -0.75
+    ]
+    l = [0.0, 0.25, -0.25, 0.0]
+    c = 0.125
+    objective(x) = dot(x, Q, x) + dot(l, x) + c
+    exact_energy, _ = brute_force(objective, Int, 4)
+
+    energy, solution = minimize(Q, l, c; backend, verbosity = 0)
+    state = sample(solution)
+
+    @test energy ≈ exact_energy atol = 1e-6
+    @test objective(state) ≈ energy atol = 1e-6
+    @test solution.metadata["backend"] == "SpinGlassPEPS"
+    @test solution.metadata["topology"] == "square"
+    @test solution.metadata["selected_transformation"] == string(SpinGlassEngine.rotation(0))
+    @test first(solution.energies) ≈ energy atol = 1e-6
+  end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 6a6b7df..4069e6d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -13,6 +13,7 @@ include(filepath("utils.jl"))
 
 include(filepath("ising_conversion.jl"))
 include(filepath("backend.jl"))
+include(filepath("peps_backend.jl"))
 
 @testset "Ill-formed input" begin
   dim = 4