The 4 $\times$ 4 orthostochastic variety

A real square matrix is called orthostochastic if it is the entrywise square of an orthogonal matrix. Explicitly, $A=(a_{ij})\in{\mathbb{R}}^{n\times n}$ is orthostochastic if there exists an orthogonal matrix $V=(v_{ij})$ with $a_{ij}=v_{ij}^{2}$ for all $i,j=1,\ldots,n$. It follows immediately that an orthostochastic matrix is doubly stochastic (i.e. has nonnegative entries and all row and column sums equal to 1), as all rows and columns of an orthogonal matrix are unit vectors.

As an interesting and special class of doubly stochastic matrices, orthostochastic (and their unitary generalization, unistochastic) matrices arise naturally in a number of contexts, including spectral theory [13, 14], convex analysis [1, 11], and physics [3, 7]. More recently – and of interest in algebraic geometry – it has been shown that orthostochastic matrices are deeply connected to definite symmetric determinantal representations of real polynomials. Indeed, for hyperbolic plane curves, every monic symmetric determinantal representation arises from certain associated orthostochastic matrices, which yields an effective algorithm [5, 8] for computing such representations for cubic curves.

It is thus of interest to find intrinsic characterizations of the set of orthostochastic matrices. One approach to this is: by definition, the set of orthostochastic matrices is the image of an algebraic variety (the real orthogonal group) under a polynomial map (coordinate-wise squaring) – thus the Zariski closure of the image is a real algebraic variety. Our goal then is to find equations for this Zariski closure, which we refer to as finding equations for the orthostochastic variety. The equations of the $3\times 3$ orthostochastic variety are known, which made the computation of determinantal representations of cubic curves possible, but for matrices of size $\geq 4$ no set of equations which define the orthostochastic variety were previously known. In view of this, our main result is:

We outline the remainder of the article: in Section 2 we formalize the set-up and notation, and review some basic invariants of the varieties under consideration. In Section 3 we give a procedure for obtaining a naive set of equations which always cut out a superset of the orthostochastic variety, and see how dimension counts give a simple proof of equality in the case $n=3$. In Section 4 we consider the case $n=4$, and detail the techniques used to find the equations in Theorem 1. Section 5 concludes with some remarks and remaining questions.

We begin by setting up some notation which will be used in the remainder of the article. We reserve $n$ to denote the size of a square matrix, and all matrices are considered to have real entries. For $n\in{\mathbb{N}}$, let $\operatorname{O}(n)$ (resp. $\operatorname{SO}(n)$) denote the group of $n\times n$ orthogonal (resp. special orthogonal) matrices, which is a real algebraic variety.

Next, the set of doubly stochastic $n\times n$ matrices is defined by the linear conditions $\sum_{i}a_{ij}=\sum_{j}a_{ij}=1$, $a_{ij}\geq 0$, for all $i,j=1,\ldots,n$. This can be interpreted as saying that an $n\times n$ doubly stochastic matrix is uniquely determined by any $(n-1)\times(n-1)$ submatrix obtained by deleting a row and column. Moreover, recall that the set of $n\times n$ doubly stochastic matrices equals the convex hull of all $n\times n$ permutation matrices, which is the so-called Birkhoff polytope ${\mathcal{B}}_{n}$, of dimension $(n-1)^{2}$.

We now introduce the varieties in question. Identifying the space of $n\times n$ real matrices with $n^{2}$-dimensional affine space ${\mathbb{A}}^{n^{2}}_{{\mathbb{R}}}$, we have the coordinate-wise squaring map

Restriction to $\operatorname{O}(n)$ gives a map $\operatorname{O}(n)\to{\mathbb{A}}^{n^{2}}$, whose image lies in ${\mathcal{B}}_{n}$ (note that the image of $\phi$ is contained in the non-negative orthant ${\mathbb{A}}^{n^{2}}_{\geq 0}$). Next, the coordinate projection $\pi:{\mathbb{A}}^{n^{2}}\to{\mathbb{A}}^{(n-1)^{2}}$, given by projecting onto the upper-left $(n-1)\times(n-1)$ submatrix, is injective on ${\mathcal{B}}_{n}$, so we may compose with $\pi$ to obtain $\pi\circ\phi\big{|}_{\operatorname{O}(n)}:\operatorname{O}(n)\to{\mathbb{A}}^{(n-1)^{2}}$. Finally, for both theoretical and practical reasons it is convenient to work in projective space, so taking the projective closure of the image yields the map

We set $Z_{n}:=\overline{\varphi(\operatorname{O}(n))}$, the Zariski-closure of the image of $\varphi$, which is a projective variety in ${\mathbb{P}}^{(n-1)^{2}}_{{\mathbb{R}}}$. Concretely, we view $Z_{n}$ as the projective closure of the set of $(n-1)\times(n-1)$ matrices which are the upper-left submatrix of an $n\times n$ orthostochastic matrix. In this way the linear equations which define (the linear span of) ${\mathcal{B}}_{n}$ are already accounted for in $Z_{n}$, and furthermore, the reduction of variables from $n^{2}$ to $(n-1)^{2}+1$ will be valuable for computation.

Introducing coordinates, we have that the map $\varphi$ above corresponds algebraically to the ring map

We next give the dimension and degree of $Z_{n}$. Note that $\phi$ is a finite map, with general fibers of size $2^{n^{2}}$, corresponding to sign choices on each of the $n^{2}$ entries of a potential preimage. This implies that $\dim Z_{n}=\dim\operatorname{O}(n)$, which we now recall: a matrix is orthogonal iff for all $i,j$, the dot product of the $i^{\text{th}}$ and $j^{\text{th}}$ rows equals $\delta_{ij}$, so

In fact, these ${n\choose 2}+n={n+1\choose 2}$ quadrics form a regular sequence, so $\dim\operatorname{O}(n)=n^{2}-\operatorname{codim}I_{\operatorname{O}(n)}=n^{2}-{n+1\choose 2}={n\choose 2}$.

The degree of $Z_{n}$ is also known (albeit much more recently): the degree of $\operatorname{SO}(n)$ was computed in [4], and in [9], a general formula for the degree of a coordinate-wise power of a variety is given, and combining these yields:

For reference, we tabulate the values of $\dim Z_{n},\deg Z_{n}$ for small values of $n$:

We now review what is known about the defining equations of $Z_{n}$, for $n\leq 3$. For $n=2$, every doubly stochastic matrix is orthostochastic: indeed, a $2\times 2$ doubly stochastic matrix is of the form $A=\begin{pmatrix}a&1-a\\ 1-a&a\end{pmatrix}$ for some $0\leq a\leq 1$, and writing $a=:\cos^{2}\theta$ gives $1-a=\sin^{2}\theta$, so $A$ is the coordinate-wise square of the rotation matrix $\begin{pmatrix}\cos\theta&-\sin\theta\\ \sin\theta&\cos\theta\end{pmatrix}\in\operatorname{SO}(2)$. Thus no equations are needed to define $Z_{2}={\mathbb{P}}^{1}$.

For $n=3$, the variety $Z_{3}$ is a $3$-dimensional variety in ${\mathbb{P}}^{4}$, hence is a hypersurface, and is defined by a single equation. We now show how to find this equation, first found in [12], following the presentation in Section 3 of [7], which will in fact give a set of “naive” equations for any $n$. Consider a $3\times 3$ doubly stochastic matrix

In order for $A$ to be orthostochastic, there must exist sign choices $\epsilon_{1},\epsilon_{2}\in\{\pm 1\}$ such that with these sign choices, the entrywise square roots of the first two columns are orthogonal:

By rearranging the above equation and squaring, we can eliminate one sign choice and all but one square root:

and another rearrangement and squaring produces a polynomial relation without $\epsilon_{i}$:

and finally, using the fact that $A$ is doubly stochastic, we may express $y_{(3,1)}$ (resp. $y_{(3,2)}$) in terms of $y_{(1,1)},y_{(2,1)}$ (resp. $y_{(1,2)},y_{(2,2)}$) and homogenize with respect to $s$ to obtain

This defines a degree $4$ hypersurface in ${\mathbb{P}}^{4}$ which contains $Z_{3}$, and therefore equals $Z_{3}$ by dimension and degree considerations.

In general, one can perform the same procedure for any pair of columns or rows of an $n\times n$ doubly stochastic matrix, to obtain a set of equations which any $n\times n$ orthostochastic matrix must satisfy:

It follows from the discussion above that every $n\times n$ orthostochastic matrix must satisfy the naive equations $C_{i,j},R_{i,j}$. More precisely, a doubly stochastic matrix $A$ satisfies $C_{i,j}=0$ if and only if there exist choices of signs $\epsilon_{1},\ldots,\epsilon_{n-1}$ such that with these choices of signs, the entrywise square roots of the $i^{\text{th}}$ and $j^{\text{th}}$ columns of $A$ become orthogonal. A natural question is:

Stated another way: satisfying a single equation $C_{i,j}$ only guarantees sign choices which work for a single pair of columns. 5 asks whether existence of local sign choices for each pair of columns (and rows), implies existence of a single global sign choice which simultaneously makes all pairs of columns orthogonal.

As an example of what could go wrong, it is conceivable that (after fixing sign choices on the first two columns as required by $C_{1,2}$) the sign choices imposed on the third column by $C_{1,3}$ could differ from those imposed by $C_{2,3}$. This turns out not to happen in the case of $n=3$, as $Z_{3}$ is defined by the vanishing of (the homogenization of) $C_{1,2}$. However, as noted in [7], this good behavior is indeed special only to small values of $n$, as evidenced by the following proposition (which justifies our use of the term naive):

We now arrive at the main focus of this paper, the case $n=4$. In this section, most of the methods and arguments we use will be primarily computational in nature, so we first remark on the techniques used.

Computations were performed in Macaulay2 [10] and Bertini [2], and involve a mix of symbolic and numerical methods. All symbolic computations were done in exact arithmetic, using Gröbner bases over the rational numbers – this is possible since all the varieties in question are defined over ${\mathbb{Q}}$. For the purposes of this article, we regard results obtained by symbolic computation to be on par with those obtained by theoretical proof.

On the other hand, numerical computations were done in floating point, to 53 bits of precision (i.e. standard double-precision floating-point). Results obtained by numerical methods are regarded as correct “with probability 1”: the fact that these programs give reproducible results agreeing with theory in thousands of cases allows us to say with overwhelming confidence that the results obtained in this way are correct. Throughout the section we indicate when a computation was used, as well as the type (symbolic or numerical). We encourage the interested reader to confirm the results of the computations themselves, included in the appendix/auxiliary files.

Returning to the variety at hand: $Z_{4}$ is a $6$-dimensional variety in ${\mathbb{P}}^{9}$ of degree $40$ (we remark that by coincidence, this is the unique value for $n\geq 2$ for which $\deg Z_{n}=\deg\operatorname{SO}(n)$). As before, we have the naive equations $C_{i,j},R_{i,j}$ for $1\leq i<j\leq 4$, which constitute $12$ octics whose vanishing locus (after taking projective closure) contains $Z_{4}$.

We know that $Z_{4}\subseteq Y$, and in fact the two agree generically:

Proposition 9 shows that $Y$ has the same dimension and degree as $Z_{4}$. In particular, $Z_{4}$ (being irreducible) is the unique top-dimensional component of $Y$, and so $Y\neq Z_{4}$ if and only if $Y$ contains additional components of lower dimension. Thus, if $I(Y)=(C_{i,j},R_{i,j})$ denotes the ideal generated by the 12 naive octics, then showing that $I(Y)$ has only one minimal prime (or equivalently, that the radical of $I(Y)$ is prime) would imply $Y=Z_{4}$.

However, the basic issue is that $I(Y)$ is too large to handle directly: half of the octics have 967 terms, and the other half have 6760 terms. Neither symbolic nor numerical methods (via numerical irreducible decomposition) were able to determine the number of minimal primes of $I(Y)$. This motivates our search for other, “smaller” equations of $Z_{4}$.

How does one find equations for a parametrized variety? This procedure is known as implicitization, and amounts to computing the kernel of a ring map. Classically, this is accomplished with Gröbner bases, but in the $n=4$ case, this is infeasible (as of yet). Thus, we try a numerical approach, using interpolation (cf. Theorem 3 in [6]). The method is as follows: since $Z_{4}$ is a projective variety, its defining ideal $I(Z_{4})$ has homogeous generators, so we may search degree by degree. Fixing a degree $d$, we may find a basis of $I(Z_{4})_{d}$ (the degree $d$ part of $I(Z_{4})$) by sampling a large number of points on $Z_{4}$, evaluating all monomials of degree $d$ in ${\mathbb{P}}^{9}$ at these points, and then computing the numerical kernel of the resulting matrix (with rows indexed by points, and columns indexed by monomials). The kernel vectors are then coefficients of degree $d$ forms (linearly independent over ${\mathbb{R}}$) which vanish at all of the sampled points, and if the number of points is large, one expects such a form to vanish on all of $Z_{4}$. There are no linear forms in $I(Z_{4})$ (as these are already accounted for by restricting to ${\mathcal{B}}_{4}$), and similarly we find no forms of degrees $2,3$, and $4$. However, in degree $5$ we find:

Although still too large to comfortably reproduce here, the quintics in $J$ are significantly more manageable than the octics $C_{i,j},R_{i,j}$: 4 of the quintics have 284 terms, and the most complicated involves 454 terms. In fact a Gröbner basis of $J$ can be computed (taking around 8 hours), which yields $\dim X=6$, $\deg X=40$ (the same as $Z_{4}$). Since we also know $Z_{4}\subseteq X$ as well, we may again test whether equality holds by computing the number of minimal primes of $J$, i.e. irreducible components of $X$. Although symbolic computation of minimal primes of $J$ is still too slow to be practical, the key difference from the previous case with $Y$ is that $X$ is small enough for a numerical irreducible decomposition to be feasible. The result, though initially negative, is immediately promising and leads to a resolution of the problem.

The data of a numerical irreducible decomposition of a variety (given by defining equations) consists of a collection of witness sets, each representing an irreducible component of the variety, along with linear equations defining a complementary-dimensional slice of each witness set, and the finitely many intersection points of each witness set with its linear slice, thus encoding the dimension and degree of each irreducible component. Since all the lower-dimensional components of $X$ obtained in the numerical irreducible decomposition were themselves linear spaces (being degree 1), one should expect that equations for them can be found and moreover interpreted as certain classes of $4\times 4$ matrices. This indeed turns out to be the case:

1. (1)

   The 4-dimensional components correspond to the 16 ways of specifying that one entry of a $4\times 4$ doubly stochastic matrix be equal to $1$: this is a codimension $5$ constraint in ${\mathcal{B}}_{4}$, as it forces the remaining entries in that row and column to be $0$. For instance, one such component is defined by the ideal $(y_{(2,1)},y_{(2,2)},y_{(1,3)},y_{(3,3)},y_{(2,3)}-s)$, which (after dehomogenizing $s=1$) requires that the $(2,3)$-entry is $1$. Up to row and column permutations, such matrices are a direct sum of a $3\times 3$ matrix with the $1\times 1$ identity.

2. (2)

   The 5-dimensional components are slightly more complicated: of these, 6 are contained in the hyperplane at infinity $s=0$, and thus are irrelevant for the orthostochastic variety. The 9 relevant components are as follows: fix a pair of indices $(a,a^{\prime})$ with $a,a^{\prime}\in\{2,3,4\}$, and write $\{2,3,4\}\setminus\{a\}=:\{b,c\},\{2,3,4\}\setminus\{a^{\prime}\}=:\{b^{\prime},c^{\prime}\}$. Then the ideal $(y_{(1,1)}-y_{(a,a^{\prime})},y_{(a,1)}-y_{(1,a^{\prime})},y_{(b,1)}-y_{(c,a^{\prime})},y_{(1,b^{\prime})}-y_{(a,c^{\prime})})$ defines a 5-dimensional component of $X$. Note that these components can be described by partitioning a $4\times 4$ doubly stochastic matrix into 4 disjoint $2\times 2$ submatrices, and requiring each submatrix to be circulant (i.e. has equal cross terms).

These results were obtained by following the same procedure as in the proof of Proposition 10 (sampling points, obtaining approximate linear equations, and then extracting exact linear equations via LLL), and the resulting equations are symbolically verified to define subvarieties of $X$. In the case of the 5-dimensional components, an additional, carefully chosen, change-of-basis was necessary to obtain binomial linear generators, and with it the ensuing description as matrices. The explicit linear ideals are recorded (as o22 and o24) in the appendix below.

Each of the relevant lower-dimensional components of $X$ can thus be defined over ${\mathbb{Z}}$ and interpreted as a particular class of matrices. With these descriptions via integer equations at hand, it is also straightforward to confirm that each of these classes of matrices defines an actual component of $X$, via checking dimensions of tangent spaces. Finally, we may use this to obtain Theorem 1:

A few remarks on the methodology used in Section 4 are in order. First, the only part of the proof of Theorem which relies on numerical computation (i.e. for which there is no symbolic verification) is the correctness of the numerical irreducible decomposition Proposition 12. Specifically, what is required is that there are no other irreducible components of $X$ other than those listed in Proposition 12. For readers concerned about the random choices involved in path tracking and homotopy continuation, we mention that the numerical irreducible decomposition of $X$ has been calculated multiple times, each time giving the same consistent result. Moreover, the numerical irreducible decomposition of the affine variety $X_{s=1}$ (the dehomogenization of $X$) has also been computed, and is consistent with Proposition 12 (e.g. among the $5$-dimensional components, only the $9$ relevant components not at infinity are present).

One may also ask about the utility of set-theoretic equations for $Z_{4}$. After all, given a particular $4\times 4$ matrix, the brute-force check to determine if it is orthostochastic is still feasible (namely checking all possible sign patterns in a fiber of $\phi$. Since $({\mathbb{Z}}/2{\mathbb{Z}})^{n}\oplus({\mathbb{Z}}/2{\mathbb{Z}})^{n}$ acts on fibers of $\phi$ by row and column sign changes, one may assume the first row and column are all nonnegative, so there are only $2^{(4-1)^{2}}=512$ sign patterns to check). In response to this, we note that in addition to their intrinsic interest, knowledge of set-theoretic equations is extremely useful for many other purposes. For instance, given a variety which meets the set of $4\times 4$ orthostochastic matrices, one can now describe their common intersection locus (which was previously not possible to do). We plan to use this in a forthcoming upgrade to [5], to compute monic symmetric determinantal representations of hyperbolic plane quartic curves.

As for the equations themselves, it is natural to ask about minimality. For degree reasons, any set of minimal generators of $I(Z_{4})$ must include the $6$ quintics $J$ in Definition 11. It has been checked that the spaces of sextics and 7-forms in $I(Z_{4})$ are $60$- and $330$-dimensional respectively, and moreover that the quintics in $J$ have no quadratic (or linear) syzygies. It follows that there are no forms of degree $6$ or $7$ in $I(Z_{4})$ other than those already in $J$. For degree $8$, it is not difficult to see that $2$ of the naive octics (along with $J$) would not be sufficient to cut out $Z_{4}$, so the generating set in Theorem 1 is both inclusion-wise and degree-wise minimal.

Finally, we list some questions that remain unanswered with this work. First, it would be interesting to have a combinatorial description of the $6$ quintics $J$, as well as a theoretical proof for why they vanish on $4\times 4$ orthostochastic matrices. Next, for $n=4,5$ it remains open whether the naive equations cut out the orthostochastic variety inside the Birkhoff polytope set-theoretically. Lastly, we know that for $n\geq 6$ additional equations are needed beyond the naive ones. However, finding where these come from is related to the first question listed here, and even a computational proof of sufficiency may only be possible with some advances in computing technology.

The first author would like to thank Anton Leykin for helpful discussions, especially related to verifying components via tangent spaces. The second author would like to thank Paul Breiding, Mateusz Michalek and Bernd Sturmfels for their encouragement towards work on this project.