Squared Tensor Networks in Theory and Applications

• Squared tensor networks are doubled tensor-network constructions that fuse a tensor with its complex conjugate to compute quadratic observables like norms, overlaps, and density matrices.
• They underpin methods such as MPS, PEPS, and MERA, with the double-layer structure encoding correlations, entanglement, and phase information.
• These networks enable efficient computation in numerical linear algebra and probabilistic modeling by leveraging canonical forms and orthonormal constraints to reduce contraction complexity.

Squared tensor networks are doubled, or double-layer, tensor-network constructions obtained by taking a network for a state or operator and forming a new network from the tensor together with its complex conjugate, typically glued along physical legs. In this form, the objects computed are norms, overlaps, density matrices, transfer matrices, expectation values, and, in probabilistic settings, Born-style densities of the form $p(\mathbf{x}) = Z^{-1}|\psi(\mathbf{x})|^2$. The phrase “squared tensor network” is explicit in recent machine-learning work, whereas earlier tensor-network expositions systematically used the same structure without naming it: a single-layer network is replaced by a doubled network of local composites such as $A \otimes \bar A$, and the resulting contraction defines scalar quantities, reduced states, or effective channels whose spectra encode correlations, entanglement, and phase data (Bridgeman et al., 2016, Loconte et al., 2024).

1. Definition and diagrammatic structure

In tensor-network notation, a tensor is a node with one leg per index, and contraction is the joining of legs together with summation over the associated indices. The elementary operations are tensor product, trace or partial trace, contraction, and grouping or splitting indices so that a tensor can be treated as a matrix and subjected to linear-algebra operations such as SVD. A many-body pure state

$$|\psi\rangle = \sum_{j_1,\ldots,j_N} C_{j_1\ldots j_N}\,|j_1\ldots j_N\rangle$$

is represented by a rank-$N$ tensor $C$, and an MPS, PEPS, or MERA is precisely a factorization of that tensor into a network of smaller tensors with contracted virtual legs (Bridgeman et al., 2016).

A squared tensor network arises when a tensor network is paired with itself or with another network through an inner-product or outer-product construction. For a state $|\psi\rangle$, the norm

$\langle \psi|\psi\rangle$

is obtained by contracting the ket network with its bra, that is, its complex-conjugate network, along all physical indices. More generally, overlaps $\langle \phi|\psi\rangle$, density matrices $|\psi\rangle\langle\psi|$, reduced density matrices, Gram matrices, and quadratic objectives all have this doubled form. In the diagrammatic calculus, the same pattern yields Frobenius norms,

$$\|\mathcal{T}\|_F^2 = \langle \mathcal{T},\mathcal{T}\rangle,$$

inner products of tensors, and trace loops. The “Tensor Cookbook” makes these operations explicit as general tensor-network identities even though it does not name squared TNs as a separate class (Rakhshan et al., 15 May 2026).

This double-layer viewpoint is not merely notational. In a squared network, each local pair of tensors can be fused into a local double tensor, and the global contraction becomes the propagation of these composite objects through the network. This is the basis for transfer matrices in one dimension, effective environments in variational algorithms, and doubled operators in probability models. A common misconception is that squared TNs are a distinct architecture; in the sources considered here, they are better understood as a derived construction that recurs whenever one computes a quadratic quantity from an underlying tensor network (Bridgeman et al., 2016).

2. Manifestations in MPS, PEPS, MERA, TT, and HT

For a translationally invariant MPS on a periodic chain,

$A \otimes \bar A$0

squaring produces the transfer matrix

$A \otimes \bar A$1

and the norm becomes

$A \otimes \bar A$2

Expectation values are computed by replacing $A \otimes \bar A$3 at operator insertion sites with

$A \otimes \bar A$4

If the leading eigenvalue of $A \otimes \bar A$5 is normalized to $A \otimes \bar A$6, then connected two-point functions decay as $A \otimes \bar A$7, so exponential correlation decay is encoded directly in the spectrum of the squared transfer matrix. The density matrix $A \otimes \bar A$8 is likewise a double layer, and tracing out part of the system leaves an MPO whose virtual legs are doubled as ket-virtual $A \otimes \bar A$9 bra-virtual (Bridgeman et al., 2016).

For PEPS, the same construction gives a genuinely two-dimensional doubled network. A PEPS norm is obtained from local double tensors

$$|\psi\rangle = \sum_{j_1,\ldots,j_N} C_{j_1\ldots j_N}\,|j_1\ldots j_N\rangle$$0

one per lattice site, connected through virtual legs. Expectation values are computed by replacing bulk tensors with modified tensors $$|\psi\rangle = \sum_{j_1,\ldots,j_N} C_{j_1\ldots j_N}\,|j_1\ldots j_N\rangle$$1 where operators act. The contraction of this 2D double-layer network is generically $$|\psi\rangle = \sum_{j_1,\ldots,j_N} C_{j_1\ldots j_N}\,|j_1\ldots j_N\rangle$$2-hard in 2D, and the notes illustrate this through a mapping to counting graph colorings. Practical contractions therefore rely on approximations such as boundary MPS methods and corner transfer matrix methods, both of which act on the doubled PEPS rather than on the original single-layer state (Bridgeman et al., 2016).

For MERA, the squared construction is structurally different because unitary disentanglers $$|\psi\rangle = \sum_{j_1,\ldots,j_N} C_{j_1\ldots j_N}\,|j_1\ldots j_N\rangle$$3 and isometries $$|\psi\rangle = \sum_{j_1,\ldots,j_N} C_{j_1\ldots j_N}\,|j_1\ldots j_N\rangle$$4 satisfy

$$|\psi\rangle = \sum_{j_1,\ldots,j_N} C_{j_1\ldots j_N}\,|j_1\ldots j_N\rangle$$5

When expectation values are computed in the double-layer MERA, most of the network cancels outside the causal cone of the operator insertion. The surviving ladder of doubled tensors defines a channel

$$|\psi\rangle = \sum_{j_1,\ldots,j_N} C_{j_1\ldots j_N}\,|j_1\ldots j_N\rangle$$6

and its eigenvalues determine scaling dimensions through

$$|\psi\rangle = \sum_{j_1,\ldots,j_N} C_{j_1\ldots j_N}\,|j_1\ldots j_N\rangle$$7

Two-point functions then decay as a power law, with exponents determined by the spectral data of this squared local channel (Bridgeman et al., 2016).

The TT/HT literature presents the same pattern in numerical linear algebra. In TT format,

$$|\psi\rangle = \sum_{j_1,\ldots,j_N} C_{j_1\ldots j_N}\,|j_1\ldots j_N\rangle$$8

the squared norm $$|\psi\rangle = \sum_{j_1,\ldots,j_N} C_{j_1\ldots j_N}\,|j_1\ldots j_N\rangle$$9 is a double-layer TT contraction, while quadratic forms such as $N$0 become triple-layer networks with an MPO for $N$1 between two copies of the TT/MPS for $N$2. In this setting, squared TNs underlie least-squares solvers, eigenvalue computations, SVD, generalized eigenvalue problems, and trace-ratio objectives (Cichocki et al., 2017).

3. Probabilistic interpretation and squared circuits

Recent work recasts squared TNs as probabilistic models. Starting from a complex-valued tensor network or tensorized circuit $N$3 that represents an amplitude $N$4, one defines

$N$5

This is the Born-machine viewpoint: positivity is obtained by modulus-squaring a generally non-positive amplitude. The 2024 and 2025 papers generalize from TNs to tensorized circuits, a DAG formalism with input layers, Hadamard or Kronecker product layers, and sum layers, together with scope annotations that control tractable inference (Loconte et al., 2024, Loconte et al., 18 Dec 2025).

In the tensorized-circuit formalism, an input layer $N$6 squares to

$N$7

Hadamard layers square componentwise, Kronecker layers square up to a permutation, and sum layers square as

$N$8

If $N$9 is structured-decomposable, then $C$0 can itself be represented as a decomposable circuit. The crucial computational fact is that every layer’s size squares when passing from $C$1 to $C$2, so if $C$3 is the maximum layer size and $C$4 the number of layers, evaluation costs $C$5. In the scalar circuit analysis, the same blow-up appears as $C$6, with normalization or marginals costing $C$7 (Loconte et al., 2024, Loconte et al., 18 Dec 2025).

The computational bottleneck comes from the fact that the product of two arbitrary smooth and decomposable circuits is $C$8-hard to represent as another decomposable circuit unless the circuits are compatible, and compatibility with oneself is structured decomposability. This provides the circuit-theoretic analogue of the doubled-network contraction problem familiar from PEPS and other tensor-network models. The significance of the circuit framework is that it includes classical TN architectures such as MPS and tree tensor networks while also permitting more general DAG factorizations that do not correspond to a standard TN (Loconte et al., 2024, Loconte et al., 18 Dec 2025).

This probabilistic reinterpretation also clarifies the status of squared TNs in machine learning. They are expressive distribution estimators in high dimensions, but the act of squaring introduces exactly the normalization and marginalization overhead that canonical forms in physics were designed to avoid. Much of the recent theory therefore asks how to retain the probabilistic semantics of $C$9 without paying the full quadratic cost of explicitly constructing the square (Loconte et al., 2024).

4. Canonical forms, orthonormality, and unitarity

The principal remedy for the squaring overhead is to impose canonical or orthonormal structure on the underlying network. In MPS and TTN settings, canonical forms use isometric tensors so that contractions of $|\psi\rangle$0 collapse to products of identities. The circuit generalization makes this explicit. An orthonormal tensorized circuit requires, first, orthonormal input functions

$|\psi\rangle$1

and, second, semi-unitary sum weights

$|\psi\rangle$2

For structured-decomposable circuits, these conditions imply

$|\psi\rangle$3

so the squared model is normalized by construction (Loconte et al., 2024).

The 2025 treatment frames the same idea through unitary tensorized circuits. It imposes orthonormal input functions, basis decomposability, and semi-unitary sum weights, again with

$|\psi\rangle$4

Under these conditions, for every layer $|\psi\rangle$5 of output dimension $|\psi\rangle$6,

$|\psi\rangle$7

and for the scalar output layer this yields $|\psi\rangle$8. This is a direct abstraction of upper-canonical TTN and canonical MPS logic, but it applies to circuit factorizations that have no direct TN counterpart (Loconte et al., 18 Dec 2025).

A standard two-site example makes the mechanism concrete. For

$|\psi\rangle$9

naïve marginalization over $\langle \psi|\psi\rangle$0 requires $\langle \psi|\psi\rangle$1 overlap terms. If the left factors satisfy

$\langle \psi|\psi\rangle$2

then the marginal simplifies to

$\langle \psi|\psi\rangle$3

reducing the contraction from $\langle \psi|\psi\rangle$4 to $\langle \psi|\psi\rangle$5. This is the local form of what canonical gauges do globally in MPS and TTN models (Loconte et al., 18 Dec 2025).

An alternative route is determinism. For deterministic sum nodes with disjoint supports, squaring eliminates cross terms through

$\langle \psi|\psi\rangle$6

However, the 2025 paper emphasizes that this forces activations to become non-negative after squaring and thereby loses the benefits of complex or signed parameters. Orthogonality therefore generalizes determinism to the non-monotone setting: cross terms vanish in $\langle \psi|\psi\rangle$7 without requiring disjoint support in function space (Loconte et al., 18 Dec 2025).

5. Algorithms and computational roles

Squared TNs are the computational substrate of many standard algorithms. In DMRG, the variational energy

$\langle \psi|\psi\rangle$8

is expressed as a ratio of tensor networks whose denominator is a double-layer norm network and whose numerator is a bra-MPS–MPO–ket-MPS contraction. Fixing all tensors except one or two defines environment tensors $\langle \psi|\psi\rangle$9 and $\langle \phi|\psi\rangle$0, and after gauge fixing, the update reduces to an eigenvalue problem because the effective norm environment becomes the identity. TEBD is often implemented on the single-layer state, but norm computation, error analysis, and expectation values revert to the doubled form. MERA optimization likewise uses double-layer environments, simplified by causal-cone cancellations (Bridgeman et al., 2016).

In TT-based numerical linear algebra, squared TNs appear as effective local Hessians, Gram matrices, and Rayleigh quotients. Huge linear systems with cost

$\langle \phi|\psi\rangle$1

lead to quadratic forms involving $\langle \phi|\psi\rangle$2, and ALS/AMEn reduce the problem to local micro-problems

$\langle \phi|\psi\rangle$3

Eigenvalue, SVD, generalized eigenvalue, CCA, PLS, and regression objectives are all expressed as double-layer or triple-layer TT networks, often without explicitly materializing the squared object even though the underlying optimization is quadratic in the network variables (Cichocki et al., 2017).

For probabilistic squared circuits, the main algorithmic question is marginalization. In the orthonormal tensorized-circuit formalism, if $\langle \phi|\psi\rangle$4 are marginalized variables and $\langle \phi|\psi\rangle$5 are retained, then

$\langle \phi|\psi\rangle$6

can be computed in time

$\langle \phi|\psi\rangle$7

where $\langle \phi|\psi\rangle$8 are layers depending only on retained variables and $\langle \phi|\psi\rangle$9 are mixed-scope layers. Layers depending only on marginalized variables collapse to identity and never need to be explicitly squared (Loconte et al., 2024).

The scalar orthogonal-circuit analysis gives a complementary result: for a smooth, decomposable, orthogonal circuit, the partition function $|\psi\rangle\langle\psi|$0 can be computed in time $|\psi\rangle\langle\psi|$1. The corresponding marginalization algorithm propagates integrated squared values bottom-up and uses orthogonality to remove cross terms at sum nodes. In the tensorized unitary setting, mixed-scope sum layers are the only places where quadratic work remains, producing the complexity bound

$|\psi\rangle\langle\psi|$2

This expresses the same principle as MPS or TTN canonicalization: subcircuits wholly inside the marginalized region contract to identities, while only the interface between observed and marginalized variables needs a squared treatment (Loconte et al., 18 Dec 2025).

6. Expressiveness, applications, and open questions

A central question is whether orthonormal or unitary constraints reduce model expressiveness. The 2024 paper answers negatively for many circuit classes: if the input layers encode orthonormal functions, then there exists a polynomial-time orthonormalization algorithm returning an orthonormal circuit $|\psi\rangle\langle\psi|$3 and a matrix $|\psi\rangle\langle\psi|$4 such that

$|\psi\rangle\langle\psi|$5

Thus orthonormal parameterization is as expressive as arbitrary weights up to global normalization. The construction uses QR decomposition at sum layers and transforms Hadamard products into Kronecker structures when needed (Loconte et al., 2024).

The 2025 paper proves a related unitarization result: there is a polynomial-time algorithm that transforms a tensorized circuit into one with semi-unitary sum weights and output $|\psi\rangle\langle\psi|$6 for some scalar $|\psi\rangle\langle\psi|$7. It also shows that enforcing orthogonality on arbitrary decomposable circuits is $|\psi\rangle\langle\psi|$8-hard, which distinguishes “designing in” orthogonality from retrofitting it after the fact. This suggests that the tractable squared models are not obtained by a trivial reparameterization of all decomposable circuits, even though semi-unitarity alone does not reduce expressive efficiency (Loconte et al., 18 Dec 2025).

The main application domain in recent work is probabilistic modeling. Squared TNs and squared circuits are proposed as expressive distribution estimators supporting closed-form marginalization, with potential uses in lossless compression, sampling, symbolic reasoning, and tractable reasoning. For continuous variables, the papers list Fourier series, orthonormal polynomials, and Hermite functions as suitable orthonormal input bases. For non-tree or non-structured-decomposable architectures, unitary tensorized circuits extend canonical-form ideas beyond standard MPS and TTN settings (Loconte et al., 2024, Loconte et al., 18 Dec 2025).

Empirically, the 2025 study evaluates squared unitary probabilistic circuits on MNIST and FashionMNIST. It reports that squared unitary PCs avoid materializing the squared circuit and therefore scale more favorably in training time and GPU memory. One explicit example states that “a Kronecker-layer unitary model with 357M params uses 12 GiB and 0.29 ms/iter vs. 18 GiB and 0.52 ms for the baseline.” The same study reports that squared unitary PCs achieve essentially the same bits-per-dimension as unconstrained squared PCs and that non-structured-decomposable squared unitary PCs are competitive at large scale (Loconte et al., 18 Dec 2025).

Several limitations remain explicit. The complexity improvement

$|\psi\rangle\langle\psi|$9

is most beneficial when the mixed region is small; in densely entangled circuits, $$\|\mathcal{T}\|_F^2 = \langle \mathcal{T},\mathcal{T}\rangle,$$0 can approach the total number of layers. Orthogonalization can increase layer sizes, especially when Hadamard layers are replaced by Kronecker layers. Enforcing orthonormal constraints during training raises optimization questions on Stiefel manifolds. For continuous domains, the choice of orthonormal basis affects approximation quality and efficiency. More broadly, the present normalization and marginalization theory relies on structured decomposability or related orthogonality assumptions, and extension beyond these settings remains open (Loconte et al., 2024, Loconte et al., 18 Dec 2025).

Across condensed-matter tensor networks, tensorized numerical linear algebra, and probabilistic circuits, squared TNs therefore play a unifying role. They are the double-layer objects behind norms, expectation values, reduced states, transfer matrices, effective Hessians, and Born-style densities. Their computational difficulty is the difficulty of contraction after squaring; their tractability comes from canonical gauges, orthogonality, decomposability, and semi-unitary structure.