1 Introduction↩︎

Since its derivation in 1926 by Erwin Schrödinger, the electronic Schrödinger equation has been studied extensively, mainly in order to predict molecular bindings between atoms or the structure of larger molecules. In numerical mathematics and theoretical physics/chemistry, the main task is to compute the ground state energy of the equation, that is, the smallest eigenvalue of the Hamiltonian operator. While this might seem straightforward, the problem size scales exponentially with the number of electrons involved, and therefore, standard eigensolvers fail even for small systems of \(10\)-\(20\) electrons. Before the advent of crypto-currencies and AI, computations aiming at solving the Schrödinger problem commonly used about a third of the computation power of some large supercomputers [1].

Many approximation methods for energy minimization of the electronic Schrödinger equation have been proposed by physicists and mathematicians, among them the Hartree-Fock method, which treats the involved electrons as uncorrelated, and density functional theory, where the two-particle term of the Hamiltonian is replaced by a one-particle potential. See [2] for a mathematical introduction into both. For correlated many-body systems, widely used methods are the coupled cluster method [3]–[5], where a surrogate problem is solved on a lower-dimensional algebraic variety, and the density matrix renormalization group (DMRG) algorithm, which approximates the eigenstates by a low-rank tensor [6]. The latter was proposed by theoretical physicists along with matrix product states to approximate solutions of the electronic Schrödinger equation in second quantization [7]. The realization that this approach is actually a low-rank tensor method emerged slowly and separately for physicists [8] and mathematicians [9]. In mathematics, the matrix product states have been dubbed tensor trains (TT) corresponding to the tensor train rank (TT rank) of a tensor, and the DMRG algorithm has been recognized as a two-site version of the alternating linear scheme (ALS) [10]. However, while this method often yields accurate results in practice (empirically validated using chemical measurements), theoretical guarantees for convergence and accuracy are lacking. In particular, the DMRG algorithm converges at best to local minima of the energy function. The number of such local minima on the set of low-rank tensors as well as their approximation quality is unknown.

In this paper, our main focus is the minimization of the Rayleigh quotient \[R_H(\psi) \;=\; \frac{\psi^T H \psi}{\psi^T \psi}\] for a Hamiltonian \(H\), where \(\psi\) is constrained to lie in a given tensor train variety. We first study the constrained Rayleigh quotient optimization problem from the perspective of algebraic geometry over arbitrary projective varieties. We define the Rayleigh–Ritz degree (RR degree) of a projective variety in Section 2 as the number of complex critical points of this problem for a generic symmetric matrix \(H\). Here, we present various ways of computing these critical points (Propositions 2 and 7). Proposition 6 ([11]) gives an upper bound on the RR degree. In Section 3, we define the RR correspondence of a projective variety (Theorem 9) along with its parametric analog (Theorem 10). We then define the RR discriminant, which is the set of symmetric matrices \(H\) for which the critical point set does not have RR degree many isolated points. In Section 4, we follow [11] and interpret our optimization problem as a Euclidean distance optimization problem with respect to the Bombieri-Weyl inner product. This leads us to define the BW correspondence. We show that it is rational whenever the original variety is rational. This is our Corollary 4.

Next, we turn to tensor train varieties in Section 5 and identify instances when they are Segre products of projective spaces (Theorem 23). We also report what we know about the defining ideals of tensor trains. A birational parametrization of TT varieties from a product of Grassmannians will be presented in Section 6. From this parametrization, we obtain the formula for the dimension of a TT variety (Corollary 8).

We present our numerical experiments in Section 7. We compute all critical points of our optimization problem for small cases using the numerical algebraic geometry package HomotopyContinuation.jl [12]. This enables us to test the correctness of the ALS and DMRG algorithms on these examples. In small cases, we observe that the ALS method can converge to any of the local minima and that the corresponding energy value can be far from optimal.

2 Rayleigh-Ritz Optimization and the RR degree↩︎

We start with the energy minimization problem arising in quantum chemistry. Given a matrix \(H\in {\rm Sym}^2(\mathbb{R}^{n})\) and a complex projective algebraic variety \(V \subseteq \mathbb{P}^{n-1}\) not contained in the isotropic quadric \(Q = \mathcal{V}(\psi^T\psi)\), we seek to solve the following optimization problem: \[\begin{align} \label{eqn:optproblem} &\textrm{minimize} \,\, \frac{\psi^T H \psi}{\psi^T\psi} \quad \, \textrm{subject to} \,\, \psi \in V_\mathbb{R}. \end{align}\tag{1}\] Here \(V_\mathbb{R}\) denotes the real locus of the complex projective variety \(V\). This is a relaxation of the problem “Find eigenvectors of \(H\) in a variety \(V\)”. This problem is over-constrained, as the set of matrices that have eigenvectors in a fixed variety \(V\) is a closed subvariety of \({\rm Sym}^2(\mathbb{C}^n)\). Such a set is an example of a nonlinear Kalman variety; see [13].

When \(V = \mathbb{P}^{n-1}\), the optimal value for 1 is the smallest eigenvalue of \(H\), and a corresponding eigenvector is an optimal solution. The rational function \(\frac{\psi^T H \psi}{\psi^T\psi}\) is known as the Rayleigh quotient. When \(V\) is a linear space, then 1 is a Rayleigh-Ritz optimization problem over \(V\) [14]. We maintain this terminology in the nonlinear case.

A straightforward way of solving 1 is by computing all of its complex critical points and identifying among them real critical points which minimize the value of the Rayleigh quotient. For us, a critical point is a nonsingular \(\hat{\psi} \in V\) with \(\hat{\psi}^T\hat{\psi} \neq 0\) such that the gradient \(\nabla_\psi\left( \frac{\psi^T H \psi}{\psi^T\psi}\right)\) evaluated at \(\hat{\psi}\) is in \((T_{\hat{\psi}}V)^\perp\).

In Theorem 9, we will prove that the critical ideal \(I_{\mathrm{crit}}(V,H)\) is a zero-dimensional ideal for generic \(H\) and the number of critical points is a constant for all such \(H\). Therefore, we define the Rayleigh-Ritz degree (RR degree) of the variety \(V\) to be the number of complex critical points of the optimization problem 1 for general symmetric matrices \(H\).

The RR degrees of Segre-Veronese varieties were recently derived in [11], where, in particular, explicit formulae for the RR degrees of rank-one matrices of any size and of rank-one binary tensors were obtained.

2.1 The Lagrangian Locus↩︎

We introduce the Lagrangian locus of a projective variety and a symmetric matrix \(H\). This yields another method for computing the critical points of 1 , and gives an upper bound on the RR degree of a variety.

In the above, the critical locus consists of points where the gradient of the level set \({\frac{\psi^T H \psi}{\psi^T \psi} = \lambda}\) is in the span of the gradients of the forms defining \(V\). These level sets correspond to a pencil of quadrics \(\{\mu_1 \psi^T H \psi - \mu_2 \psi^T\psi =0\}\) with \([\mu_1 \, : \, \mu_2] \in \mathbb{P}^1\); see [15]. The gradient vector of this pencil at the point \(\hat{\psi} \in \mathbb{P}^{n-1}\) is the normal vector to a hyperplane \(\Gamma_{\hat{\psi}}\). This way we get a pencil of hyperplanes \(\{\mu_1 (H\hat{\psi})^T\psi - \mu_2 (\hat{\psi})^T \psi = 0\}\). Equivalently, this pencil can be characterized as a family of hyperplanes containing the codimension two linear subspace \(\Lambda_{\hat{\psi}}\) given by \((H\hat{\psi})^T \psi = (\hat{\psi})^T \psi = 0\). Let \(V_{\rm reg} = V\backslash V_{\rm sing}\) be the regular locus of \(V\).

We now specify the genericity condition that ensures that the critical and Lagrangian loci are equal.

For completeness, we include a sketch of the proof leaving out the details about intersection theory; for these, see [16] and [17].

2.2 Parametric Critical Locus↩︎

We now consider the case when \(V \subseteq \mathbb{P}^{n-1}\) admits a birational parametrization \(\psi \, : U \dashrightarrow V\) where \(U\) is either an affine space or a product of projective spaces. Note that this includes tensor train varieties with the parametrization described in Section 6. We consider a parametric version of the Rayleigh-Ritz optimization problem 1 , namely, \[\label{eqn:opt-param} \textrm{minimize } R_H(x) = \frac{\psi(x)^TH\psi(x)}{\psi(x)^T\psi(x)}.\tag{2}\] This is an unconstrained optimization problem, and its critical points are precisely the values of \(x \in \psi^{-1}(V_{\rm reg} \backslash Q)\) for which \(\nabla_x R_H(x) = 0\). If \(x\) is critical for 2 , then its image \(\psi(x)\) is critical for 1 if the Jacobian of the parametrization \(\mathrm{Jac}(\psi)\) defines a surjective map onto the tangent space \(T_{\psi(x)} V\) of \(V\) at \(\psi(x)\); see [18]. We therefore fix the notation \(I_{\psi, s}\) for the ideal generated by \(s\)-minors of the Jacobian \(\mathrm{Jac}(\psi)\).

Similarly, the Lagrangian locus has a parametric analog. As above, let \(\psi: U \dashrightarrow V\) be a birational parametrization of a variety \(V\). Suppose that the parametrization space \(U\) has (projective or affine) coordinates \(x_1, \ldots, x_m\). Write \(f(x)\) for the numerator \(\psi(x)^TH\psi(x)\) of \(R_H(x)\) and \(g(x) = \psi(x)^T\psi(x)\) for the denominator. By the quotient rule, the derivative with respect to \(x_i\) vanishes precisely when \(g(x)\frac{\partial f}{\partial x_i}(x) - f(x)\frac{\partial g}{\partial x_i}(x)\) equals zero. Therefore the critical points of 2 are cut out by the \(2\)-minors of \[\label{eqn:2-row-mat} \begin{pmatrix} f & \frac{\partial f}{\partial x_1} & \cdots & \frac{\partial f}{\partial x_m}\\ g & \frac{\partial g}{\partial x_1} & \cdots & \frac{\partial g}{\partial x_m} \end{pmatrix} = \begin{pmatrix} \psi(x)^T\\ \psi(x)^TH \end{pmatrix} \cdot \begin{pmatrix} \psi(x) & \mathrm{Jac}(\psi)(x) \end{pmatrix}.\tag{3}\] We remark that if \(U\) is a projective space, then the first column of 3 can be omitted using Euler’s relation.

As in the implicit case, we obtain the critical locus after saturating by isotropic quadric.

3 Rayleigh-Ritz Correspondence and Discriminant↩︎

In this section, we allow the symmetric data matrix \(H\) to vary. Because the Rayleigh quotient \(\frac{\psi^TH\psi}{\psi^T\psi}\) is linear in the entries of \(H\), the problem 1 has the same critical points if \(H\) is replaced with any nonzero scalar multiple \(\lambda H\). We therefore only consider the data matrices up to scaling, and accordingly work in \(\mathbb{P}({\rm Sym}^2(\mathbb{C}^n))\). We first define the Rayleigh-Ritz correspondence as the incidence variety of symmetric matrices \(H\) and the corresponding critical points \(\psi\). We define this correspondence for general projective varieties. We illustrate this construction with small examples, namely projective spaces, rational normal curves and the smallest tensor train variety, \(\mathbb{P}^1 \times \mathbb{P}^1\). Similar correspondences have been introduced for maximum likelihood estimation [19], [20] and for Euclidean distance optimization [21].

We then define the Rayleigh-Ritz discriminant, which describes the matrices for which the critical locus of 1 exhibits nongeneral behavior. Analogous discriminants were defined in [22] for the maximum likelihood estimation problem and in [21] for the Euclidean distance optimization problem.

3.1 Rayleigh-Ritz Correspondence↩︎

Throughout this section, \(V \subseteq \mathbb{P}^{n-1}\) represents an irreducible projective variety not contained in the isotropic quadric \(Q = V(\psi^T\psi)\). We define the Rayleigh-Ritz correspondence (RR correspondence) as the pairs of smooth points \(\psi \in V_{\rm reg}\) and symmetric matrices \(H\) such that \(\psi\) is a critical point of 1 for that matrix \(H\): \[\label{eqn:RR-cor} \mathcal{R}_V = \overline{\{(\psi, H) \in V_{\rm reg} \times \mathbb{P}({\rm Sym}^2(\mathbb{C}^n)): \psi \textrm{ is a critical point of \eqref{eqn:optproblem}} \}}.\tag{4}\] Here both the vector \(\psi\) and the data matrix \(H\) are indeterminates. If \(V\) has codimension \(c\), then by Proposition 4 and Theorem 5 the correspondence is the zero locus of the ideal \[\label{eqn:IRV} \big \langle (c + 2)\textrm{-minors of \eqref{eqn:lagrange-mat}} \big \rangle : (I_{V_{\rm sing}} \cdot \langle \psi^T\psi \rangle)^\infty.\tag{5}\] If \(V\) is smooth and intersects the quadric \(Q\) transversally, then this ideal can be computed without saturation by Theorem 5.

We turn now to the parametric situation introduced in Section 2. Suppose \(V\) admits a birational parametrization \(\psi: U \dashrightarrow V\) from an affine space or product of projective spaces \(U\). We define the parametric RR correspondence as \[\mathcal{P}_V = \overline{\{(x, H) \in \psi^{-1}(V_{\rm reg} \backslash Q) \times \mathbb{P}({\rm Sym}^2(\mathbb{C}^n)) : x \textrm{ is a critical point of \eqref{eqn:opt-param}} \}}.\] Since the parametric and implicit critical points agree, \[\psi \times {\rm Id}: U \times \mathbb{P}( {\rm Sym}^2(\mathbb{C}^n) )\dashrightarrow \mathbb{P}^{n-1} \times \mathbb{P}( {\rm Sym}^2(\mathbb{C}^n) )\] is a birational map from the parametric correspondence \(\mathcal{P}_V\) to the usual correspondence \(\mathcal{R}_V\). By Corollary 2, the parametric correspondence is the zero set of the ideal \[\big \langle 2\textrm{-minors of } {\rm Jac}(\psi)^T \cdot \begin{pmatrix} \psi(x) & H\cdot \psi(x) \end{pmatrix} \big \rangle : (\langle (n-c)\textrm{-minors of } {\rm Jac}(\psi) \rangle \cdot \langle \psi(x)^T\psi(x)\rangle \cdot I_B)^\infty\] where \(I_B\) is the ideal of the base locus \(B\) of \(\psi\).

3.2 Rayleigh-Ritz Discriminant↩︎

In Theorem 9, we proved that for a general choice of data matrix \(H\), the fiber \(\pi_2^{-1}(H)\) has RR degree many points. By “general,” we mean that this statement holds outside of a subvariety of the data space \(\mathbb{P}({\rm Sym}^2(\mathbb{C}^n))\). We call this subvariety, which is typically a hypersurface, the Rayleigh-Ritz discriminant (RR discriminant): \[\Sigma_V = \overline{\{ H \in \mathbb{P}({\rm Sym}^2(\mathbb{C}^n)) : \textrm{ the critical locus of \eqref{eqn:optproblem} is infinite or nonreduced} \}}.\]

We use the Jacobian criterion to determine when the critical locus of a matrix \(H\) is singular by treating the variables \(h_{11}, h_{12}, \ldots, h_{nn}\) as parameters. Given the ideal \(I_{\mathcal{R}_V} = \langle g_1, \ldots, g_\ell\rangle \subseteq \mathbb{C}[\psi_1, \ldots,\psi_n, h_{11}, h_{12}, \ldots, h_{nn}]\) of the correspondence, a point in the critical locus is nonreduced or lies on a positive dimensional component precisely when the following \(\ell \times n\) matrix drops rank: \[\label{eqn:JpsiRV} {\rm Jac}_\psi(\mathcal{R}_V) = \begin{pmatrix} \frac{\partial g_1}{\partial \psi_1} & \frac{\partial g_1}{\partial \psi_2} & \cdots & \frac{\partial g_1}{\partial \psi_n}\\ \frac{\partial g_2}{\partial \psi_1} & \frac{\partial g_2}{\partial \psi_2} & \cdots & \frac{\partial g_2}{\partial \psi_n}\\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial g_\ell}{\partial \psi_1} & \frac{\partial g_\ell}{\partial \psi_2} & \cdots & \frac{\partial g_\ell}{\partial \psi_n}\\ \end{pmatrix}.\tag{6}\] The RR discriminant can therefore be written as \[\Sigma_V = \pi_2(\{ (\psi, H) \in \mathcal{R}_V : {\rm rank}\,\, {\rm Jac}_\psi(\mathcal{R}_V)(\psi, H) \leq n-2 \}).\]

The RR discriminant can also be computed from the parametric RR correspondence. Suppose that the parametrization space \(U\) has (projective or affine) coordinates \(x_1, \ldots, x_m\). If \(I_{\mathcal{P}_V} = \langle f_1, \ldots, f_p \rangle \subseteq \mathbb{C}[x_1, \ldots, x_m, h_{11}, h_{12}, \ldots, h_{nn}]\), we define the \(p \times m\) Jacobian \[\label{eqn:JxPV} {\rm Jac}_x(\mathcal{P}_V) = \begin{pmatrix} \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_m}\\ \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots & \frac{\partial f_2}{\partial x_m}\\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial f_p}{\partial x_1} & \frac{\partial f_p}{\partial x_2} & \cdots & \frac{\partial f_p}{\partial x_m}\\ \end{pmatrix}.\tag{7}\]

We introduced \(\Sigma_V\) as the set of matrices \(H\) whose critical locus is infinite or singular. This is the perspective taken in [22] when defining the discriminant of maximum likelihood estimation. The authors of [21] take a different perspective, defining the discriminant of the Euclidean distance optimization problem as the branch locus of the projection \(\pi_2\) onto the second factor. We now argue that these two definitions agree for RR discriminants.

We now examine the RR discriminants of projective spaces and rational normal curves, before turning to the components of RR discriminants.

We now turn to the structure of the discriminant. The RR discriminant is generally not irreducible. Let \(\mathcal{R}^{\rm ram}_V = \{ (\psi, H) \in \mathcal{R}_V : {\rm rank}\,\, J_\psi(\mathcal{R}_V) \leq n-2\}\) denote the ramification locus of the projection \(\pi_2\). We can write \(\mathcal{R}^{\rm ram}_V\) as the union \[\mathcal{R}^{\rm ram}_V = (\mathcal{R}^{\rm ram}_V \cap Q) \cup \overline{(\mathcal{R}^{\rm ram}_V \setminus Q)}\] where \(Q\) denotes \(\mathcal{V}(\psi^T\psi) \subseteq \mathbb{P}^{n-1} \times \mathbb{P}({\rm Sym}(\mathbb{C}^n))\). The projections of these components \(\pi_2(\mathcal{R}^{\rm ram}_V \cap Q)\) and \(\pi_2(\overline{(\mathcal{R}^{\rm ram}_V \setminus Q)})\) are often, but not always, hypersurfaces. We call the set \(\pi_2(\mathcal{R}^{\rm ram}_V \cap Q)\) the isotropic part of the discriminant and we call \(\pi_2(\overline{(\mathcal{R}^{\rm ram}_V \setminus Q)})\) the nonisotropic part of the discriminant. We illustrate this phenomenon with some small examples.

In the following examples, we were not able to compute the full discriminant polynomials. The degrees were computed symbolically by intersecting the discriminant hypersurface with a line. We can give a lower bound on the number of regions in the discriminant complement by computing all critical points for 2 for random real matrices and counting the numbers of real solutions that arise.

We have verified Conjecture 16 numerically for \(d \leq 10\); see [23] for computations.

It is not clear whether there can exist \((\psi, H) \in \mathcal{R}^{\rm ram}_V \cap Q\) with real \(H\). The vector \(\psi\) cannot be real, since \(\psi^T\psi = 0\) is a sum of squares. In the examples we have computed, the variety \(\pi_2(\mathcal{R}^{\rm ram}_V \cap Q)\) has no real points and its defining polynomial factors as a sum of squares over \(\mathbb{Q}\). We therefore make the following conjecture, which we have verified for the rational quadric (Example 13), the twisted cubic (Example 14), and projective space (Example 12).

One motivation for studying the RR discriminant of a variety is to compute the average RR degree, that is, the expected number of real critical points. The role of the discriminant in this computation is explained in the Euclidean distance case in [21]. In particular, in this application we are only interested in the real locus of the discriminant. Thus, if Conjecture 17 holds, then only the nonisotropic part of the discriminant is relevant. We therefore present degrees of the nonisotropic parts of RR discriminants of various varieties in Tables 2 and 3, and make the following conjecture.

4 Bombieri-Weyl Correspondence↩︎

In this section, we will introduce yet another correspondence that essentially parametrizes pairs \((\psi, H) \in V \times \mathrm{Sym}^2(\mathbb{C}^n)\) where \(\psi\) is a critical point of the Rayleigh-Ritz optimization problem 1 . For this, we will consider \(V\) in its second Veronese embedding \(X = \nu_2(V)\) and transform the Rayleigh-Ritz optimization problem 1 to a distance optimization problem using Bombieri-Weyl distance. The equivalence of these two optimization problems has been introduced and fully developed in [11]. We give a very brief introduction to the Bombieri-Weyl distance optimization and cite the main result we will use. This optimization problem is a variant of the Euclidean distance optimization problem. We will show that, as in the Euclidean distance case, the Bombieri-Weyl correspondence (BW correspondence) of a parametrized variety is itself parametrized. We include this with the hope that it will be useful for computing the average number of real critical points of 1 ; see [21].

First, note that the Rayleigh-Ritz optimization problem 1 is equivalent to the following optimization problem \[\begin{align} \label{eqn:optsphere} &\textrm{minimize} \,\, f_H(\psi) = \psi^T H \psi \quad \, \textrm{subject to} \,\, \psi \in C(V_\mathbb{R}) \cap S^{n-1} \end{align}\tag{8}\] where \(C(V_\mathbb{R}) \subset \mathbb{R}^n\) is the cone over the projective variety \(V \subset \mathbb{P}^{n-1}\) and \(S^{n-1}\) is the \((n-1)\)-dimensional sphere in \(\mathbb{R}^n\). We also point out that each critical point \(\hat{\psi}\) of 1 corresponds to two critical points \(\pm \frac{1}{(\hat{\psi})^T \hat{\psi}} \hat{\psi}\) of 8 .

Recall that any homogeneous form \(f\) of degree \(d\) in \(n\) variables with real coefficients can be represented by an element of \(\mathop{\mathrm{Sym}}^d(\mathbb{R}^n)\). This representation as a symmetric tensor can be written as \[f(x) = \sum_{|\alpha| = d} \binom{d}{\alpha}f_\alpha x^\alpha,\] where \(\binom{d}{\alpha}\) is the multinomial coefficient \(\frac{d!}{\alpha_1 ! \cdots \alpha_n !}\). Hence \(f\) is identified with \((f_\alpha)_{|\alpha| = d}\). The case we will work with is \(d=2\). A homogeneous quadratic polynomial in \(n\) variables and the corresponding symmetric matrix of this quadratic form are \[\label{eq:quadratic} \sum_{i=1}^n h_{ii} x_i^2 + \sum_{i<j} 2 h_{ij}x_ix_j, \quad \quad \quad H = \begin{footnotesize} \begin{pmatrix} h_{11} & h_{12} & \cdots & h_{1n} \\ h_{12} & h_{22} & \cdots & h_{2n} \\ \vdots & & \ddots & \vdots \\ h_{1n} & h_{2n} & \cdots & h_{nn} \end{pmatrix} \end{footnotesize}.\tag{9}\] From now on we will identify a symmetric matrix \(H\) with the quadratic form \(f_H = x^T H x\).

Let \(I_X = \langle g_1, \ldots, g_s\rangle\) be the prime ideal in \(\mathbb{R}[Z] = \mathbb{R}[z_{ij} \, : \, 1\leq i \leq j \leq n]\) defining \(X\). Suppose the codimension of \(I_X\) is \(c\). We let \(\mathrm{Jac}(I_X)\) to be the Jacobian of the ideal \(I_X\). Then the critical ideal is \[\label{eqn:critical-ideal} \left( I_X + \bigg \langle (c+1) -minors of\left( \begin{array}{c} \nabla_Z(\mathrm{trace}((Z-H)^2)) \\ \mathrm{Jac}(I_X) \end{array} \right)\bigg \rangle \right) \, : \, I_{X_{\rm sing}}^\infty.\tag{10}\]

5 Geometry of Tensor Train Varieties↩︎

In this section, we introduce tensor train varieties (TT varieties) – projective varieties of tensors in a tensor train format. These varieties can be described by the combinatorial structure of a graph following the tensor network viewpoint as, for example, in [6], [30], [31].

Let \(\Gamma_n\) be an undirected path graph on vertices \(\{0,1,\dots,n\}\), with edges \(\{i-1,i\}\) for \(i=1,\dots,n\). We assign positive integers \({\mathbf{k}} = (k_0, k_1, \ldots, k_n)\) to the vertices, known as physical dimensions, and \({\mathbf{r}} = (r_1, r_2, \ldots, r_n)\) to the edges, called bond dimensions; by convention, we set \(r_0=r_{n+1}=1\). We identify tensors in the ambient projective space \(\mathbb{P}\big(\bigotimes_{i=0}^n \mathbb{C}^{k_i}\big)\) with points having homogeneous coordinates \(\psi_{j_0,\dots,j_n}\), where \(j_i\in[k_i]:=\{1,\dots,k_i\}\). For each index \(i\in\{0,1,\ldots,n\}\) and each \(j_i\in[k_i]\), let \(A_i^{(j_i)}\) be a general \(r_i\times r_{i+1}\) matrix of parameters.

It is often convenient to characterize \(V_{\mathbf{k},\mathbf{r}}\) by rank constraints on certain matrix flattenings of a tensor. For \(i \in \{0,1,\dots,n-1\}\), let \(\psi^i\) denote the flattening of \(\psi\) according to the partition \(\{0,1,\dots,i\}\sqcup\{i+1,\dots,n\}\), viewed as a \((k_0\cdots k_i)\times (k_{i+1}\cdots k_n)\) matrix.

Example 17. Take \(\mathbf{k}=(2,2,2,2)\) and \(\mathbf{r}=(1,2,1)\). The parametrization 12 yields \[\psi_{ijkl}=A_0^{(i)}A_1^{(j)}A_2^{(k)}A_3^{(\ell)}, \qquad i,j,k,\ell\in\{1,2\},\] where \(A_0^{(i)}\) and \(A_3^{(\ell)}\) are scalars, \(A_1^{(j)}\in\mathbb{C}^{1\times 2}\), and \(A_2^{(k)}\in\mathbb{C}^{2\times 1}\). Writing \[A_0^{(i)}=(a_0^{i}),\quad A_1^{(j)}= \begin{pmatrix} a_{11}^{j} & a_{12}^{j}\end{pmatrix},\quad A_2^{(k)}=\begin{pmatrix} a_{21}^{k} \\ a_{22}^{k}\end{pmatrix},\quad A_3^{(\ell)}=(a_3^{\ell}),\] one checks that \(V_{\mathbf{k},\mathbf{r}}\) is isomorphic to the Segre product \(\mathbb{P}^3\times (\mathbb{P}^1)^2\) whose vanishing ideal is minimally generated by the \(2\)-minors of the flattenings \(\psi^0\) and \(\psi^2\). Minors of \(\psi^1\) don’t appear since \(r_2 = 2 = 1\cdot2 = r_{1}k_{1}\). The corresponding tensor network is shown in Figure 1.

From now on, we use the shorthand \((1)_n := (\underbrace{1,\dots,1}_{n})\) for sequences of repeated indices.

We next isolate a large, explicit family of tensor train varieties that are again Segre products. We say that indices \(i\) with \(r_i>1\) are separated if no two such indices are adjacent; likewise, we call a separated block a consecutive sequence \((i_1,\dots,i_t)\) with \({r_{i_1},\dots,r_{i_t}>1}\) such that no other \(r_j>1\) is adjacent to the block.

We are mostly interested in binary tensors (\(k_i=2\)) due to their prevalence in applications.

The dimensions of the varieties above also follow from the general formula in Corollary 8.

We now turn to equations. Lemma 2 shows that \(V_{\mathbf{k},\mathbf{r}}\) is set-theoretically cut out by determinantal equations. The following conjecture, first suggested to us by Bernd Sturmfels in personal communication, asserts that these suffice scheme-theoretically as well.

We collect some positive evidence in the case of binary tensors, in particular, in the case of Segre varieties whose Gröbner bases have been studied widely. We start with the following lemma; a proof is available in [32].

We record a conjecture that represents a key step towards the proof of the Conjecture 24.

6 Birational Parametrization and TT manifolds↩︎

In this section, we present a one-to-one parametrization of a dense open set of the tensor train variety \(V_{\mathbf{k}, \mathbf{r}}\) from a product of Grassmannians defined by \[\mathcal{Z} \coloneq \mathrm{Gr}(r_1,k_0) \times \mathrm{Gr}(r_2,r_1 k_1) \times \cdots \times \mathrm{Gr}(r_{n},r_{n-1} k_{n-1}) \times \mathbb{P}(\mathbb{C}^{r_{n} \times k_n}).\] This open set is contained in the tensor train manifold, denoted by \(V_{\mathbf{k},\mathbf{r}}^=\), which is the set of tensors in \(V_{\mathbf{k},\mathbf{r}}\) whose \(i\)-th flattening has rank exactly \(r_{i+1}\) for all \(i = 0, \ldots, n-1\). The construction is inspired by the TT-cross algorithm [34], itself an adaptation of the TT-SVD [9], [35] in which successive SVDs are replaced by so‑called skeleton decompositions of matrices [36]. Skeleton (or \(CUR\)) decompositions factor an \(m\times n\) matrix \(A\) by selecting \(r\) linearly independent columns and rows of \(A\) as columns and rows of \(C \in \mathbb{C}^{m \times r}\) and \(R \in \mathbb{C}^{r \times n}\), respectively. In the literature, this decomposition is often generalized to the case where only rows or columns are chosen, resulting in an \(XR\)-decomposition or \(CX\)-decomposition, respectively [37].

In our algorithm, we use a special case of the \(XR\)-decomposition for a rank-\(r\) matrix \(A\), under the following assumption.

If Assumption 1 holds, the \(XR\)‑decomposition is defined as: \[(X,R) \coloneq \mathtt{XR}(A,r),\] where \(R \in \mathbb{C}^{r \times n}\) is formed from the first \(r\) rows of \(A\), and \(X = A R^\dagger (R R^\dagger)^{-1} \in \mathbb{C}^{m \times r}\). Invertibility of \(RR^\dagger\) is guaranteed by Assumption 1. This decomposition satisfies \(XR = AR^\dagger(RR^\dagger)^{-1}R = AP\), where \(P\) is the orthogonal projection onto the rowspan of \(A\). Since all rows of \(A\) are in the rowspan of \(A\), projecting onto the rowspan of \(A\) leaves them unchanged, and hence \(XR = AP = A\). Moreover, \(X\) has the form \[\label{eq:Xgr} X = \begin{bmatrix} {\rm Id}_r \\ B \end{bmatrix},\tag{15}\] with \(B \in \mathbb{C}^{(m-r) \times r}\). Such matrices parametrize a large open subset of \(\mathrm{Gr}(r,m)\) via their column spans, with the first block fixed to the identity to ensure a one‑to‑one parametrization.

Our algorithm gives a bijection from an open set \(\mathcal{U} \subseteq V_{\mathbf{k}, \mathbf{r}}^=\) to an open set \(\mathcal{W} \subseteq \mathcal{Z}\). In order to write down our algorithm, we will first need to identify these open subsets. The description of the open set \(\mathcal{W}\) is straightforward. Let \(W_i\) denote the open set \[W_i := \{\mathrm{span}(A) : A \in \mathbb{C}^{r_{i-1}k_{i-1} \times r_i},\;\text{first r_i rows linearly independent}\} \subseteq \mathrm{Gr}(r_i,r_{i-1}k_{i-1})\] for \(i = 1,\ldots,n\) and let \(\mathbb{C}_\ast^{r_{n} \times k_n} \subseteq \mathbb{C}^{r_{n} \times k_n}\) denote the open set of full rank matrices. Then the following set is open in \(\mathcal{Z}\): \[\mathcal{W} \coloneq W_1 \times \cdots \times W_{n} \times \mathbb{C}_\ast^{r_{n} \times k_n}.\] The set \(\mathcal{U}\) consists of tensors for which Assumption 1 holds in every step of the algorithm; its precise description appears in the proof of Theorem 28, where we show it is Zariski dense.

Similar to Section 5, we denote the \(i\)-th flattening of a tensor \(\psi \in \bigotimes_{i=0}^n \mathbb{C}^{k_i}\) by \[\psi^i = \mathtt{unfold}(\psi,[k_0 \cdots k_i,k_{i+1} \cdots k_n]) \in \mathbb{C}^{(k_0 \cdots k_i) \times (k_{i+1} \cdots k_n)}\] and the tensorization of a matrix \(\psi^i \in \mathbb{C}^{(k_0 \cdots k_i) \times (k_{i+1} \cdots k_n)}\) is denoted by \[\psi = \mathtt{tensorize}(\psi^i,[k_0,\ldots,k_n]) \in \bigotimes_{i=0}^n \mathbb{C}^{k_i}.\] With this, we are equipped to present the parametrization algorithm for a tensor with given TT-rank \(\mathbf{r} = (r_1,\ldots,r_{n})\). The procedure is similar to the one in Lemma 2; see Algorithm 2.

We now describe the inverse of Algorithm 2, namely the decompression procedure that reconstructs a TT tensor from its parametrized form; see Algorithm 3.

The dimension of the manifold \(V_{\mathbf{k},\mathbf{r}}^=\) is known in the tensor decomposition literature; see, for example, [39]. It can also be derived from [40], where more general tensor network varieties were studied. We present it here as a corollary of Theorem 28.

7 Numerical Experiments↩︎

In our last section, we return to TT varieties and perform numerical computations for some small instances. The primary numerical task is the computation of all complex critical points of 1 over a TT variety. For this, we use homotopy continuation and monodromy methods. We report the results of our computations for small determinantal varieties as well as for the first few instances of TT varieties that are Segre products by Theorem 23; see Tables 1 and 2. In many of these cases, the RR degrees themselves can be derived from the formulae obtained in [11]; see, in particular, Corollaries 5.2 and 5.3. However, we compute not only the degrees, but the actual critical points. One of our observations is that the number of local extrema is typically much smaller than the corresponding RR degree.

We also report the degree of the nonisotropic part of the RR discriminant (see Section 3.2) for some determinantal and TT varieties; see Tables 2 and 3. The take-away message there is that these degrees are high, implying that the computation of the actual discriminant polynomials are currently out of reach.

Finally, we try to gauge the quality of the low energy state solutions produced by the DMRG and ALS algorithms by comparing to the complete list of local and global minima produced by our numerical methods. We conclude that ALS frequently gets stuck in local minima with suboptimal energy, although this happens less often for physical Hamiltonians; see Table 4. The DMRG method does not converge to critical points of 1 at all.

7.1 Hamiltonians↩︎

In our experiments, whenever we look at a specific Hamiltonian \(H\), we consider two types:

1. Random symmetric real matrices.

   We generate a symmetric real matrix \(H \in \mathbb{R}^{(k_0 \cdots k_n) \times (k_0 \cdots k_n)}\) by drawing entries from the standard normal distribution and symmetrizing. For the ALS and DMRG experiments, \(H\) is then expressed in tensor operator form, see [9].

2. Random Hamiltonians in second quantization.

   In second-quantized form, the Schrödinger operator has the structure (see [6]): \[H = \sum_{i,j} t_{ij} \mathbf{a}_i^\dagger \mathbf{a}_j + \sum_{i,j,k,l} v_{ijkl} \mathbf{a}_i^\dagger \mathbf{a}_j^\dagger \mathbf{a}_k \mathbf{a}_l,\] where \(\mathbf{a}_i^\dagger = \mathbf{s} \otimes \cdots \otimes \mathbf{s} \otimes \mathbf{a}^\dagger \otimes \mathbf{I} \otimes \cdots \otimes \mathbf{I}\) and \(\mathbf{a}_i = \mathbf{s} \otimes \cdots \otimes \mathbf{s} \otimes \mathbf{a} \otimes \mathbf{I} \otimes \cdots \otimes \mathbf{I}\) are the creation and annihilation operators, respectively. Here, we have \[\mathbf{s} = \begin{footnotesize}\begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \end{footnotesize}, \qquad \mathbf{a} = \begin{footnotesize}\begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} \end{footnotesize}, \qquad \mathbf{I} = \begin{footnotesize}\begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \end{footnotesize}.\] The coefficients \(t_{ij}, v_{ijkl} \in \mathbb{R}\) for \(i,j,k,l \in \{0,1,\dots,n\}\) are chosen randomly from the standard normal distribution. In this case, \(\mathbf{k} = (2)_{n+1}\), and \(H\) is symmetric.

7.2 Homotopy Continuation Computations↩︎

A central task in our setting is to solve polynomial systems that encode the optimization problem 2 on tensor trains. There are many ways to write down such a system, for example by using Proposition 2, Proposition 7, or 10 .

Our main computational tool is numerical homotopy continuation, see [41]. This is a numerical method for finding all isolated solutions to a system of polynomial equations. Concretely, let \(\mathbf{x}=(x_1,\dots,x_n)\) be variables and \(\mathbf{p}=(p_1,\dots,p_k)\) parameters, and consider a family of polynomial systems \[\mathcal{F}(\mathbf{x};\mathbf{p})=\big(f_1(\mathbf{x};\mathbf{p}),\ldots,f_n(\mathbf{x};\mathbf{p})\big)\in\mathbb{C}[\mathbf{x}]^n,\quad \text{ with } \; \mathbf{p}\in\mathbb{C}^k.\] Choose two parameter values \(\mathbf{q}_1,\mathbf{q}_2\in\mathbb{C}^k\) such that the target system is \(F(\mathbf{x})=\mathcal{F}(\mathbf{x};\mathbf{q}_1)\), and the start system \(G(\mathbf{x})=\mathcal{F}(\mathbf{x};\mathbf{q}_2)\) has easily computable solutions. The parameter homotopy \[H(\mathbf{x},t):=\mathcal{F}\big(\mathbf{x};(1-t)\mathbf{q}_1+t\mathbf{q}_2\big),\quad t\in[0,1],\] deforms \(G(\mathbf{x})\) at \(t=1\) to \(F(\mathbf{x})\) at \(t=0\). The parameter continuation theorem [42], [43] guarantees that the number of regular isolated solutions in \(\mathbb{C}^n\) is constant along \(t\in(0,1]\). Starting from solutions of \(G(\mathbf{x})\), we track solution paths to \(t=0\) using predictor–corrector schemes that numerically integrate the Davidenko ODE \[\big(\partial_{\mathbf{x}} H(\mathbf{x},t)\big)\,\dot{\mathbf{x}} + \partial_t H(\mathbf{x},t)=0.\] If the start system \(G(\mathbf{x})\) has at least as many regular isolated solutions as the target \(F(\mathbf{x})\) and all solution paths arrive at \(t=0\), then we obtain all solutions of \(F(\mathbf{x})\). When \(G(\mathbf{x})\) has exactly as many solutions as \(F(\mathbf{x})\), the homotopy is called optimal.

Another numerical method for solving polynomial systems is monodromy, introduced in [44]. Starting from a general parameter value \(\mathbf{p}_0\) and a set of known solutions to \({\mathcal{F}(\mathbf{x};\mathbf{p}_0)=0}\), one traces each solution along closed loops in the parameter space based at \(\mathbf{p}_0\) using parameter homotopy. These loops induce permutations of the solution set, and running random loops often uncovers previously unknown solutions. When the total number of solutions to a general system in the family is known, monodromy has a stopping criterion. Otherwise, suitable heuristics may be applied. Both methods are implemented in the Julia [45] package HomotopyContinuation.jl [12].

We combine the two homotopy continuation methods introduced above to solve the optimization problem 2 , in particular to compute RR degrees and to numerically determine the degrees of RR discriminants. All algorithms are implemented in our Julia package TensorTrainOptimization.jl using HomotopyContinuation.jl and Oscar.jl [46], which is provided together with the complete dataset of our experiments at:

get_RRdegree_monodromy([2,2], [1])

TT_optimization_monodromy([2,2,2,2], [1,2,1]; H = M1, RRdeg = 352)   

Application-driven Hamiltonians are typically very sparse, with many zero entries. Specializing parameters within the family of systems can only decrease the number of complex solutions, so the number of complex critical points for physical Hamiltonians is typically smaller than the RR degree.

7.3 DMRG and ALS↩︎

In practice, the ground state energy of the Hamiltonian \(H\) is computed by solving 1 using the DMRG method [6]. For its derivation, we note that Algorithm 3 can be naturally generalized to take any tensor train \((A_0,\ldots,A_n) \in \prod_{i=0}^n \mathbb{C}^{r_i k_i \times r_{i+1}}\) as an input, where \(A_i\) is the collection of the \((r_i\times r_{i+1})\) parameter matrices \(A_i^{(j_i)}\). We denote this generalized map by \[\tau : \prod_{i=0}^n \mathbb{C}^{r_i k_i \times r_{i+1}} \rightarrow \mathbb{C}^{k_0 \times \cdots \times k_n}.\] The DMRG is a two-site version of the ALS [10], which we introduce first. A random TT tensor is chosen as an initial guess. In each outer iteration loop (often called a sweep), every component \(A_i\) is updated independently and successively, starting at \(A_0\) and ending at \(A_n\). Typically, the sweep is then performed backwards (from \(A_n\) to \(A_0\)). These forward and backward sweeps are repeated until convergence.

In the \(i\)-th step of each sweep, only the component \(A_i\) is updated. The algorithm ensures that the left components \(A_0,\ldots,A_{i-1}\) are left-unitary, that is, \(A_j^\dagger A_j = I_{r_{j+1}}\) for \(j = 0,\ldots,i-1\), and the right components \(A_{i+1},\ldots,A_n\) are right-unitary: \(B_jB_j^\dagger = I_{r_j}\), where \(B_j = \mathtt{unfold}(A_j,[r_j,k_jr_{j+1}])\) for \(j = i+1,\ldots,n\). The solution of the subproblem \[\label{eqn:ALSsubproblem} \tilde{A}_i = \operatorname*{argmin}_{\substack{Z \in \mathbb{C}^{r_i k_i \times r_{i+1}}}} \frac{ \langle \tau(A_0,\ldots,Z,\ldots,A_n),\, H\,\tau(A_0,\ldots,Z,\ldots,A_n) \rangle }{ \langle \tau(A_0,\ldots,Z,\ldots,A_n),\, \tau(A_0,\ldots,Z,\ldots,A_n) \rangle }.\tag{16}\] can be found by computing the eigenvectors of a significantly reduced Hamiltonian using conventional methods. In a forward sweep, \(\tilde{A}_i\) is subsequently left-orthogonalized (for example, via a QR decomposition, setting \(A_i \leftarrow Q\)) to prepare for the update of \(A_{i+1}\).

Since the initial guess determines the size of each component, the ALS algorithm is automatically constrained to the TT manifold \(V_{\mathbf{k}, \mathbf{r}}^=\) and the TT rank remains fixed. In contrast, the DMRG method is rank adaptive. The initial guess is typically chosen to have TT rank \(\boldsymbol{1}\), i.e., \(r_0 = \cdots = r_n = 1\). Here, all but two neighboring components (say \(A_i\) and \(A_{i+1}\)) are fixed and orthogonalized as above. The subproblem 16 is then solved for these two components simultaneously and without imposing the any constraint on \(r_{i+1}\). This yields a merged double component \(\tilde{A}_{i,i+1} \in \mathbb{C}^{r_ik_i \times k_{i+1}r_{i+2}}\) that is subsequently decomposed to produce new \(A_i\) (left-orthogonal) and \(A_{i+1}\). This is typically done using a truncated SVD. The updated rank \(r_{i+1}\) is chosen to meet prescribed error tolerances.

In theory, for a general initial guess and without rank restrictions, the DMRG method will eventually find the global optimum of the unconstrained Rayleigh quotient. However, for large systems, a maximum rank must be chosen to maintain computational feasibility. Once this maximum rank is reached, the DMRG suffers from the same issues as ALS: it might get stuck in local minima whose quality (in terms of the resulting energy) is unknown.

In our experiments, we have observed another issue with the DMRG. Since we want to compare the output of the DMRG with critical points on the manifold \(V_{\mathbf{k}, \mathbf{r}}^=\), we must truncate the ranks to the maximum rank \(\mathbf{r}\) in each step. While the solution to the subproblem is optimal in terms of energy, the rank truncation of the solution is not. The truncated SVD yields the optimal approximation in the Frobenius norm but not in energy. Hence, the DMRG does not converge to critical points of the Rayleigh quotient on the manifold, but to other points with suboptimal energy. To address this comparison, we ran DMRG 1000 times for five random symmetric matrices \(D_1,\ldots,D_5\) and compared the resulting energies with the local minima computed with HomotopyContinuation.jl. The DMRG analog of Table 4 can be found under the link [23]. The results are summarized in Table 5.

It is possible to replace the truncated SVD with a more advanced rank truncation that finds the best low-rank approximation in terms of the energy in each step (see [47]), but this exceeds the scope of this article.

We are grateful to Bernd Sturmfels for introducing the authors. We thank Flavio Salizzoni, Luca Sodomaco, and Julian Weigert for many helpful conversations and for sharing their notes on RR degrees with us. We thank Otto Schmidt for interesting discussions and for providing several physical Hamiltonians for our experiments. Finally, we thank anonymous MEGA referees for their helpful comments, which helped to improve the presentation of the paper.

VB acknowledges support from the European Research Council

\(\!\!\) Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

Authors’ addresses:

Max Planck Institute for Mathematics in the Sciences Leipzig, Germany borovik@mis.mpg.de

University of California Berkeley, USA hannahfriedman@berkeley.edu

San Francisco State University, USA serkan@sfsu.edu

Georg-August-Universität Göttingen, Germany m.pfeffer@math.uni-goettingen.de

References↩︎