diff --git a/src/ITensorNetworksNext.jl b/src/ITensorNetworksNext.jl
index fe1b2540..ac62144b 100644
--- a/src/ITensorNetworksNext.jl
+++ b/src/ITensorNetworksNext.jl
@@ -13,8 +13,10 @@ include("select_algorithm.jl")
 include("AlgorithmsInterfaceExtensions/AlgorithmsInterfaceExtensions.jl")
 include("abstracttensornetwork.jl")
 include("tensornetwork.jl")
-include("normnetwork.jl")
-include("normnetworkview.jl")
+include("itensornetworkoperator.jl")
+include("bilinearforms/abstractbilinearformnetwork.jl")
+include("bilinearforms/normnetwork.jl")
+include("bilinearforms/quadraticformnetwork.jl")
 include("ITensorNetworkGenerators/ITensorNetworkGenerators.jl")
 include("contract_network.jl")
 
diff --git a/src/abstracttensornetwork.jl b/src/abstracttensornetwork.jl
index 341d7a6d..a8d959fa 100644
--- a/src/abstracttensornetwork.jl
+++ b/src/abstracttensornetwork.jl
@@ -133,9 +133,15 @@ function insertlink!(tn::AbstractGraph, e)
 end
 
 function operator_support(tn::AbstractGraph, op::ITensorOperator)
+    return operator_support_names(tn, inputnames(op))
+end
+
+# Shared with the `ITensorNetworkOperator` method, which is defined alongside that type
+# because it is not yet known at this point in the include order.
+function operator_support_names(tn::AbstractGraph, opnames)
     support = Indices{vertextype(tn)}()
 
-    for name in inputnames(op)
+    for name in opnames
         vertices = dimnamevertices(tn, name)
 
         if length(vertices) > 1
diff --git a/src/beliefpropagation/beliefpropagation.jl b/src/beliefpropagation/beliefpropagation.jl
index 3741fb3a..75f118f4 100644
--- a/src/beliefpropagation/beliefpropagation.jl
+++ b/src/beliefpropagation/beliefpropagation.jl
@@ -237,24 +237,13 @@ end
     contraction_alg::ContractionAlg = Exact()
 end
 
-# The tensors making up the factor at `vertex`, as separate operands for `contract_network`. A
-# `NormNetwork`'s factor is a lazy `ket * conj(bra)` product, and the contraction order sees each
-# operand as one node carrying only its outer axes — which hides the physical index the two layers
-# share, forcing the doubled vertex to be formed before any message is absorbed (χ^(2 * degree)
-# rather than the χ^(degree + 1) an interleaved order reaches).
-factor_tensors(factors, vertex) = [factors[vertex]]
-function factor_tensors(factors::NormNetwork, vertex)
-    return [kettensor(factors, vertex), bratensor(factors, vertex)]
-end
-
 # Contract the incoming messages into the source factor to form the (unnormalized) new message on
 # `edge`.
 function updated_message(algorithm::SimpleMessageUpdate, cache, factors, edge)
     messages = collect(incoming_messages(cache, edge))
-    # TODO: Remove `factor_tensors` once `contract_network` handles lazy tensors in
-    # contraction sequences properly.
     return contract_network(
-        [messages; factor_tensors(factors, src(edge))]; alg = algorithm.contraction_alg
+        [messages; [factors[src(edge)]]];
+        alg = algorithm.contraction_alg
     )
 end
 
@@ -275,9 +264,9 @@ end
 # the message is positive semidefinite and its trace is a positive normalization.
 function message_update!(algorithm::SimpleMessageUpdate, cache, factors::NormNetwork, edge)
     new_tensor = updated_message(algorithm, cache, factors, edge)
-    new_message = operator(
-        new_tensor, linknames(BraView(factors), edge), linknames(KetView(factors), edge)
-    )
+    branames = linknames(branetwork(factors), edge)
+    ketnames = linknames(ketnetwork(factors), edge)
+    new_message = operator(new_tensor, branames, ketnames)
     if algorithm.normalize
         message_norm = tr(new_message)
         iszero(message_norm) || (new_message /= message_norm)
diff --git a/src/beliefpropagation/messagecache.jl b/src/beliefpropagation/messagecache.jl
index 650a5858..be902628 100644
--- a/src/beliefpropagation/messagecache.jl
+++ b/src/beliefpropagation/messagecache.jl
@@ -132,10 +132,7 @@ end
 
 function vertex_scalar(factors, messages, vertex; kwargs...)
     in_messages = incoming_edge_data(messages, [vertex])
-    # TODO: Remove `factor_tensors` once `contract_network` handles lazy tensors in
-    # contraction sequences properly.
-    tensors = [factor_tensors(factors, vertex); collect(in_messages)]
-    return contract_network(tensors; kwargs...)[]
+    return contract_network([[factors[vertex]]; collect(in_messages)]; kwargs...)[]
 end
 
 vertex_scalars(factors, messages) = vertex_scalars(factors, messages, keys(factors))
@@ -193,16 +190,16 @@ bethe_free_energy(factors, messages) = -bethe_free_entropy(factors, messages)
 function similar_message_environment(nn::NormNetwork)
     messages = mapmany(vertices(nn)) do vertex
         return map(in_incident_edges(nn, vertex)) do edge
-            braview = BraView(nn)
-            ketview = KetView(nn)
+            bra = branetwork(nn)
+            ket = ketnetwork(nn)
 
-            ketnames = linknames(ketview, edge)
-            branames = linknames(braview, edge)
-            braaxis = unnamed.(linkaxes(braview, edge))
+            ketnames = linknames(ket, edge)
+            branames = linknames(bra, edge)
+            braaxis = unnamed.(linkaxes(bra, edge))
 
             # Bra leg = operator output, ket leg = input, the bipartition in which the message
             # is positive semidefinite.
-            message = similar_operator(ketview[vertex], braaxis, branames, ketnames)
+            message = similar_operator(ket[vertex], braaxis, branames, ketnames)
 
             return edge => message
         end
diff --git a/src/bilinearforms/abstractbilinearformnetwork.jl b/src/bilinearforms/abstractbilinearformnetwork.jl
new file mode 100644
index 00000000..f59823a9
--- /dev/null
+++ b/src/bilinearforms/abstractbilinearformnetwork.jl
@@ -0,0 +1,194 @@
+using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned
+using Dictionaries: Dictionaries, Dictionary, isinsertable, issettable
+using Graphs: Graphs, edges, vertices
+using ITensorBase: ITensorBase, conj, inds, name, 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 supplies its own graph structure, implements [`braname`](@ref), and returns an
+[`AbstractGramian`](@ref) from `getindex`; [`kettensor`](@ref), [`bratensor`](@ref) and,
+where the subtype has an operator layer, [`operatortensor`](@ref) read a vertex's layers from
+that Gramian. 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
+
+"""
+    abstract type AbstractGramian
+
+The layers of an `AbstractBilinearFormNetwork` at one vertex. A subtype implements
+[`kettensor`](@ref), [`braname`](@ref), `layertensors` and `layerinds`; the bra tensor is built
+from the ket tensor and the name map each time it is requested. Its `inds`, `names` and `axes`
+are the indices of its layers that no other layer shares, those the layer product leaves open.
+"""
+abstract type AbstractGramian 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 ======================================= #
+
+"""
+    braname(bn::AbstractBilinearFormNetwork, name)
+    braname(g::AbstractGramian, name)
+
+The bra-layer index name corresponding to the ket-layer index name `name`. The `AbstractGramian`
+form maps a name absent from its name map to itself, without checking it belongs to the network.
+"""
+function braname end
+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
+braname(g::AbstractGramian, name) = get(branamemap(g), name, name)
+
+"""
+    branamemap(bn::AbstractBilinearFormNetwork)
+    branamemap(g::AbstractGramian)
+
+The ket→bra name map, holding a bra name for each ket index name that has a separate bra copy.
+"""
+function branamemap end
+
+# A link name, or a name in `acted`, gets its bra name from `map`; every other name has none.
+function select_branames(ket::ITensorNetwork{T, V, I}, map, acted) where {T, V, I}
+    braname = Dictionary{I, I}()
+    for (name, vertices) in pairs(ket.dimname_vertices)
+        if length(vertices) == 2 || name in acted
+            insert!(braname, name, map[name])
+        end
+    end
+    return braname
+end
+
+"""
+    kettensor(g::AbstractGramian)
+
+The ket-layer tensor of the Gramian `g`.
+"""
+function kettensor end
+
+"""
+    operatortensor(g::AbstractGramian)
+
+The operator-layer tensor of the Gramian `g`, 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
+
+conj_bratensor(g::AbstractGramian) = rename(n -> braname(g, n), kettensor(g))
+
+"""
+    bratensor(g::AbstractGramian)
+
+The bra-layer tensor of the Gramian `g`.
+"""
+bratensor(g::AbstractGramian) = conj(conj_bratensor(g))
+
+# Read from `conj_bratensor`, which only renames, so the tensor data is not conjugated.
+brainds(g::AbstractGramian) = conj.(inds(conj_bratensor(g)))
+
+function ITensorBase.inds(g::AbstractGramian)
+    layer_inds = reduce(vcat, collect.(layerinds(g)))
+    layer_names = name.(layer_inds)
+    return [i for i in layer_inds if count(==(name(i)), layer_names) == 1]
+end
+ITensorBase.names(g::AbstractGramian) = name.(inds(g))
+Base.axes(g::AbstractGramian) = Tuple(inds(g))
+
+"""
+    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
+
+# ====================================== 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/bilinearforms/normnetwork.jl b/src/bilinearforms/normnetwork.jl
new file mode 100644
index 00000000..39da8153
--- /dev/null
+++ b/src/bilinearforms/normnetwork.jl
@@ -0,0 +1,65 @@
+using Dictionaries: Dictionary
+using ITensorBase: similar_operator, uniquename
+using ITensorNetworksNext
+
+"""
+    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} <: AbstractBilinearFormNetwork{T, V, I}
+    ket::ITensorNetwork{T, V, I}
+    braname::Dictionary{I, I}
+    function NormNetwork(
+            ket::ITensorNetwork{T, V, I},
+            map::Dictionary{I, I}
+        ) where {T, V, I}
+        return new{T, V, I}(ket, select_branames(ket, map, ()))
+    end
+end
+
+"""
+    struct NormGramian{T, I} <: AbstractGramian
+
+The layers of a `NormNetwork` at one vertex: the ket tensor and the network's ket→bra name map,
+from which the bra tensor is built when requested.
+"""
+struct NormGramian{T, I} <: AbstractGramian
+    ket::T
+    braname::Dictionary{I, I}
+end
+
+kettensor(g::NormGramian) = g.ket
+branamemap(g::NormGramian) = g.braname
+layertensors(g::NormGramian) = (; ket = kettensor(g), bra = bratensor(g))
+layerinds(g::NormGramian) = (inds(kettensor(g)), brainds(g))
+
+Base.eltype(::Type{<:NormNetwork{T, V, I}}) where {T, V, I} = NormGramian{T, I}
+
+function NormNetwork(tn::ITensorNetwork)
+    return NormNetwork(tn, map(uniquename, keys(tn.dimname_vertices)))
+end
+
+# ==================================== DataGraphs.jl ===================================== #
+
+function DataGraphs.get_vertex_data(nn::NormNetwork{T, V, I}, vertex) where {T, V, I}
+    return NormGramian{T, I}(nn.ket[vertex], nn.braname)
+end
+
+# ====================================== interface ======================================= #
+
+ketnetwork(nn::NormNetwork) = nn.ket
+branamemap(nn::NormNetwork) = nn.braname
+
+"""
+    normnetwork(tn::ITensorNetwork, [braname]) -> NormNetwork
+
+Build the double-layer norm network `⟨tn|tn⟩`, represented lazily as a `NomnNetwork` object.
+The optional second argument `braname` should implement `braname[ketdimname] = bradimname` for
+every link dimension name `ketdimname` in `tn`. If this is not specified, then a name is
+generated via the `ITensorBase.uniquename` function.
+"""
+normnetwork(tn::ITensorNetwork) = NormNetwork(tn)
+normnetwork(tn::ITensorNetwork, braname) = NormNetwork(tn, braname)
diff --git a/src/bilinearforms/quadraticformnetwork.jl b/src/bilinearforms/quadraticformnetwork.jl
new file mode 100644
index 00000000..ed50ae81
--- /dev/null
+++ b/src/bilinearforms/quadraticformnetwork.jl
@@ -0,0 +1,113 @@
+using Dictionaries: Dictionary
+using ITensorBase: inputnames, outputnames, rename, state, uniquename
+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}
+        if !issetequal(vertices(operator), vertices(ket))
+            error("the operator layer must be defined on every vertex of the ket layer.")
+        end
+        if !issubset(inputnames(operator), keys(ket.dimname_vertices))
+            error("every operator input name must be an index name of the ket layer.")
+        end
+        braname = select_branames(ket, map, Set{I}(inputnames(operator)))
+        return new{T, V, I, typeof(operator)}(ket, operator, braname)
+    end
+end
+
+"""
+    struct QuadraticFormGramian{T, O, I} <: AbstractGramian
+
+The layers of a `QuadraticFormNetwork` at one vertex: the ket tensor, the operator at that vertex
+and the ket→bra name map, from which the bra tensor and the renamed operator tensor are built
+when requested.
+"""
+struct QuadraticFormGramian{T, O, I} <: AbstractGramian
+    ket::T
+    operator::O
+    braname::Dictionary{I, I}
+end
+
+kettensor(g::QuadraticFormGramian) = g.ket
+branamemap(g::QuadraticFormGramian) = g.braname
+function layertensors(g::QuadraticFormGramian)
+    return (; ket = kettensor(g), operator = operatortensor(g), bra = bratensor(g))
+end
+function layerinds(g::QuadraticFormGramian)
+    return (inds(kettensor(g)), inds(operatortensor(g)), brainds(g))
+end
+
+function Base.eltype(::Type{<:QuadraticFormNetwork{T, V, I, O}}) where {T, V, I, O}
+    return QuadraticFormGramian{T, eltype(O), I}
+end
+
+function QuadraticFormNetwork(ket::ITensorNetwork, operator::ITensorNetworkOperator)
+    return QuadraticFormNetwork(ket, operator, map(uniquename, keys(ket.dimname_vertices)))
+end
+
+# ==================================== DataGraphs.jl ===================================== #
+
+function DataGraphs.get_vertex_data(
+        qf::QuadraticFormNetwork{T, V, I, O}, vertex
+    ) where {T, V, I, O}
+    return QuadraticFormGramian{T, eltype(O), I}(
+        qf.ket[vertex], qf.operator[vertex], qf.braname
+    )
+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(g::QuadraticFormGramian)
+    op = g.operator
+    replacements = [
+        output => braname(g, input) for
+            (output, input) in zip(outputnames(op), inputnames(op))
+    ]
+    return rename(state(op), replacements...)
+end
+
+"""
+    quadraticformnetwork(tn::ITensorNetwork, op::ITensorNetworkOperator, [braname])
+
+Build the triple-layer network `⟨tn|op|tn⟩`, represented lazily as a
+`QuadraticFormNetwork` object. The optional third argument `braname` should implement
+`braname[ketdimname] = bradimname` for every link dimension name `ketdimname` in `tn` and
+every dimension name of `tn` the operator acts on. If this is not specified, then a name is
+generated via the `ITensorBase.uniquename` function.
+"""
+function quadraticformnetwork(tn::ITensorNetwork, op::ITensorNetworkOperator)
+    return QuadraticFormNetwork(tn, op)
+end
+function quadraticformnetwork(tn::ITensorNetwork, op::ITensorNetworkOperator, braname)
+    return QuadraticFormNetwork(tn, op, braname)
+end
diff --git a/src/contract_network.jl b/src/contract_network.jl
index f13b36c8..c4897f3f 100644
--- a/src/contract_network.jl
+++ b/src/contract_network.jl
@@ -1,5 +1,6 @@
 using Base.Broadcast: materialize
 using Base: @kwdef
+using Dictionaries: Dictionary
 using ITensorBase: EvaluationOrderAlgorithm, Greedy, Mul, lazy, optimize_evaluation_order,
     substitute, symnamedtensor
 
@@ -30,12 +31,24 @@ function get_order(alg::Exact, tn)
         Dict(symnamedtensor(i) => symnamedtensor(i, Tuple(axes(t))) for (i, t) in pairs(tn))
     return substitute(order, subs)
 end
+# A Gramian enters as its separate layer tensors, so the order can place other operands between
+# layers; every operand gets a `(key, layer)` key so all keys share one concrete type.
+function split_gramians(tn)
+    any(t -> t isa AbstractGramian, tn) || return tn
+    pairs_split = [
+        (key, layer) => tensor for (key, t) in pairs(tn) for
+            (layer, tensor) in (t isa AbstractGramian ? pairs(layertensors(t)) : [:tensor => t])
+    ]
+    return Dictionary(first.(pairs_split), last.(pairs_split))
+end
+
 # Promote the operands to their common type before lowering to the lazy expression, so every lazy
 # operand shares one concrete type. Otherwise a network of mixed types (a plain tensor is a trivial
 # operator, so mixing operators and plain tensors is the common case) widens the symbolic `Mul`
 # container to a `UnionAll` it cannot construct. `promote_type`/`convert` keep an all-plain network
 # at the plain type (the promotion is a no-op), so its fast path is unchanged.
 function contract_network(alg::Exact, tn)
+    tn = split_gramians(tn)
     order = get_order(alg, tn)
     T = mapreduce(typeof, promote_type, tn)
     syms_to_ts = Dict(
@@ -48,7 +61,7 @@ end
 # `contraction_order`
 function contraction_order end
 function contraction_order(tn; alg = Greedy())
-    return contraction_order(alg, tn)
+    return contraction_order(alg, split_gramians(tn))
 end
 # Convert the tensor network to a flat symbolic multiplication expression.
 struct Flat end
diff --git a/src/itensornetworkoperator.jl b/src/itensornetworkoperator.jl
new file mode 100644
index 00000000..c8b68686
--- /dev/null
+++ b/src/itensornetworkoperator.jl
@@ -0,0 +1,117 @@
+using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned
+using Dictionaries: Dictionaries, Dictionary
+using Graphs: Graphs, AbstractGraph, edges, vertices
+using ITensorBase: ITensorBase, NamedTensorOperator, inputnames, names, nametype, operator,
+    outputnames, state
+using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex
+
+"""
+    struct ITensorNetworkOperator{T, V, I, P} <: AbstractITensorNetwork{T, V}
+
+The network equivalent of `ITensorBase.ITensorOperator`: a tensor network of type
+`P <: AbstractITensorNetwork{T, V}` together with a pairing of its dangling index names,
+stored as a map from each output name to its input name. Applying the operator contracts over
+the input names and leaves the output names.
+
+The output and input of each pair must sit on the same vertex. Indexing returns the vertex
+tensor wrapped as an `ITensorOperator` carrying the pairs at that vertex.
+"""
+struct ITensorNetworkOperator{T, V, I, P <: AbstractITensorNetwork{T, V}} <:
+    AbstractITensorNetwork{T, V}
+    parent::P
+    pairing::Dictionary{I, I}
+    function ITensorNetworkOperator(
+            parent::AbstractITensorNetwork{T, V}, outputnames, inputnames
+        ) where {T, V}
+        I = nametype(T)
+        outputnames = collect(I, outputnames)
+        inputnames = collect(I, inputnames)
+        if length(outputnames) != length(inputnames)
+            throw(
+                ArgumentError(
+                    "Operator `outputnames` and `inputnames` must have equal length " *
+                        "(positional pairing), got $(length(outputnames)) and " *
+                        "$(length(inputnames))."
+                )
+            )
+        end
+        if !allunique(outputnames)
+            throw(ArgumentError("each operator output name must appear only once."))
+        end
+        for opname in Iterators.flatten((outputnames, inputnames))
+            nvertices = length(dimnamevertices(parent, opname))
+            if nvertices != 1
+                throw(
+                    ArgumentError(
+                        "operator dim name $opname is associated with $nvertices vertices " *
+                            "in the tensor network; an operator leg must be a dangling index."
+                    )
+                )
+            end
+        end
+        for (output, input) in zip(outputnames, inputnames)
+            if only(dimnamevertices(parent, output)) != only(dimnamevertices(parent, input))
+                throw(
+                    ArgumentError(
+                        "operator output $output and its paired input $input must sit on " *
+                            "the same vertex."
+                    )
+                )
+            end
+        end
+        pairing = Dictionary(outputnames, inputnames)
+        return new{T, V, I, typeof(parent)}(parent, pairing)
+    end
+end
+
+function ITensorBase.operator(tn::AbstractITensorNetwork, outputnames, inputnames)
+    return ITensorNetworkOperator(tn, outputnames, inputnames)
+end
+
+function Base.eltype(::Type{<:ITensorNetworkOperator{T, V, I}}) where {T, V, I}
+    return NamedTensorOperator{I, T}
+end
+
+# ====================================== Graphs.jl ======================================= #
+
+Graphs.edges(op::ITensorNetworkOperator) = edges(state(op))
+Graphs.vertices(op::ITensorNetworkOperator) = vertices(state(op))
+
+# ==================================== NamedGraphs.jl ==================================== #
+
+function NamedGraphs.encoded_vertex(op::ITensorNetworkOperator, vertex)
+    return encoded_vertex(state(op), vertex)
+end
+function NamedGraphs.decoded_vertex(op::ITensorNetworkOperator, code::Integer)
+    return decoded_vertex(state(op), code)
+end
+NamedGraphs.encoded_graph(op::ITensorNetworkOperator) = encoded_graph(state(op))
+
+# ==================================== DataGraphs.jl ===================================== #
+
+# Both names of a pair sit on one vertex, so the pairs whose output is on `vertex` are its pairing.
+function DataGraphs.get_vertex_data(op::ITensorNetworkOperator, vertex)
+    tensor = state(op)[vertex]
+    outputs = filter(name -> haskey(op.pairing, name), names(tensor))
+    return operator(tensor, outputs, [op.pairing[output] for output in outputs])
+end
+
+function DataGraphs.is_vertex_assigned(op::ITensorNetworkOperator, vertex)
+    return is_vertex_assigned(state(op), vertex)
+end
+
+# =================================== Dictionaries.jl ==================================== #
+
+Dictionaries.issettable(::ITensorNetworkOperator) = false
+Dictionaries.isinsertable(::ITensorNetworkOperator) = false
+
+# ====================================== interface ======================================= #
+
+ITensorBase.state(op::ITensorNetworkOperator) = op.parent
+Base.parent(op::ITensorNetworkOperator) = state(op)
+ITensorBase.outputnames(op::ITensorNetworkOperator) = collect(keys(op.pairing))
+ITensorBase.inputnames(op::ITensorNetworkOperator) = collect(op.pairing)
+
+function operator_support(tn::AbstractGraph, op::ITensorNetworkOperator)
+    return operator_support_names(tn, inputnames(op))
+end
diff --git a/src/normnetwork.jl b/src/normnetwork.jl
deleted file mode 100644
index 46725507..00000000
--- a/src/normnetwork.jl
+++ /dev/null
@@ -1,92 +0,0 @@
-using Dictionaries: Dictionary
-using ITensorBase: LazyNamedTensor, lazy, rename, setname, similar_operator, uniquename
-using ITensorNetworksNext
-
-"""
-    struct NormNetwork{T, V, I} <: AbstractITensorNetwork{T, V}
-
-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}
-    ket::ITensorNetwork{T, V, I}
-    braname::Dictionary{I, I}
-    function NormNetwork(
-            ket::ITensorNetwork{T, V, I},
-            map::Dictionary{I, I}
-        ) where {T, V, I}
-        braname = Dictionary{I, I}()
-        for (name, vertices) in pairs(ket.dimname_vertices)
-            if length(vertices) == 2
-                insert!(braname, name, map[name])
-            end
-        end
-        return new{T, V, I}(ket, braname)
-    end
-end
-
-Base.eltype(::Type{<:NormNetwork{T, V, I}}) where {T, V, I} = LazyNamedTensor{I, T}
-
-function NormNetwork(tn::ITensorNetwork)
-    return NormNetwork(tn, map(uniquename, keys(tn.dimname_vertices)))
-end
-
-# ====================================== Graphs.jl ======================================= #
-
-Graphs.edges(nn::NormNetwork) = edges(nn.ket)
-Graphs.vertices(nn::NormNetwork) = vertices(nn.ket)
-
-# ==================================== NamedGraphs.jl ==================================== #
-
-NamedGraphs.encoded_vertex(nn::NormNetwork, vertex) = encoded_vertex(nn.ket, vertex)
-NamedGraphs.decoded_vertex(nn::NormNetwork, code::Integer) = decoded_vertex(nn.ket, code)
-NamedGraphs.encoded_graph(nn::NormNetwork) = encoded_graph(nn.ket)
-
-# ==================================== DataGraphs.jl ===================================== #
-
-function DataGraphs.get_vertex_data(nn::NormNetwork, vertex)
-    A = kettensor(nn, vertex)
-    B = conj_bratensor(nn, vertex)
-    # TODO: implement and use a lazy `conj` via `LazyNamedDimsArrays` here?
-    return lazy(A) * lazy(conj(B))
-end
-
-function DataGraphs.is_vertex_assigned(nn::NormNetwork, vertex)
-    return isassigned(nn.ket, vertex)
-end
-# =================================== Dictionaries.jl ==================================== #
-
-Dictionaries.issettable(::NormNetwork) = false
-Dictionaries.isinsertable(::NormNetwork) = false
-
-# ====================================== interface ======================================= #
-
-function braname(nn::NormNetwork, name)
-    if !has_dimname(nn.ket, name)
-        error("index name $name not found underlying tensor network.")
-    end
-    # The indices not stored in `nn.braname` are precisely the site indices, which
-    # get mapped to themselves.
-    return get(nn.braname, name, name)
-end
-
-indmap(nn::NormNetwork, ind) = setname(conj(ind), braname(nn, name(ind)))
-
-kettensor(nn::NormNetwork, vertex) = nn.ket[vertex]
-function conj_bratensor(nn::NormNetwork, vertex)
-    return rename(n -> braname(nn, n), kettensor(nn, vertex))
-end
-
-bratensor(nn::NormNetwork, vertex) = conj(conj_bratensor(nn, vertex))
-
-"""
-    normnetwork(tn::ITensorNetwork, [braname]) -> NormNetwork
-
-Build the double-layer norm network `⟨tn|tn⟩`, represented lazily as a `NomnNetwork` object.
-The optional second argument `braname` should implement `braname[ketdimname] = bradimname` for
-every link dimension name `ketdimname` in `tn`. If this is not specified, then a name is
-generated via the `ITensorBase.uniquename` function.
-"""
-normnetwork(tn::ITensorNetwork) = NormNetwork(tn)
-normnetwork(tn::ITensorNetwork, braname) = NormNetwork(tn, braname)
diff --git a/src/normnetworkview.jl b/src/normnetworkview.jl
deleted file mode 100644
index 75c57926..00000000
--- a/src/normnetworkview.jl
+++ /dev/null
@@ -1,58 +0,0 @@
-using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned
-using Dictionaries: Dictionaries, isinsertable, issettable
-using Graphs: Graphs, edges, vertices
-using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex
-
-struct KetView{T, V, I} <: AbstractITensorNetwork{T, V}
-    parent::NormNetwork{T, V, I}
-end
-
-struct BraView{T, V, I} <: AbstractITensorNetwork{T, V}
-    parent::NormNetwork{T, V, I}
-end
-
-# ====================================== Graphs.jl ======================================= #
-
-for View in (:KetView, :BraView)
-    @eval begin
-        Graphs.edges(nnv::$View) = edges(nnv.parent)
-        Graphs.vertices(nnv::$View) = vertices(nnv.parent)
-    end
-end
-
-# ==================================== NamedGraphs.jl ==================================== #
-
-for View in (:KetView, :BraView)
-    @eval begin
-        function NamedGraphs.encoded_vertex(nnv::$View, vertex)
-            return encoded_vertex(nnv.parent, vertex)
-        end
-        function NamedGraphs.decoded_vertex(nnv::$View, code::Integer)
-            return decoded_vertex(nnv.parent, code)
-        end
-
-        NamedGraphs.encoded_graph(nnv::$View) = encoded_graph(nnv.parent)
-    end
-end
-
-# ==================================== DataGraphs.jl ===================================== #
-
-DataGraphs.get_vertex_data(nn::KetView, vertex) = kettensor(nn.parent, vertex)
-DataGraphs.get_vertex_data(nn::BraView, vertex) = bratensor(nn.parent, vertex)
-
-for View in (:KetView, :BraView)
-    @eval begin
-        function DataGraphs.is_vertex_assigned(nnv::$View, vertex)
-            return isassigned(nnv.parent.ket, vertex)
-        end
-    end
-end
-
-# =================================== Dictionaries.jl ==================================== #
-
-for View in (:KetView, :BraView)
-    @eval begin
-        Dictionaries.issettable(nnv::$View) = issettable(nnv.parent)
-        Dictionaries.isinsertable(nnv::$View) = isinsertable(nnv.parent)
-    end
-end
diff --git a/test/test_beliefpropagation.jl b/test/test_beliefpropagation.jl
index f838f5d3..d86bb00e 100644
--- a/test/test_beliefpropagation.jl
+++ b/test/test_beliefpropagation.jl
@@ -1,5 +1,4 @@
 import AlgorithmsInterface as AI
-using Base.Broadcast: materialize
 using DataGraphs: DataGraphs, DataGraph, edge_data, edge_data_type
 using Dictionaries: Dictionary, dictionary, set!
 using GradedArrays: U1, gradedrange, isdual
@@ -10,9 +9,9 @@ using ITensorBase:
 using ITensorNetworksNext: ITensorNetworksNext, Exact, ITensorNetwork, MessageCache,
     NormNetwork, SimpleMessageUpdate, StopWhenConverged, beliefpropagation,
     bethe_free_energy, bethe_free_entropy, bratensor, contract_network, contraction_order,
-    edge_scalar, edge_scalars, factor_tensors, incoming_messages, insertlink!, kettensor,
-    linkaxes, linkinds, message_environment, messagecache, region_scalar, subgraph,
-    tensornetwork, updated_message, vertex_scalar, vertex_scalars
+    edge_scalar, edge_scalars, incoming_messages, insertlink!, kettensor, linkaxes,
+    linkinds, message_environment, messagecache, region_scalar, subgraph, tensornetwork,
+    updated_message, vertex_scalar, vertex_scalars
 using LinearAlgebra: LinearAlgebra, norm, tr
 using NamedGraphs: NamedEdge, all_edges, incident_edges, named_comb_tree, named_grid,
     named_path_graph, vertextype
@@ -344,7 +343,7 @@ end
             ones = message_environment(one, nn)
             for (edge, rest) in ((1 => 2, 2:4), (4 => 3, 1:3))
                 layers =
-                    [[kettensor(nn, v) for v in rest]; [bratensor(nn, v) for v in rest]]
+                    [[kettensor(nn[v]) for v in rest]; [bratensor(nn[v]) for v in rest]]
                 z_rest = contract_network([state(ones[edge]); layers])[]
                 @test z_rest ≈ norm(prod([network[v] for v in rest]))^2 rtol =
                     eps(real(T))^(1 / 3)
@@ -369,11 +368,10 @@ end
             nn = NormNetwork(network)
             v = (2, 2)
 
-            # A doubled vertex splits into its two layers, and the split is faithful.
-            @test length(factor_tensors(nn, v)) == 2
-            @test prod(factor_tensors(nn, v)) ≈ materialize(nn[v])
-            # A single-layer network's factor is a single operand.
-            @test factor_tensors(network, v) == [network[v]]
+            # A doubled vertex is a Gramian whose layers contract to the multiplied-out vertex.
+            gram = nn[v]
+            doubled = kettensor(gram) * bratensor(gram)
+            @test contract_network([gram]) ≈ doubled
 
             # The message update passes the layers to `contract_network` as separate operands, so
             # the contraction order can interleave the incoming messages between them, and the
@@ -388,7 +386,7 @@ end
             message = updated_message(algorithm, cache, nn, edge)
             # `v` has degree 4, so 3 incoming messages plus the ket and bra layers.
             @test only(counts) == 5
-            @test message ≈ contract_network([messages; [nn[v]]])
+            @test message ≈ contract_network([messages; [doubled]])
         end
     end
 end
diff --git a/test/test_contract_network.jl b/test/test_contract_network.jl
index 667b7d54..c5dd719f 100644
--- a/test/test_contract_network.jl
+++ b/test/test_contract_network.jl
@@ -1,8 +1,8 @@
 using Graphs: edges, vertices
 using ITensorBase:
     Greedy, Index, NamedTensorOperator, inputnames, operator, outputnames, state
-using ITensorNetworksNext: Exact, ITensorNetwork, LeftAssociative, contract_network,
-    linkinds, siteinds, tensornetwork
+using ITensorNetworksNext: ITensorNetworksNext, Exact, ITensorNetwork, LeftAssociative,
+    contract_network, linkinds, siteinds, tensornetwork
 using NamedGraphs: incident_edges, named_grid
 using OMEinsumContractionOrders: ExhaustiveSearch, GreedyMethod, TreeSA
 using Test: @test, @testset
@@ -17,6 +17,10 @@ using Test: @test, @testset
         C = [5.0, 1.0][j]
         D = [-2.0, 3.0, 4.0, 5.0, 1.0][k]
 
+        # A `Vector` of tensors holding no Gramian is returned unchanged.
+        ABCD = [A, B, C, D]
+        @test ITensorNetworksNext.split_gramians(ABCD) === ABCD
+
         ABCD_1 = contract_network([A, B, C, D]; alg = orderalg(LeftAssociative()))
         ABCD_2 = contract_network([A, B, C, D]; alg = orderalg(Greedy()))
         ABCD_3 = contract_network([A, B, C, D]; alg = orderalg(ExhaustiveSearch()))
@@ -37,6 +41,9 @@ using Test: @test, @testset
             return randn(Tuple(is))
         end
 
+        # A plain `ITensorNetwork` holds no Gramian, so `split_gramians` returns it unchanged.
+        @test ITensorNetworksNext.split_gramians(tn) === tn
+
         z1 = contract_network(tn; alg = orderalg(LeftAssociative()))[]
         z2 = contract_network(tn; alg = orderalg(Greedy()))[]
         z3 = contract_network(tn; alg = orderalg(ExhaustiveSearch()))[]
diff --git a/test/test_itensornetworkoperator.jl b/test/test_itensornetworkoperator.jl
new file mode 100644
index 00000000..ae6eb51c
--- /dev/null
+++ b/test/test_itensornetworkoperator.jl
@@ -0,0 +1,108 @@
+using DataGraphs: is_vertex_assigned
+using Dictionaries: isinsertable, issettable
+using Graphs: edges, vertices
+using ITensorBase: ITensor, Index, inputnames, name, operator, outputnames, state
+using ITensorNetworksNext:
+    ITensorNetwork, ITensorNetworkOperator, operator_support, tensornetwork
+using NamedGraphs: incident_edges, named_path_graph
+using Test: @test, @test_throws, @testset
+
+# A bondless operator network on `g`, with one input/output pair per vertex.
+function local_operator(g; d = 2)
+    inp = Dict(v => Index(d) for v in vertices(g))
+    out = Dict(v => Index(d) for v in vertices(g))
+    tn = tensornetwork(vertices(g)) do v
+        return randn((out[v], inp[v]))
+    end
+    vs = collect(vertices(g))
+    op = operator(tn, [name(out[v]) for v in vs], [name(inp[v]) for v in vs])
+    return op, inp, out
+end
+
+@testset "`ITensorNetworkOperator`" begin
+    @testset "Basics" begin
+        g = named_path_graph(3)
+        op, inp, out = local_operator(g)
+
+        @test op isa ITensorNetworkOperator
+
+        # The operator shares the graph structure of the network it wraps.
+        @test issetequal(vertices(op), vertices(state(op)))
+        @test issetequal(edges(op), edges(state(op)))
+        @test parent(op) === state(op)
+
+        # The pairing is positional: `outputnames[i]` goes with `inputnames[i]`.
+        vs = collect(vertices(g))
+        @test outputnames(op) == [name(out[v]) for v in vs]
+        @test inputnames(op) == [name(inp[v]) for v in vs]
+
+        @test is_vertex_assigned(op, first(vs))
+
+        # It is a lazy wrapper, so it is neither settable nor insertable.
+        @test !issettable(op)
+        @test !isinsertable(op)
+    end
+
+    @testset "`getindex` wraps the vertex tensor" begin
+        g = named_path_graph(3)
+        op, inp, out = local_operator(g)
+
+        # Both halves of each pair sit on the same vertex here, so the wrapper carries them.
+        for v in vertices(g)
+            @test outputnames(op[v]) == [name(out[v])]
+            @test inputnames(op[v]) == [name(inp[v])]
+            @test state(op[v]) === state(op)[v]
+        end
+
+        # `eltype` is the type of the wrapped vertex data.
+        @test eltype(op) === typeof(op[first(vertices(g))])
+    end
+
+    @testset "a pair must sit on one vertex" begin
+        # A swap pairs the output at vertex 1 with the input at vertex 2, and vice versa.
+        a, b = Index(2), Index(2)
+        a′, b′ = Index(2), Index(2)
+        tn = ITensorNetwork(Dict(1 => randn((a′, a)), 2 => randn((b′, b))))
+        @test_throws ArgumentError operator(tn, [name(a′), name(b′)], [name(b), name(a)])
+    end
+
+    @testset "constructor validation" begin
+        g = named_path_graph(2)
+        op, inp, out = local_operator(g)
+        tn = state(op)
+        vs = collect(vertices(g))
+
+        # The pairing is a bijection, so the two name lists must have equal length.
+        @test_throws ArgumentError operator(tn, [name(out[vs[1]])], [])
+        @test_throws ArgumentError operator(
+            tn, [name(out[v]) for v in vs], [name(inp[vs[1]])]
+        )
+        # An output name can be paired with only one input.
+        @test_throws ArgumentError operator(
+            tn, [name(out[vs[1]]), name(out[vs[1]])], [name(inp[vs[1]]), name(inp[vs[1]])]
+        )
+
+        # An operator leg must be dangling: a link name touches two vertices.
+        l = Index(2)
+        linked = tensornetwork(vertices(g)) do v
+            return randn((inp[v], l))
+        end
+        @test_throws ArgumentError operator(linked, [name(l)], [name(inp[vs[1]])])
+    end
+
+    @testset "`operator_support`" begin
+        g = named_path_graph(3)
+        s = Dict(v => Index(2) for v in vertices(g))
+        state_tn = tensornetwork(vertices(g)) do v
+            is = map(e -> Index(2), incident_edges(g, v))
+            return randn((s[v], is...))
+        end
+
+        # An operator acting on sites 1 and 3 is supported on exactly those vertices.
+        out1, out3 = Index(2), Index(2)
+        optn = ITensorNetwork(Dict(1 => randn((out1, s[1])), 3 => randn((out3, s[3]))))
+        op = operator(optn, [name(out1), name(out3)], [name(s[1]), name(s[3])])
+
+        @test issetequal(operator_support(state_tn, op), Set([1, 3]))
+    end
+end
diff --git a/test/test_normnetwork.jl b/test/test_normnetwork.jl
index 5007be69..ec1d64f8 100644
--- a/test/test_normnetwork.jl
+++ b/test/test_normnetwork.jl
@@ -1,11 +1,11 @@
 using DataGraphs: is_vertex_assigned
 using Dictionaries: isinsertable, issettable
 using Graphs: edges, vertices
-using ITensorBase:
-    ITensor, Index, IndexName, LazyITensor, conj, inds, name, setname, uniquename
-using ITensorNetworksNext: BraView, ITensorNetwork, KetView, NormNetwork, braname,
-    bratensor, conj_bratensor, contract_network, indmap, kettensor, normnetwork,
-    tensornetwork
+using ITensorBase: ITensor, Index, IndexName, LazyITensor, inds, name, names, uniquename
+using ITensorNetworksNext: ITensorNetworksNext, BraView, Exact, ITensorNetwork, NormGramian,
+    NormNetwork, braname, branetwork, bratensor, conj_bratensor, contract_network,
+    contraction_order, dimnamevertices, ketnetwork, kettensor, linkaxes, linkinds,
+    linknames, normnetwork, siteaxes, siteinds, sitenames, tensornetwork
 using LinearAlgebra: norm
 using NamedGraphs: NamedEdge, incident_edges, named_grid, named_path_graph
 using Test: @test, @test_throws, @testset
@@ -54,30 +54,25 @@ end
         nn = NormNetwork(tn)
 
         # `kettensor` returns the underlying tensor untouched.
-        @test kettensor(nn, 2) === tn[2]
+        @test kettensor(nn[2]) === tn[2]
 
         # Site indices appear in a single tensor, so they are *not* renamed: the ket and
         # bra layers share them (they get contracted, forming the physical overlap).
         sname = name(s[2])
         @test braname(nn, sname) == sname
-        @test sname in name.(inds(kettensor(nn, 2)))
-        @test sname in name.(inds(conj_bratensor(nn, 2)))
+        @test sname in name.(inds(kettensor(nn[2])))
+        @test sname in name.(inds(conj_bratensor(nn[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(nn, lname) != lname
-        @test lname in name.(inds(kettensor(nn, 2)))
-        @test !(lname in name.(inds(conj_bratensor(nn, 2))))
-        @test braname(nn, lname) in name.(inds(conj_bratensor(nn, 2)))
+        @test lname in name.(inds(kettensor(nn[2])))
+        @test !(lname in name.(inds(conj_bratensor(nn[2]))))
+        @test braname(nn, lname) in name.(inds(conj_bratensor(nn[2])))
 
         # `bra` is the elementwise conjugate of `conj_bratensor` and carries the same indices.
-        @test inds(bratensor(nn, 2)) == inds(conj_bratensor(nn, 2))
-
-        # `indmap` conjugates an index and renames it according to the name map.
-        ind = only(i for i in inds(kettensor(nn, 2)) if name(i) == lname)
-        @test name(indmap(nn, ind)) == braname(nn, name(ind))
-        @test indmap(nn, ind) == setname(conj(ind), braname(nn, name(ind)))
+        @test inds(bratensor(nn[2])) == inds(conj_bratensor(nn[2]))
 
         # Querying the name map with an index name absent from the network errors.
         @test_throws ErrorException braname(nn, name(Index(2)))
@@ -93,39 +88,82 @@ end
 
         lname = name(l[NamedEdge(1 => 2)])
         @test braname(nn, lname) == custom[lname]
-        @test braname(nn, lname) in name.(inds(conj_bratensor(nn, 2)))
+        @test braname(nn, lname) in name.(inds(conj_bratensor(nn[2])))
     end
 
-    @testset "`KetView` / `BraView`" begin
+    @testset "`ketnetwork` / `branetwork`" begin
         g = named_path_graph(3)
         tn, l, s = random_state(Float64, g)
         nn = NormNetwork(tn)
 
-        kv = KetView(nn)
-        bv = BraView(nn)
+        # The ket layer is the network the norm network was built from.
+        @test ketnetwork(nn) === tn
 
-        # Views share the graph structure of the underlying network.
-        @test issetequal(vertices(kv), vertices(tn))
+        # The bra layer is not stored, so it is a view sharing the ket layer's graph structure.
+        bv = branetwork(nn)
+        @test bv isa BraView
         @test issetequal(vertices(bv), vertices(tn))
-        @test issetequal(edges(kv), edges(tn))
         @test issetequal(edges(bv), edges(tn))
-
-        # The ket view exposes the bare ket tensors; the bra view exposes the bra tensors.
         for v in vertices(tn)
-            @test kv[v] === kettensor(nn, v)
-            @test inds(bv[v]) == inds(bratensor(nn, v))
+            @test inds(bv[v]) == inds(bratensor(nn[v]))
         end
-
-        @test is_vertex_assigned(kv, 1)
         @test is_vertex_assigned(bv, 1)
 
-        # Views inherit the (non-)mutability of their parent norm network.
-        @test !issettable(kv)
-        @test !isinsertable(kv)
+        # The view inherits the (non-)mutability of its parent norm network.
         @test !issettable(bv)
         @test !isinsertable(bv)
     end
 
+    @testset "`NormGramian`" begin
+        g = named_path_graph(3)
+        tn, l, s = random_state(Float64, g)
+        nn = NormNetwork(tn)
+        gram = nn[2]
+
+        @test gram isa NormGramian
+        @test eltype(nn) === typeof(gram)
+        # The Gramian holds the network's ket tensor and name map, not copies.
+        @test kettensor(gram) === tn[2]
+        @test gram.braname === nn.braname
+        @test inds(bratensor(gram)) == inds(conj_bratensor(gram))
+        @test keys(ITensorNetworksNext.layertensors(gram)) == (:ket, :bra)
+        # Contracting a Gramian contracts its layers.
+        @test contract_network([gram]) ≈ kettensor(gram) * bratensor(gram)
+        # A Gramian's indices are those its ket and bra layers do not share.
+        @test issetequal(inds(gram), inds(kettensor(gram) * bratensor(gram)))
+        @test names(gram) == name.(inds(gram))
+        @test axes(gram) == Tuple(inds(gram))
+    end
+
+    @testset "index queries on a `NormNetwork`" begin
+        g = named_path_graph(3)
+        tn, l, s = random_state(Float64, g)
+        nn = NormNetwork(tn)
+        e = NamedEdge(1 => 2)
+        lname = name(l[e])
+
+        # A link of the norm network is the ket link together with its bra-layer copy.
+        @test issetequal(linknames(nn, e), [lname, braname(nn, lname)])
+        @test issetequal(name.(linkinds(nn, e)), linknames(nn, e))
+        @test issetequal(name.(linkaxes(nn, e)), linknames(nn, e))
+        @test issetequal(dimnamevertices(nn, lname), [1, 2])
+        # The site index contracts between the layers, so no vertex has a site index.
+        @test isempty(siteinds(nn, 2))
+        @test isempty(sitenames(nn, 2))
+        @test isempty(siteaxes(nn, 2))
+    end
+
+    @testset "`contraction_order` on a `NormNetwork`" begin
+        g = named_path_graph(3)
+        tn, l, s = random_state(Float64, g)
+        nn = NormNetwork(tn)
+
+        # `contraction_order` splits the Gramians before computing an order, so it does not
+        # throw trying to call `size` on a `NormGramian`.
+        order = contraction_order(nn)
+        @test contract_network(nn; alg = Exact(; order))[] ≈ contract_network(nn)[]
+    end
+
     @testset "contraction / physics" begin
         @testset "single normalized tensor contracts to 1" begin
             s = Index(3)
diff --git a/test/test_quadraticformnetwork.jl b/test/test_quadraticformnetwork.jl
new file mode 100644
index 00000000..4dd3122e
--- /dev/null
+++ b/test/test_quadraticformnetwork.jl
@@ -0,0 +1,227 @@
+using DataGraphs: is_vertex_assigned
+using Dictionaries: isinsertable, issettable
+using Graphs: edges, vertices
+using ITensorBase: ITensor, Index, IndexName, conj, inds, inputnames, name, names, operator,
+    outputnames, rename, state, uniquename
+using ITensorNetworksNext: ITensorNetworksNext, BraView, ITensorNetwork, NormNetwork,
+    QuadraticFormGramian, QuadraticFormNetwork, braname, branetwork, bratensor,
+    conj_bratensor, contract_network, ketnetwork, kettensor, operatornetwork,
+    operatortensor, 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)
+
+@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
+
+        # The quadratic form shares the graph structure of the ket layer.
+        @test issetequal(vertices(qf), vertices(tn))
+        @test issetequal(edges(qf), edges(tn))
+
+        # `eltype` is the type of the (lazy triple-layer) vertex data.
+        @test eltype(qf) === typeof(qf[1])
+
+        # 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 "`QuadraticFormGramian`" begin
+        g = named_path_graph(3)
+        tn, l, s = random_state(Float64, g)
+        op = identity_operator(g, s)
+        qf = QuadraticFormNetwork(tn, op)
+        gram = qf[2]
+
+        @test gram isa QuadraticFormGramian
+        @test eltype(qf) === typeof(gram)
+        @test kettensor(gram) === tn[2]
+        # The Gramian holds the operator at its vertex, with that vertex's pairing.
+        @test state(gram.operator) === state(op)[2]
+        @test outputnames(gram.operator) == outputnames(op[2])
+        @test inputnames(gram.operator) == [name(s[2])]
+        @test keys(ITensorNetworksNext.layertensors(gram)) == (:ket, :operator, :bra)
+        @test contract_network([gram]) ≈
+            kettensor(gram) * operatortensor(gram) * bratensor(gram)
+        # A Gramian's indices are those no two of its layers share.
+        @test issetequal(
+            inds(gram), inds(kettensor(gram) * operatortensor(gram) * bratensor(gram))
+        )
+        @test names(gram) == name.(inds(gram))
+        @test axes(gram) == Tuple(inds(gram))
+    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 "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)))
+    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_network(QuadraticFormNetwork(tn, op))[] ≈
+                contract_network(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_network(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_network(QuadraticFormNetwork(tn, op))[] ≈
+                8 * contract_network(NormNetwork(tn))[]
+        end
+    end
+end