Computing symmetric determinantal representations

TL;DR

This paper presents a Macaulay2 package that efficiently computes definite symmetric determinantal representations for real polynomials, focusing on quadrics and low-degree plane curves, using linear algebra and numerical algebraic geometry.

Contribution

Introduction of a new Macaulay2 package that computes symmetric determinantal representations without genericity assumptions, emphasizing speed and robustness.

Findings

Successfully computes representations for quadrics and low-degree curves.

Employs linear algebra and numerical methods for robustness.

No assumptions on polynomial genericity.

Abstract

We introduce the DeterminantalRepresentations package for Macaulay2, which computes definite symmetric determinantal representations of real polynomials. We focus on quadrics and plane curves of low degree (i.e. cubics and quartics). Our algorithms are geared towards speed and robustness, employing linear algebra and numerical algebraic geometry, without genericity assumptions on the polynomials.

Equations

$$\hat{D}(\underset{k_1}{\underbrace{A^{(1)},\dots ,A^{(1)}}},\underset{k_2}{\underbrace{A^{(2)},\dots ,A^{(2)}}},\dots ,\underset{k_n}{\underbrace{A^{(n)},\dots ,A^{(n)}}}):=\sum _{\alpha \in S[m]}\sum _{\sigma \in \tilde{S}}\det \left[\begin{array}{ccc}{a}_{{\alpha }_1{\alpha }_1}^{(\sigma (1))}& \dots & a_{{\alpha }_1{\alpha }_k}^{(\sigma (1))}\\ ⋮& & \\ a_{{\alpha }_k{\alpha }_1}^{(\sigma (k))}& \dots & a_{{\alpha }_k{\alpha }_k}^{(\sigma (k))}\end{array}\right]$$

Full text

1. Introduction

The problem of representing a polynomial as the determinant of a linear matrix pencil is classical, cf. [1, 3, 7, 8, 9]. A polynomial $f\in{\mathbb{R}}[x_{1},\ldots,x_{n}]$ of degree $d$ (not necessarily homogeneous) is called determinantal if $f$ is the determinant of a matrix with linear entries – i.e. there exist matrices $A_{0},\ldots,A_{n}\in{\mathbb{R}}^{d\times d}$ such that $f(x_{1},\ldots,x_{n})=\det(A_{0}+x_{1}A_{1}+\ldots+x_{n}A_{n})$. The matrix $A_{0}+x_{1}A_{1}+\ldots+x_{n}A_{n}$ is said to give a determinantal representation of $f$ of size $d$. If the matrices $A_{i}$ can be chosen to be all symmetric (resp. Hermitian), then the determinantal representation is called symmetric (resp. Hermitian). The determinantal representation is called definite if $A_{0}$ is definite, and monic if $A_{0}=I_{d}$ is the identity matrix.

Computing definite symmetric (resp. Hermitian) determinantal representations of a polynomial is known as the determinantal representation problem in convex algebraic geometry [2]. It has generated interest to the optimization community due to its connection with the problem of determining (definite) linear matrix inequality (LMI) representable sets ([9, 12]). The problem of characterizing the LMI-representable subsets of ${\mathbb{R}}^{n}$ (i.e. spectrahedra) can be answered by characterizing determinantal polynomials, which leads to the generalized Lax conjecture [10].

Throughout we focus mainly on homogeneous polynomials, typically in $3$ variables, corresponding to projective plane curves (though internally via de-homogenization, it suffices to compute determinantal representations for bivariate polynomials). By a celebrated theorem of Helton-Vinnikov ([9, 10]), all hyperbolic polynomials in $3$ variables admit symmetric determinantal representations. When $n\geq 4$, a general homogeneous polynomial of degree $d$ in $n$ variables does not admit any determinantal representation of size $d$ (except for $(n,d)=(4,3)$). We abbreviate monic symmetric (resp. Hermitian) determinantal representation to MSDR (resp. MHDR).

2. Quadratic determinantal polynomials

For a quadratic polynomial $f({\mathbf{x}})={\mathbf{x}}^{T}A{\mathbf{x}}+b^{T}{\mathbf{x}}+1\in{\mathbb{R}}[{\mathbf{x}}]$ in $n$ (= any number of) variables, one can give a necessary and sufficient condition for a MSDR (resp. MHDR) to exist, via an associated matrix $W:=A-\frac{1}{4}bb^{T}$; namely one (or both) of the following conditions must hold:

1. (1)

$W$ is negative semidefinite with $\operatorname{rank}W\leq 3$ 2. (2)

$A$ is negative semidefinite.

One can explicitly find a MSDR (resp. MHDR) of size $2$ in case (1), and one of size at most $n+1$ in case (2), cf. ([6], Theorem 3.4) for an algorithm. This is implemented as follows:

Macaulay2, version 1.13 i1 : needsPackage "DeterminantalRepresentations" i2 : R = CC[x_1..x_4] i3 : f = -25x_1^2 + 254x_1x_2 + 243x_2^2 + 234x_1x_3 + 494x_2x_3 + 247x_3^2 + 198x_1x_4 + 378x_2x_4 + 378x_3x_4 + 143x_4^2 + 18x_1 + 32x_2 + 32x_3 + 24x_4 + 1 i4 : quadraticDetRep f o4 = | 8.06514x_1+19.0421x_2+17.2128x_3+12.4043x_4+1 10.2531x_1+1.93541x_2+2.74393x_3+... | | 10.2531x_1+1.93541x_2+2.74393x_3+.914643x_4 9.93486x_1+12.9579x_2+14.7872x_3+... |

3. Generalized mixed discriminant

The relations between the coefficients of a determinantal polynomial and the entries of the coefficient matrices are captured by the generalized mixed discriminant:

Definition 1**.**

Let $\{A^{(1)},\ldots,A^{(n)}\}=\{(a_{ij}^{(1)}),\ldots,(a_{ij}^{(n)})\}$ be a set of $n\times n$ matrices of size $n$. The generalized mixed discriminant of a tuple of matrices $(\underbrace{A^{(1)},\dots,A^{(1)}}_{k_{1}},\underbrace{A^{(2)},\dots,A^{(2)}}_{k_{2}},\dots,\underbrace{A^{(n)},\dots,A^{(n)}}_{k_{n}})$ is defined as

[TABLE]

where $m$ is the number of distinct matrices, $S[m]$ is the set of order-preserving $m$-cycles in $S_{n}$ (i.e. $\alpha=(\alpha_{1},\dots,\alpha_{m})\in S[m]\Rightarrow\alpha_{1}<\alpha_{2}<\dots<\alpha_{m}$), and $\widetilde{S}$ is the set of all distinct permutations of $\{\underbrace{1,\dots,1}_{k_{1}},\dots,\underbrace{n,\dots,n}_{k_{n}}\}$.

Theorem 3.1**.**

*(cf. [5], Theorem 2.5) If $f=\det(I_{d}+\sum x_{i}A_{i})$ is determinantal of degree $d$, then the coefficients of $f$ are given by the generalized mixed discriminants of the matrices $A_{i}$: namely the coefficient of $x_{1}^{k_{1}}\dots x_{n}^{k_{n}}$ in $f$ is equal to $\widehat{D}(\underbrace{A_{1},\dots,A_{1}}_{k_{1}},\underbrace{A_{2},\dots,A_{2}}_{k_{2}},\dots,\underbrace{A_{n},\dots,A_{n}}_{k_{n}}).$ *

Generalized mixed discriminants can be computed with the function generalizedMixedDiscriminant:

i5 : n = 3; S = QQ[a_(1,1)..a_(n,n),b_(1,1)..b_(n,n),c_(1,1)..c_(n,n)]; R = S[x_1..x_n] i8 : A = sub(genericMatrix(S, n, n), R) i9 : B = sub(genericMatrix(S, b_(1,1), n, n), R) i10 : C = sub(genericMatrix(S, c_(1,1), n, n), R) i11 : G = generalizedMixedDiscriminant({A, B, C}) o11 = a b c - a b c - a b c + a b c - a b c + ... 3,3 2,2 1,1 3,2 2,3 1,1 2,3 3,2 1,1 2,2 3,3 1,1 3,3 2,1 1,2 i12 : P = det(id_(R^n) + x_1A + x_2B + x_3C); i13 : G == (last coefficients(P, Monomials => {x_1x_2*x_3}))_(0,0) o13 = true

4. Higher degree determinantal polynomials

We now consider polynomials of degree $>2$. Let $f\in{\mathbb{R}}[x_{0},x_{1},x_{2}]_{d}$ be a homogeneous polynomial of degree $d$ in $3$ variables. The homogenization of a MSDR for f\Big{|}_{x_{0}=1} gives a MSDR for $f$, so we reduce to the case of $2$ variables. Given a (not necessarily homogeneous) polynomial $f\in{\mathbb{R}}[x_{1},x_{2}]$ of total degree $d$, we seek to compute $A_{1},A_{2}\in\operatorname{Sym}_{d\times d}({\mathbb{R}})$ such that $f=\det(I_{d}+x_{1}A_{1}+x_{2}A_{2})$. It is easy to obtain the eigenvalues of the unknown matrices $A_{i}$: for instance, the eigenvalues of $A_{1}$ are the negative reciprocals of the roots of the univariate polynomial f\Big{|}_{x_{2}=0}=\det(I_{d}+x_{1}A_{1}), and similarly for $A_{2}$ (note that this polynomial has nonzero roots).

Now, if the $A_{i}$ are symmetric, then by the spectral theorem there are orthogonal matrices $V_{i}$ such that $V_{i}^{T}A_{i}V_{i}=D_{i}$ is diagonal, with entries equal to the known eigenvalues of $A_{i}$, for $i=1,2$. Setting $V:=V_{2}^{T}V_{1}$ (which is orthogonal), one has $f=\det(I_{d}+x_{1}D_{1}+x_{2}V^{T}D_{2}V)=\det(I_{d}+x_{1}VD_{1}V^{T}+x_{2}D_{2})$.

With this, one can obtain the diagonal entries of $V^{T}D_{2}V$: it follows from Theorem 3.1 that these can be obtained by solving a linear system involving only $D_{1}$ (i.e. eigenvalues of $A_{1}$ – note the subscript) as well as coefficients of monomials in $f$ which are linear in $x_{2}$. We note that the linear system giving rise to diagonal entries of $V^{T}D_{2}V$ has unique solutions iff $A_{1}$ has distinct eigenvalues – if the solutions are nonunique, then (as we will see) we must choose a solution which is majorized by (the diagonal entries of) $D_{2}$. By a symmetrical argument, we may henceforth assume that diagonal entries of $V^{T}D_{2}V$ and $VD_{1}V^{T}$ are known.

It thus suffices to compute the symmetric matrix $V^{T}D_{2}V$: we propose two methods to do so. The first, brute-force, method is to simply find the $d\choose 2$ off-diagonal entries of $V^{T}D_{2}V$ by solving a square polynomial system arising from Theorem 3.1. Although symbolic methods are slow, a numerical method for solving the system is implemented in the method trivariateDetRep, with the (default) option Strategy => DirectSystem:

i14 : R = RR[x,y,z] i15 : f = det(xid_(R^4) + ydiagonalMatrix {4,3,2,1_R} + z*randomIntegerSymmetric(4,R)) i16 : sols = trivariateDetRep f; Solving 6 x 6 polynomial system ... i17 : #sols o17 = 64 i18 : all(sols, M -> clean(1e-7, f - det M) == 0) o18 = true

This method works well up to degree $4$, but for higher degrees, even numerical methods take too long to finish.

An alternative, more theoretical, method for cubics is to note that since $D_{2}$ is known, finding $V^{T}D_{2}V$ is equivalent to finding $V$, and in this case (of degree $3$), it turns out that the Hadamard square of $V$ can be determined essentially by linear algebra (cf. [4], Theorem 2.14).

Definition 2**.**

If $A,B$ are matrices of the same size, we denote their Hadamard product by $A\odot B$, i.e. $(A\odot B)_{ij}=A_{ij}B_{ij}$. We say that a square matrix $A$ is orthostochastic if it is the Hadamard square of an orthogonal matrix, i.e. $A=V\odot V$ for some orthogonal matrix $V$. If $v,w\in{\mathbb{R}}^{n}$, we say that $v$ is majorized by $w$ if $\sum_{j=1}^{n}v_{j}=\sum_{j=1}^{n}w_{j}$ and $\sum_{j=1}^{i}(\widetilde{v})_{j}\leq\sum_{j=1}^{i}(\widetilde{w})_{j}$ for all $i=1,\ldots,n$, where $\widetilde{v},\widetilde{w}$ are the decreasing rearrangements of $v,w$.

To elaborate: in the cubic case, given the majorization conditions mentioned above, one can set up a zero-dimensional polynomial system to find the unknown entries of $V\odot V$, which consists of one cubic equation with all other equations linear, and only involves the (known) diagonal entries of $D_{1},D_{2},V^{T}D_{2}V,VD_{1}V^{T}$.

Thus to recover $V$, it suffices to determine all orthogonal matrices with Hadamard square equal to a given orthostochastic matrix. Given a $n\times n$ orthostochastic matrix $A$, there are $2^{n^{2}}$ possible matrices whose Hadamard square is $A$ (not all of which will be orthogonal in general though). Let $G\cong({\mathbb{Z}}/2{\mathbb{Z}})^{n}$ be the group of diagonal matrices with diagonal entries equal to $\pm 1$. Then $G\times G$ acts on the set of orthogonal matrices whose Hadamard square is $A$, via $(g_{1},g_{2})\cdot O:=g_{1}Og_{2}$. The method orthogonalFromOrthostochastic computes all such orthogonal matrices, modulo the action of $G\times G$. We note the following:

Proposition 4.1**.**

For a general orthostochastic matrix $A$, there is exactly one $G\times G$-orbit of orthogonal matrices with Hadamard square equal to $A$.

Proof.

Since the action of $G\times G$ amounts to performing sign changes in each row and column, each $G\times G$ orbit contains a unique element with nonnegative entries in the first row and first column. If there were distinct $G\times G$ orbits for $A$, then there would be orthogonal matrices $U,V$ with identical first column $u_{(i,1)}=v_{(i,1)}=\sqrt{a_{(i,1)}}$, but $u_{j}\neq\pm v_{j}$ for some column $j>1$ (since $u_{(1,j)}=v_{(1,j)}=\sqrt{a_{(1,j)}}\geq 0$), which would impose a Zariski-closed condition on the entries of $U$ (namely $u_{1}^{T}(u_{j}-v_{j})=0$), hence does not hold for general $A$. ∎

We illustrate this with some examples: note that when using floating point inputs, it may be necessary to specify the value of the option Tolerance (default value $10^{-5}$) in order to obtain useful results.

i19 : (A1, A2) = (randomIntegerSymmetric(3, R), randomIntegerSymmetric(3, R)) o19 = (| 12 3 18 |, | 18 16 9 |) | 3 12 14 | | 16 4 11 | | 18 14 6 | | 9 11 16 | i20 : f = det(xid_(R^3) + yA1 + zA2) 3 2 2 3 2 2 2 2 3 o20 = x + 30x y - 241xy - 3918y + 38x z + 52xyz + 768y z - 34xz + 3282yz - 2278z i21 : reps = trivariateDetRep f o21 = {| x+33.7014y+36.8578z .607983z 3.51039z |, | .607983z x+9.08918y-6.37807z -3.92965z | | 3.51039z -3.92965z x-12.7906y+7.52028z | ------------------------------------------------------------------ | x+33.7014y+36.8578z 1.11917z 3.2642z |} | 1.11917z x+9.08918y-6.37807z 4.02828z | | 3.2642z 4.02828z x-12.7906y+7.52028z | i22 : all(reps, M -> clean(1e-9, f - det M) == 0) o22 = true i23 : g = x^3+7x^2y+16xy^2+12y^3+3x^2z+22xyz+32y^2z-45xz^2-65yz^2-175*z^3 o24 : reps = trivariateDetRep g o24 = {| x+3y-7z 0 0 |} | 0 x+2y+5z 0 | | 0 0 x+2y+5z |

As seen above, the output of trivariateDetRep is a list of matrices $M$, whose entries are linear forms, such that the input polynomial $f$ is equal to $\det M$. As the second example shows, no genericity assumptions are made on the polynomial $f$ – that is, the eigenvalues of the coefficient matrices $A_{1}$ and $A_{2}$ need not be distinct (as opposed to the treatment in e.g. [11]). This in turn reveals information which is typically hidden to numerical methods (such as a numericalIrreducibleDecomposition), e.g. that one of the component lines of the plane curve defined by $g$ above has multiplicity $2$.

5. Additional methods

In the course of creating this package, various functions for working with matrices were needed which (to the best of our knowledge) were not available in Macaulay2. Thus a number of helper functions are included in this package, which may be of general interest to users beyond the scope of computing determinantal representations. These include: hadamard (for computing Hadamard products), cholesky (for computing the Cholesky decomposition of a PSD matrix), companionMatrix (which returns a matrix whose characteristic polynomial is any given univariate monic polynomial), isOrthogonal, isDoublyStochastic (for checking properties of a given matrix), randomIntegerSymmetric, randomUnipotent, randomOrthogonal, randomPSD (for generating various types of random matrices), and liftRealMatrix, roundMatrix (for lifting matrices from ${\mathbb{C}}$ to ${\mathbb{R}}$ and ${\mathbb{Q}}$).

Acknowledgements. Both authors gratefully acknowledge the support of ICERM, and the Fall 2018 Nonlinear Algebra program in particular, where this project began. We thank Ritvik Ramkumar for helpful comments and software testing. The second author would like to thank her Ph.D. thesis supervisor Harish K. Pillai for his helpful suggestions in this topic and Deepak Patil for implementing the notion of generalized mixed discriminant of matrices in Matlab during her Ph.D.