diff --git a/src/ITensorNetworksNext.jl b/src/ITensorNetworksNext.jl
index a8513d17..838d7339 100644
--- a/src/ITensorNetworksNext.jl
+++ b/src/ITensorNetworksNext.jl
@@ -14,8 +14,9 @@ include("AlgorithmsInterfaceExtensions/AlgorithmsInterfaceExtensions.jl")
 include("abstracttensornetwork.jl")
 include("tensornetwork.jl")
 include("itensornetworkoperator.jl")
-include("normnetwork.jl")
-include("normnetworkview.jl")
+include("bilinearforms/abstractbilinearformnetwork.jl")
+include("bilinearforms/normnetwork.jl")
+include("bilinearforms/quadraticformnetwork.jl")
 include("ITensorNetworkGenerators/ITensorNetworkGenerators.jl")
 include("contraction_tree.jl")
 include("prod_tensors.jl")
diff --git a/src/beliefpropagation/beliefpropagation.jl b/src/beliefpropagation/beliefpropagation.jl
index df39650a..fbb8f625 100644
--- a/src/beliefpropagation/beliefpropagation.jl
+++ b/src/beliefpropagation/beliefpropagation.jl
@@ -246,6 +246,12 @@ factor_tensors(factors, vertex) = [factors[vertex]]
 function factor_tensors(factors::NormNetwork, vertex)
     return [kettensor(factors, vertex), bratensor(factors, vertex)]
 end
+function factor_tensors(factors::QuadraticFormNetwork, vertex)
+    return [
+        kettensor(factors, vertex), operatortensor(factors, vertex),
+        bratensor(factors, vertex),
+    ]
+end
 
 # Contract the incoming messages into the source factor to form the (unnormalized) new message on
 # `edge`.
diff --git a/src/bilinearforms/abstractbilinearformnetwork.jl b/src/bilinearforms/abstractbilinearformnetwork.jl
new file mode 100644
index 00000000..0d4441b0
--- /dev/null
+++ b/src/bilinearforms/abstractbilinearformnetwork.jl
@@ -0,0 +1,189 @@
+using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned
+using Dictionaries: Dictionaries, isinsertable, issettable
+using Graphs: Graphs, edges, vertices
+using ITensorBase: ITensorBase, conj, inds, nametype, rename
+using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex
+
+"""
+    abstract type AbstractBilinearFormNetwork{T, V, I} <: AbstractITensorNetwork{T, V}
+
+Supertype of the lazy multi-layer networks built from a ket layer of type
+`ITensorNetwork{T, V, I}` and a ket→bra index name mapping.
+
+A subtype implements [`ketnetwork`](@ref) and [`branamemap`](@ref), from which the graph
+structure, [`braname`](@ref), [`kettensor`](@ref) and [`bratensor`](@ref) are defined, and,
+where it has an operator layer, [`operatornetwork`](@ref) and [`operatortensor`](@ref). The layers as whole networks are returned by [`ketnetwork`](@ref),
+[`branetwork`](@ref) and [`operatornetwork`](@ref).
+"""
+abstract type AbstractBilinearFormNetwork{T, V, I} <: AbstractITensorNetwork{T, V} end
+
+# The checks on a ket→bra name map shared by the bilinear-form constructors; `layers` are the
+# networks other than `ket` whose index names a bra name must not reuse.
+function check_branamemap(ket, braname, layers...)
+    if !issubset(linknames(ket), keys(braname))
+        error("every link name of the ket layer must have a bra name.")
+    end
+    if !allunique(braname)
+        error("the bra names of a bilinear-form network must be distinct.")
+    end
+    for name in braname
+        if has_dimname(ket, name) || any(layer -> has_dimname(layer, name), layers)
+            error("bra name $name is already an index name of another layer.")
+        end
+    end
+    return nothing
+end
+
+# ====================================== Graphs.jl ======================================= #
+
+Graphs.edges(bn::AbstractBilinearFormNetwork) = edges(ketnetwork(bn))
+Graphs.vertices(bn::AbstractBilinearFormNetwork) = vertices(ketnetwork(bn))
+
+# ==================================== NamedGraphs.jl ==================================== #
+
+function NamedGraphs.encoded_vertex(bn::AbstractBilinearFormNetwork, vertex)
+    return encoded_vertex(ketnetwork(bn), vertex)
+end
+function NamedGraphs.decoded_vertex(bn::AbstractBilinearFormNetwork, code::Integer)
+    return decoded_vertex(ketnetwork(bn), code)
+end
+NamedGraphs.encoded_graph(bn::AbstractBilinearFormNetwork) = encoded_graph(ketnetwork(bn))
+
+# ==================================== DataGraphs.jl ===================================== #
+
+function DataGraphs.is_vertex_assigned(bn::AbstractBilinearFormNetwork, vertex)
+    return isassigned(ketnetwork(bn), vertex)
+end
+
+# =================================== Dictionaries.jl ==================================== #
+
+Dictionaries.issettable(::AbstractBilinearFormNetwork) = false
+Dictionaries.isinsertable(::AbstractBilinearFormNetwork) = false
+
+# ====================================== interface ======================================= #
+
+# The name type of the indices, the same as the ket network's.
+function ITensorBase.nametype(
+        ::Type{<:AbstractBilinearFormNetwork{T, V, I}}
+    ) where {T, V, I}
+    return I
+end
+ITensorBase.nametype(bn::AbstractBilinearFormNetwork) = nametype(typeof(bn))
+
+"""
+    braname(bn::AbstractBilinearFormNetwork, name)
+
+The bra-layer index name corresponding to the ket-layer index name `name`.
+"""
+function braname(bn::AbstractBilinearFormNetwork, name)
+    if !has_dimname(ketnetwork(bn), name)
+        error("index name $name not found underlying tensor network.")
+    end
+    # A name absent from the map has no separate bra copy and maps to itself: a site index of a
+    # norm network, or a site index a quadratic form's operator does not act on.
+    return get(branamemap(bn), name, name)
+end
+
+"""
+    branamemap(bn::AbstractBilinearFormNetwork)
+
+The ket→bra name map, holding a bra name for each ket index name that has a separate bra copy.
+"""
+function branamemap end
+
+"""
+    kettensor(bn::AbstractBilinearFormNetwork, vertex)
+
+The ket-layer tensor at `vertex`.
+"""
+kettensor(bn::AbstractBilinearFormNetwork, vertex) = ketnetwork(bn)[vertex]
+
+"""
+    operatortensor(bn::AbstractBilinearFormNetwork, vertex)
+
+The operator-layer tensor at `vertex`, with its index names renamed so that its input legs
+meet the ket layer and its output legs meet the bra layer.
+"""
+function operatortensor end
+
+function conj_bratensor(bn::AbstractBilinearFormNetwork, vertex)
+    return rename(n -> braname(bn, n), kettensor(bn, vertex))
+end
+
+"""
+    bratensor(bn::AbstractBilinearFormNetwork, vertex)
+
+The bra-layer tensor at `vertex`.
+"""
+bratensor(bn::AbstractBilinearFormNetwork, vertex) = conj(conj_bratensor(bn, vertex))
+
+"""
+    ketnetwork(bn::AbstractBilinearFormNetwork)
+
+The ket-layer network of `bn`.
+"""
+function ketnetwork end
+
+"""
+    operatornetwork(bn::AbstractBilinearFormNetwork)
+
+The operator-layer network of `bn`, for a subtype that has an operator layer.
+"""
+function operatornetwork end
+
+"""
+    branetwork(bn::AbstractBilinearFormNetwork)
+
+The bra-layer network of `bn`. Unless a subtype stores its bra layer as a network, this is a
+`BraView`, whose tensors are built by [`bratensor`](@ref) when accessed.
+"""
+branetwork(bn::AbstractBilinearFormNetwork) = BraView(bn)
+
+"""
+    struct BraView{T, V, I, P <: AbstractBilinearFormNetwork{T, V, I}} <: AbstractITensorNetwork{T, V}
+
+The bra layer of the bilinear-form network `parent(view)`, with each vertex tensor built by
+[`bratensor`](@ref) when accessed. Its graph structure and mutability are those of the parent.
+"""
+struct BraView{T, V, I, P <: AbstractBilinearFormNetwork{T, V, I}} <:
+    AbstractITensorNetwork{T, V}
+    parent::P
+    function BraView(parent::AbstractBilinearFormNetwork{T, V, I}) where {T, V, I}
+        return new{T, V, I, typeof(parent)}(parent)
+    end
+end
+
+Base.parent(nnv::BraView) = nnv.parent
+
+# ==================================== DataGraphs.jl ===================================== #
+
+DataGraphs.get_vertex_data(nnv::BraView, vertex) = bratensor(parent(nnv), vertex)
+function DataGraphs.is_vertex_assigned(nnv::BraView, vertex)
+    return is_vertex_assigned(parent(nnv), vertex)
+end
+
+# ==================================== ITensorBase.jl ==================================== #
+
+function ITensorBase.inds(nnv::BraView, vertex)
+    return conj.(inds(conj_bratensor(parent(nnv), vertex)))
+end
+
+# ====================================== Graphs.jl ======================================= #
+
+Graphs.edges(nnv::BraView) = edges(parent(nnv))
+Graphs.vertices(nnv::BraView) = vertices(parent(nnv))
+
+# ==================================== NamedGraphs.jl ==================================== #
+
+function NamedGraphs.encoded_vertex(nnv::BraView, vertex)
+    return encoded_vertex(parent(nnv), vertex)
+end
+function NamedGraphs.decoded_vertex(nnv::BraView, code::Integer)
+    return decoded_vertex(parent(nnv), code)
+end
+NamedGraphs.encoded_graph(nnv::BraView) = encoded_graph(parent(nnv))
+
+# =================================== Dictionaries.jl ==================================== #
+
+Dictionaries.issettable(nnv::BraView) = issettable(parent(nnv))
+Dictionaries.isinsertable(nnv::BraView) = isinsertable(parent(nnv))
diff --git a/src/normnetwork.jl b/src/bilinearforms/normnetwork.jl
similarity index 51%
rename from src/normnetwork.jl
rename to src/bilinearforms/normnetwork.jl
index 6736adcf..8410ae49 100644
--- a/src/normnetwork.jl
+++ b/src/bilinearforms/normnetwork.jl
@@ -1,15 +1,15 @@
 using Dictionaries: Dictionary, Indices, dictionary
-using ITensorBase: ITensorBase, conj, nametype, rename, similar_operator
+using ITensorBase: similar_operator
 using ITensorNetworksNext
 
 """
-    struct NormNetwork{T, V, I} <: AbstractITensorNetwork{T, V}
+    struct NormNetwork{T, V, I} <: AbstractBilinearFormNetwork{T, V, I}
 
 Lazy wrapper representing the norm `⟨tn|tn⟩` of `tn::ITensorNetwork{T, V, I}`,
 together with a per-edge ket→bra name mapping that, for each index in the ket layer, defines
 the name of the corresponding index in the bra layer.
 """
-struct NormNetwork{T, V, I} <: AbstractITensorNetwork{T, V}
+struct NormNetwork{T, V, I} <: AbstractBilinearFormNetwork{T, V, I}
     ket::ITensorNetwork{T, V, I}
     braname::Dictionary{I, I}
     function NormNetwork(
@@ -31,18 +31,7 @@ function NormNetwork(tn::ITensorNetwork, linknames = ITensorNetworksNext.linknam
 end
 
 function check_input(::Type{<:NormNetwork}, ket, braname)
-    if !issubset(linknames(ket), keys(braname))
-        error("every link name of the ket layer must have a bra name.")
-    end
-    if !allunique(braname)
-        error("the bra names of a norm network must be distinct.")
-    end
-    for name in braname
-        if has_dimname(ket, name)
-            error("bra name $name is already an index name of the ket layer.")
-        end
-    end
-    return nothing
+    return check_branamemap(ket, braname)
 end
 
 function Base.eltype(::Type{<:NormNetwork})
@@ -51,21 +40,6 @@ function Base.eltype(::Type{<:NormNetwork})
     )
 end
 
-# ====================================== Graphs.jl ======================================= #
-
-Graphs.edges(nn::NormNetwork) = edges(ketnetwork(nn))
-Graphs.vertices(nn::NormNetwork) = vertices(ketnetwork(nn))
-
-# ==================================== NamedGraphs.jl ==================================== #
-
-function NamedGraphs.encoded_vertex(nn::NormNetwork, vertex)
-    return encoded_vertex(ketnetwork(nn), vertex)
-end
-function NamedGraphs.decoded_vertex(nn::NormNetwork, code::Integer)
-    return decoded_vertex(ketnetwork(nn), code)
-end
-NamedGraphs.encoded_graph(nn::NormNetwork) = encoded_graph(ketnetwork(nn))
-
 # ==================================== DataGraphs.jl ===================================== #
 
 function DataGraphs.get_vertex_data(nn::NormNetwork, vertex)
@@ -74,74 +48,11 @@ function DataGraphs.get_vertex_data(nn::NormNetwork, vertex)
     )
 end
 
-function DataGraphs.is_vertex_assigned(nn::NormNetwork, vertex)
-    return isassigned(ketnetwork(nn), vertex)
-end
-
-# =================================== Dictionaries.jl ==================================== #
-
-Dictionaries.issettable(::NormNetwork) = false
-Dictionaries.isinsertable(::NormNetwork) = false
-
 # ====================================== interface ======================================= #
 
-# The name type of the indices, the same as the ket network's.
-ITensorBase.nametype(::Type{<:NormNetwork{T, V, I}}) where {T, V, I} = I
-ITensorBase.nametype(nn::NormNetwork) = nametype(typeof(nn))
-
-"""
-    ketnetwork(nn::NormNetwork)
-
-The ket-layer network of `nn`.
-"""
 ketnetwork(nn::NormNetwork) = nn.ket
-
-"""
-    branetwork(nn::NormNetwork)
-
-The bra-layer network of `nn`, a `BraView` whose tensors are built by [`bratensor`](@ref) when
-accessed.
-"""
-branetwork(nn::NormNetwork) = BraView(nn)
-
-"""
-    branamemap(nn::NormNetwork)
-
-The ket→bra name map, holding a bra name for each ket index name that has a separate bra copy.
-"""
 branamemap(nn::NormNetwork) = nn.braname
 
-"""
-    braname(nn::NormNetwork, name)
-
-The bra-layer index name corresponding to the ket-layer index name `name`.
-"""
-function braname(nn::NormNetwork, name)
-    if !has_dimname(ketnetwork(nn), name)
-        error("index name $name not found underlying tensor network.")
-    end
-    # The names absent from the map are the site indices, which map to themselves.
-    return get(branamemap(nn), name, name)
-end
-
-"""
-    kettensor(nn::NormNetwork, vertex)
-
-The ket-layer tensor at `vertex`.
-"""
-kettensor(nn::NormNetwork, vertex) = ketnetwork(nn)[vertex]
-
-function conj_bratensor(nn::NormNetwork, vertex)
-    return rename(n -> braname(nn, n), kettensor(nn, vertex))
-end
-
-"""
-    bratensor(nn::NormNetwork, vertex)
-
-The bra-layer tensor at `vertex`.
-"""
-bratensor(nn::NormNetwork, vertex) = conj(conj_bratensor(nn, vertex))
-
 """
     flatten_network(nn::NormNetwork) -> ITensorNetwork
 
diff --git a/src/bilinearforms/quadraticformnetwork.jl b/src/bilinearforms/quadraticformnetwork.jl
new file mode 100644
index 00000000..169b0b18
--- /dev/null
+++ b/src/bilinearforms/quadraticformnetwork.jl
@@ -0,0 +1,124 @@
+using Dictionaries: Dictionary, Indices
+using ITensorBase: inputnames, outputnames, rename, state
+using ITensorNetworksNext
+
+"""
+    struct QuadraticFormNetwork{T, V, I, O} <: AbstractBilinearFormNetwork{T, V, I}
+
+Lazy wrapper representing the quadratic form `⟨tn|op|tn⟩` of `tn::ITensorNetwork{T, V, I}`
+sandwiched around the operator layer `op::ITensorNetworkOperator`, together with a per-index
+ket→bra name mapping that, for each index in the ket layer, defines the name of the
+corresponding index in the bra layer.
+
+The operator layer must be defined on every vertex of `tn`. Its input names are ket index
+names, and each paired output name is renamed to the bra layer, so the operator's remaining
+index names must be distinct from every index name of `tn`.
+"""
+struct QuadraticFormNetwork{T, V, I, O <: ITensorNetworkOperator} <:
+    AbstractBilinearFormNetwork{T, V, I}
+    ket::ITensorNetwork{T, V, I}
+    operator::O
+    braname::Dictionary{I, I}
+    function QuadraticFormNetwork(
+            ket::ITensorNetwork{T, V, I},
+            operator::ITensorNetworkOperator,
+            map::Dictionary{I, I}
+        ) where {T, V, I}
+        check_input(QuadraticFormNetwork, ket, operator, map)
+        return new{T, V, I, typeof(operator)}(ket, operator, map)
+    end
+end
+
+function QuadraticFormNetwork(
+        f::Function,
+        ket::ITensorNetwork,
+        operator::ITensorNetworkOperator,
+        linknames = ITensorNetworksNext.linknames(ket)
+    )
+    names = union(linknames, inputnames(operator))
+    return QuadraticFormNetwork(ket, operator, narrow_map(f, Indices(names)))
+end
+function QuadraticFormNetwork(
+        ket::ITensorNetwork,
+        operator::ITensorNetworkOperator,
+        linknames = ITensorNetworksNext.linknames(ket)
+    )
+    return QuadraticFormNetwork(name -> derivedname(name, "bra"), ket, operator, linknames)
+end
+
+function check_input(::Type{<:QuadraticFormNetwork}, ket, operator, braname)
+    if !issetequal(vertices(operator), vertices(ket))
+        error("the operator layer must be defined on every vertex of the ket layer.")
+    end
+    inputs = inputnames(operator)
+    for input in inputs
+        if !has_dimname(ket, input)
+            error("operator input name $input is not an index name of the ket layer.")
+        end
+        if dimnamevertices(ket, input) != dimnamevertices(state(operator), input)
+            error("operator input name $input must be a ket site index on the same vertex.")
+        end
+    end
+    for vertex in vertices(operator), name in names(state(operator), vertex)
+        if name ∉ inputs && has_dimname(ket, name)
+            error("operator index name $name is not an input but is a ket index name.")
+        end
+    end
+    if !issubset(inputs, keys(braname))
+        error("every operator input name must have a bra name.")
+    end
+    return check_branamemap(ket, braname, state(operator))
+end
+
+function Base.eltype(::Type{<:QuadraticFormNetwork})
+    return error(
+        "`eltype` of a `QuadraticFormNetwork` is not defined, since the triple-layer tensor \
+        at a vertex has no representation of its own. Use `kettensor`, `operatortensor` and \
+        `bratensor` to reach the individual layers."
+    )
+end
+
+# ==================================== DataGraphs.jl ===================================== #
+
+function DataGraphs.get_vertex_data(qf::QuadraticFormNetwork, vertex)
+    return error(
+        "Indexing a `QuadraticFormNetwork` is not defined, since the triple-layer tensor at \
+        a vertex has no representation of its own. Use `kettensor`, `operatortensor` and \
+        `bratensor` to reach the individual layers."
+    )
+end
+
+function DataGraphs.is_vertex_assigned(qf::QuadraticFormNetwork, vertex)
+    return isassigned(ketnetwork(qf), vertex) && isassigned(operatornetwork(qf), vertex)
+end
+
+# ====================================== interface ======================================= #
+
+ketnetwork(qf::QuadraticFormNetwork) = qf.ket
+branamemap(qf::QuadraticFormNetwork) = qf.braname
+operatornetwork(qf::QuadraticFormNetwork) = qf.operator
+
+# Each output name is renamed to the bra name of the input it is paired with, so the output legs
+# meet the bra layer and the input legs meet the ket layer.
+function operatortensor(qf::QuadraticFormNetwork, vertex)
+    op = qf.operator[vertex]
+    replacements = [
+        output => braname(qf, input) for
+            (output, input) in zip(outputnames(op), inputnames(op))
+    ]
+    return rename(state(op), replacements...)
+end
+
+"""
+    quadraticformnetwork([f::Function], tn::ITensorNetwork, op::ITensorNetworkOperator, [linknames])
+    quadraticformnetwork(tn::ITensorNetwork, op::ITensorNetworkOperator, braname::Dictionary)
+
+Build the triple-layer network `⟨tn|op|tn⟩`, represented lazily as a `QuadraticFormNetwork`
+object. Each ket index name in `linknames` (by default `linknames(tn)`) and each input name of
+`op` gets the bra name `f(name)`, where `f` defaults to `name -> derivedname(name, "bra")`.
+Alternatively pass the ket→bra name map `braname` directly. Either way the names mapped must
+include every link name of `tn` and every input name of `op` and may include others, and an
+error is thrown if one is missing, two bra names coincide or a bra name is already an index
+name of `tn` or `op`.
+"""
+quadraticformnetwork(args...) = QuadraticFormNetwork(args...)
diff --git a/src/normnetworkview.jl b/src/normnetworkview.jl
deleted file mode 100644
index 8b675089..00000000
--- a/src/normnetworkview.jl
+++ /dev/null
@@ -1,50 +0,0 @@
-using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned
-using Dictionaries: Dictionaries, isinsertable, issettable
-using Graphs: Graphs, edges, vertices
-using ITensorBase: ITensorBase, conj, inds
-using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex
-
-"""
-    struct BraView{T, V, I} <: AbstractITensorNetwork{T, V}
-
-The bra layer of the norm network `parent(view)`, with each vertex tensor built by
-[`bratensor`](@ref) when accessed. Its graph structure and mutability are those of the parent.
-"""
-struct BraView{T, V, I} <: AbstractITensorNetwork{T, V}
-    parent::NormNetwork{T, V, I}
-end
-
-Base.parent(nnv::BraView) = nnv.parent
-
-# ==================================== DataGraphs.jl ===================================== #
-
-DataGraphs.get_vertex_data(nnv::BraView, vertex) = bratensor(parent(nnv), vertex)
-function DataGraphs.is_vertex_assigned(nnv::BraView, vertex)
-    return is_vertex_assigned(parent(nnv), vertex)
-end
-
-# ==================================== ITensorBase.jl ==================================== #
-
-function ITensorBase.inds(nnv::BraView, vertex)
-    return conj.(inds(conj_bratensor(parent(nnv), vertex)))
-end
-
-# ====================================== Graphs.jl ======================================= #
-
-Graphs.edges(nnv::BraView) = edges(parent(nnv))
-Graphs.vertices(nnv::BraView) = vertices(parent(nnv))
-
-# ==================================== NamedGraphs.jl ==================================== #
-
-function NamedGraphs.encoded_vertex(nnv::BraView, vertex)
-    return encoded_vertex(parent(nnv), vertex)
-end
-function NamedGraphs.decoded_vertex(nnv::BraView, code::Integer)
-    return decoded_vertex(parent(nnv), code)
-end
-NamedGraphs.encoded_graph(nnv::BraView) = encoded_graph(parent(nnv))
-
-# =================================== Dictionaries.jl ==================================== #
-
-Dictionaries.issettable(nnv::BraView) = issettable(parent(nnv))
-Dictionaries.isinsertable(nnv::BraView) = isinsertable(parent(nnv))
diff --git a/test/test_quadraticformnetwork.jl b/test/test_quadraticformnetwork.jl
new file mode 100644
index 00000000..cfcbdc5d
--- /dev/null
+++ b/test/test_quadraticformnetwork.jl
@@ -0,0 +1,258 @@
+using DataGraphs: is_vertex_assigned
+using Dictionaries: isinsertable, issettable
+using Graphs: edges, vertices
+using ITensorBase: ITensor, Index, IndexName, conj, inds, inputnames, name, names, nametype,
+    operator, outputnames, rename, setname, uniquename
+using ITensorNetworksNext: BraView, Greedy, ITensorNetwork, NormNetwork,
+    QuadraticFormNetwork, braname, branetwork, bratensor, conj_bratensor, derivedname,
+    factor_tensors, ketnetwork, kettensor, operatornetwork, operatortensor, prod_tensors,
+    quadraticformnetwork, tensornetwork
+using LinearAlgebra: I, norm
+using NamedGraphs: NamedEdge, incident_edges, named_grid, named_path_graph
+using Test: @test, @test_throws, @testset
+
+# Build a random `ITensorNetwork` state on the graph `g` with site dimension `d` and
+# bond dimension `χ`.
+function random_state(::Type{T}, g; d = 2, χ = 2) where {T}
+    l = Dict(e => Index(χ) for e in edges(g))
+    l = merge(l, Dict(reverse(e) => l[e] for e in edges(g)))
+    s = Dict(v => Index(d) for v in vertices(g))
+    tn = tensornetwork(vertices(g)) do v
+        is = map(e -> l[e], incident_edges(g, v))
+        return randn(T, (s[v], is...))
+    end
+    return tn, l, s
+end
+
+# Build a bondless operator layer on the vertices of `g`, with `f(v)` supplying the matrix
+# acting on the site index `s[v]`.
+function product_operator(f, g, s; d = 2)
+    out = Dict(v => Index(d) for v in vertices(g))
+    tn = tensornetwork(vertices(g)) do v
+        return ITensor(f(v), (out[v], s[v]))
+    end
+    vs = collect(vertices(g))
+    return operator(tn, [name(out[v]) for v in vs], [name(s[v]) for v in vs])
+end
+
+identity_operator(g, s; d = 2) = product_operator(v -> Matrix(1.0I, d, d), g, s; d)
+
+function contract(bn)
+    return prod_tensors(reduce(vcat, factor_tensors(bn, v) for v in vertices(bn)), Greedy())
+end
+
+@testset "`QuadraticFormNetwork`" begin
+    @testset "Basics" begin
+        g = named_path_graph(3)
+        tn, l, s = random_state(Float64, g)
+        op = identity_operator(g, s)
+        qf = QuadraticFormNetwork(tn, op)
+
+        # `quadraticformnetwork` is the public constructor and agrees with the type.
+        @test quadraticformnetwork(tn, op) isa QuadraticFormNetwork
+        @test qf isa QuadraticFormNetwork
+        @test nametype(qf) === nametype(tn) === IndexName
+
+        # The quadratic form shares the graph structure of the ket layer.
+        @test issetequal(vertices(qf), vertices(tn))
+        @test issetequal(edges(qf), edges(tn))
+
+        # The triple-layer tensor at a vertex has no representation, so neither indexing nor
+        # `eltype` is defined.
+        @test_throws ErrorException qf[1]
+        @test_throws ErrorException eltype(qf)
+
+        # Vertex data is assigned wherever both layers are.
+        @test is_vertex_assigned(qf, 1)
+
+        # The quadratic form is neither settable nor insertable (it is a lazy view).
+        @test !issettable(qf)
+        @test !isinsertable(qf)
+
+        # The operator layer must cover every vertex of the ket layer.
+        g4 = named_path_graph(4)
+        _, _, s4 = random_state(Float64, g4)
+        @test_throws ErrorException QuadraticFormNetwork(tn, identity_operator(g4, s4))
+    end
+
+    @testset "operator input outside the ket network" begin
+        g = named_path_graph(2)
+        tn, l, s = random_state(Float64, g)
+        stray = Index(2)
+        out1, out2 = Index(2), Index(2)
+        optn = ITensorNetwork(Dict(1 => randn((out1, stray)), 2 => randn((out2, s[2]))))
+        op = operator(optn, [name(out1), name(out2)], [name(stray), name(s[2])])
+        @test_throws ErrorException QuadraticFormNetwork(tn, op)
+    end
+
+    @testset "operator input on a different vertex of the ket network" begin
+        g = named_path_graph(2)
+        tn, l, s = random_state(Float64, g)
+        out1, out2 = Index(2), Index(2)
+        optn = ITensorNetwork(Dict(1 => randn((out1, s[2])), 2 => randn((out2, s[1]))))
+        op = operator(optn, [name(out1), name(out2)], [name(s[2]), name(s[1])])
+        @test_throws ErrorException QuadraticFormNetwork(tn, op)
+    end
+
+    @testset "operator input equal to a ket link name" begin
+        g = named_path_graph(2)
+        tn, l, s = random_state(Float64, g)
+        e = only(edges(g))
+        out1, out2 = Index(2), Index(2)
+        optn = ITensorNetwork(Dict(1 => randn((out1, l[e])), 2 => randn((out2, s[2]))))
+        op = operator(optn, [name(out1), name(out2)], [name(l[e]), name(s[2])])
+        @test_throws ErrorException QuadraticFormNetwork(tn, op)
+    end
+
+    @testset "operator link name equal to a ket link name" begin
+        g = named_path_graph(2)
+        tn, l, s = random_state(Float64, g)
+        e = only(edges(g))
+        out1, out2 = Index(2), Index(2)
+        optn = ITensorNetwork(
+            Dict(1 => randn((out1, s[1], l[e])), 2 => randn((out2, s[2], l[e])))
+        )
+        op = operator(optn, [name(out1), name(out2)], [name(s[1]), name(s[2])])
+        @test_throws ErrorException QuadraticFormNetwork(tn, op)
+    end
+
+    @testset "a default bra name equal to an operator name" begin
+        g = named_path_graph(3)
+        tn, l, s = random_state(Float64, g)
+
+        # Each operator output takes the name the default map would derive for its input.
+        out =
+            Dict(v => setname(Index(2), derivedname(name(s[v]), :bra)) for v in vertices(g))
+        optn = tensornetwork(v -> ITensor(Matrix(1.0I, 2, 2), (out[v], s[v])), vertices(g))
+        vs = collect(vertices(g))
+        op = operator(optn, [name(out[v]) for v in vs], [name(s[v]) for v in vs])
+
+        @test_throws ErrorException QuadraticFormNetwork(tn, op)
+    end
+
+    @testset "layer tensors and the name map" begin
+        g = named_path_graph(3)
+        tn, l, s = random_state(Float64, g)
+        op = identity_operator(g, s)
+        qf = QuadraticFormNetwork(tn, op)
+
+        # `kettensor` returns the underlying tensor untouched.
+        @test kettensor(qf, 2) === tn[2]
+
+        # Unlike the norm network, the site indices *are* renamed in the bra layer: the
+        # operator sits between the two layers, so they no longer contract directly.
+        sname = name(s[2])
+        @test braname(qf, sname) != sname
+        @test sname in name.(inds(kettensor(qf, 2)))
+        @test !(sname in name.(inds(conj_bratensor(qf, 2))))
+        @test braname(qf, sname) in name.(inds(conj_bratensor(qf, 2)))
+
+        # Link indices are shared by two tensors, so they are renamed in the bra layer to
+        # keep the two layers' bonds distinct.
+        lname = name(l[NamedEdge(1 => 2)])
+        @test braname(qf, lname) != lname
+        @test lname in name.(inds(kettensor(qf, 2)))
+        @test !(lname in name.(inds(conj_bratensor(qf, 2))))
+        @test braname(qf, lname) in name.(inds(conj_bratensor(qf, 2)))
+
+        # The operator's input name meets the ket and its output name is renamed to meet
+        # the bra.
+        o = operatortensor(qf, 2)
+        @test sname in names(o)
+        @test braname(qf, sname) in names(o)
+        @test inputnames(op) == [name(s[v]) for v in vertices(g)]
+
+        # `bratensor` is the elementwise conjugate of `conj_bratensor` and carries the same
+        # indices.
+        @test inds(bratensor(qf, 2)) == inds(conj_bratensor(qf, 2))
+
+        # Querying the name map with an index name absent from the ket layer errors.
+        @test_throws ErrorException braname(qf, name(Index(2)))
+
+        # `factor_tensors` returns the three layers as separate operands.
+        @test factor_tensors(qf, 2) ==
+            [kettensor(qf, 2), operatortensor(qf, 2), bratensor(qf, 2)]
+    end
+
+    @testset "custom name map" begin
+        g = named_path_graph(3)
+        tn, l, s = random_state(Float64, g)
+        op = identity_operator(g, s)
+
+        # A user-supplied map dictates the bra-layer name for each renamed index.
+        custom = map(uniquename, keys(tn.dimname_vertices))
+        qf = quadraticformnetwork(tn, op, custom)
+
+        lname = name(l[NamedEdge(1 => 2)])
+        @test braname(qf, lname) == custom[lname]
+        @test braname(qf, lname) in name.(inds(conj_bratensor(qf, 2)))
+
+        sname = name(s[2])
+        @test braname(qf, sname) == custom[sname]
+    end
+
+    @testset "`ketnetwork` / `operatornetwork` / `branetwork`" begin
+        g = named_path_graph(3)
+        tn, l, s = random_state(Float64, g)
+        op = identity_operator(g, s)
+        qf = QuadraticFormNetwork(tn, op)
+
+        # The ket and operator layers are the networks the quadratic form was built from.
+        @test ketnetwork(qf) === tn
+        @test operatornetwork(qf) === op
+
+        # The bra layer is not stored, so it is a view sharing the ket layer's graph structure.
+        bv = branetwork(qf)
+        @test bv isa BraView
+        @test issetequal(vertices(bv), vertices(tn))
+        @test issetequal(edges(bv), edges(tn))
+        @test !issettable(bv)
+        @test !isinsertable(bv)
+        @test is_vertex_assigned(bv, 1)
+        for v in vertices(tn)
+            @test inds(bv[v]) == inds(bratensor(qf, v))
+        end
+    end
+
+    @testset "contraction / physics" begin
+        @testset "identity operator layer reproduces the norm network" begin
+            g = named_grid((2, 2))
+            tn, l, s = random_state(Float64, g)
+            op = identity_operator(g, s)
+
+            @test contract(QuadraticFormNetwork(tn, op))[] ≈
+                contract(NormNetwork(tn))[]
+        end
+
+        @testset "$T" for T in (Float64, ComplexF64)
+            g = named_grid((2, 2))
+            tn, l, s = random_state(T, g)
+
+            # A product of on-site matrices, contracted densely for the reference value.
+            mats = Dict(v => randn(T, 2, 2) for v in vertices(g))
+            op = product_operator(v -> mats[v], g, s)
+            qf = QuadraticFormNetwork(tn, op)
+
+            # ⟨ψ|O|ψ⟩ built by applying each on-site matrix to the dense ket and
+            # overlapping with the dense bra.
+            psi = prod(tn)
+            ket = psi
+            for v in vertices(g)
+                out = Index(2)
+                gate = ITensor(mats[v], (out, s[v]))
+                ket = rename(gate * ket, name(out) => name(s[v]))
+            end
+            @test contract(qf)[] ≈ (conj(psi) * ket)[]
+        end
+
+        @testset "scaled identity" begin
+            g = named_path_graph(3)
+            tn, l, s = random_state(Float64, g)
+            op = product_operator(v -> 2.0 * Matrix(1.0I, 2, 2), g, s)
+
+            # A factor of 2 at each of the three vertices scales the norm by 2³.
+            @test contract(QuadraticFormNetwork(tn, op))[] ≈
+                8 * contract(NormNetwork(tn))[]
+        end
+    end
+end