This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

Unlabeled Sensing Using Rank-One Moment Matrix Completion

Hao Liang KLMM, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China University of Chinese Academy of Sciences [email protected] Jingyu Lu KLMM, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China University of Chinese Academy of Sciences [email protected] Manolis C. Tsakiris KLMM, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China University of Chinese Academy of Sciences [email protected]  and  Lihong Zhi KLMM, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China University of Chinese Academy of Sciences [email protected]
Abstract.

We study the unlabeled sensing problem that aims to solve a linear system of equations Ax=π(y)Ax=\pi(y) for an unknown permutation π\pi. For a generic matrix AA and a generic vector yy, we construct a system of polynomial equations whose unique solution satisfies Aξ=π(y)A\xi^{*}=\pi(y). In particular, ξ\xi^{*} can be recovered by solving the rank-one moment matrix completion problem. We propose symbolic and numeric algorithms to compute the unique solution. Some numerical experiments are conducted to show the efficiency and robustness of the proposed algorithms.

Key words and phrases:
unlabeled sensing, birational morphism, matrix completion, determinantal variety, semidefinite relaxation

1. Introduction

Given a full-rank matrix Am×nA^{*}\in\mathbb{R}^{m\times n} and a vector ysy^{*}\in\mathbb{R}^{s} satisfying ms>n>0m\geqslant s>n>0, the unlabeled sensing problem [63, 64] asks that for an unknown vector ξn\xi^{*}\in\mathbb{R}^{n}, if one only knows the vector ysy^{*}\in\mathbb{R}^{s} consisting of ss shuffled entries of AξA^{*}\xi^{*}, whether the vector ξ\xi^{*} is unique and how to recover it efficiently. This problem emerges from various fields of natural science and engineering, such as biology [24, 55, 1, 38], neuroscience [45], computer vision [17, 41, 61, 26, 35] and communication networks [44, 27, 57].

Theorem 1 in [63] asserts that the solution of the unlabeled sensing problem is indeed unique if s2ns\geqslant 2n and AA^{*} is generic. For the case m=sm=s, Song, Choi and Shi [57] proposed a new method for recovering the vector ξ\xi^{*} as follows: let

qi(x)=pi(Ax)pi(y),q_{i}(x)=p_{i}(A^{*}x)-p_{i}(y^{*}),

where x=[x1,,xn]x=[x_{1},\dots,x_{n}] and

pk(y)=i=1myik[y]p_{k}(y)=\sum_{i=1}^{m}{y_{i}^{k}}\in\mathbb{R}[y]

is the kk-th power sum of the variables y=[y1,,ym]y=[y_{1},\dots,y_{m}]. As pkp_{k} is a symmetric polynomial in yy, its value is independent of the order of the variables y1,,ymy_{1},\ldots,y_{m}. Suppose π\pi is an element of the permutation group Σm\Sigma_{m}, if ξ\xi^{*} is a solution of

Ax=π(y),A^{*}x=\pi(y^{*}), (1)

then it is a root of the polynomials qi(x)q_{i}(x), i.e.,

qi(ξ)=pi(Aξ)pi(y)=pi(π(y))pi(y)=0q_{i}(\xi^{*})=p_{i}(A^{*}\xi^{*})-p_{i}(y^{*})=p_{i}(\pi(y^{*}))-p_{i}(y^{*})=0

for i=1,,mi=1,\ldots,m.

By Theorem 1 in [57] and by [63, 64], with Am×nA^{*}\in\mathbb{R}^{m\times n} a given matrix with i.i.d random entries drawn from an arbitrary continuous probability distribution over \mathbb{R}, if m2nm\geqslant 2n, then with probability 11, ξ\xi^{*} is the unique solution of the polynomial system

Qm={q1(x)=0,,qm(x)=0}.Q_{m}=\{q_{1}(x)=0,\ldots,~{}q_{m}(x)=0\}. (2)

Numerical experiments showed that solving the first n+1n+1 equations is sufficient for recovering the solution ξ\xi^{*}; i.e. ξ\xi^{*} is the unique solution of Qn+1Q_{n+1}. Hence, [57] pointed out the following open problem (see also Conjecture 6 in [42]):

Open Problem.

For a generic matrix Am×nA^{*}\in\mathbb{R}^{m\times n}, mn+1m\geqslant n+1, and a permutation yy^{*} of a generic vector in the column space of AA^{*}, the (n+1)(n+1)-by-nn polynomial system

Qn+1={q1(x)=0,,qn(x)=0,qn+1(x)=0}Q_{n+1}=\{q_{1}(x)=0,\ldots,~{}q_{n}(x)=0,~{}q_{n+1}(x)=0\} (3)

has a unique solution ξ\xi^{*}, which satisfies Aξ=π(y)A^{*}\xi^{*}=\pi(y^{*}).

In [62], it is shown that for a generic matrix Am×nA^{*}\in\mathbb{C}^{m\times n} and any vector yy^{*}, the nn-by-nn polynomial system

Qn={q1(x)=0,,qn(x)=0}Q_{n}=\{q_{1}(x)=0,\ldots,~{}q_{n}(x)=0\} (4)

in the variables x=x1,,xnx=x_{1},\dots,x_{n} has at most n!n! solutions. Furthermore, if yy^{*} is a permutation of a generic vector in the column space of AA^{*}, then among the solutions of QnQ_{n} there is only one vector ξ\xi^{*} satisfying Aξ=π(y)A^{*}\xi^{*}=\pi(y^{*}). They can recover ξ\xi^{*} by solving QnQ_{n} via symbolic or homotopy methods and then select ξ\xi^{*} from n!n! solutions via numerical optimization methods.

Example 1.

Given a matrix AA^{*} and a vector yy^{*}

A:=[12430220],y=[51024],A^{*}:=\left[\begin{array}[]{cc}1&2\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 4&3\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&-2\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr-2&0\end{array}\right],~{}~{}y^{*}=\left[\begin{array}[]{c}-5\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr-10\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 2\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 4\end{array}\right],

find a solution ξ\xi^{*} such that

Aξ=π(y)A^{*}\xi^{*}=\pi(y^{*})

for an unknown permutation π\pi of the coordinates of yy^{*}.

From the matrix AA^{*} and the vector yy^{*}, we compute the polynomials

q1(x)\displaystyle q_{1}(x) =\displaystyle= 3x1+3x2+9,\displaystyle 3\,x_{{1}}+3\,x_{{2}}+9,
q2(x)\displaystyle q_{2}(x) =\displaystyle= 21x12+28x1x2+17x22145,\displaystyle 21\,{x_{{1}}}^{2}+28\,x_{{1}}x_{{2}}+17\,{x_{{2}}}^{2}-145,
q3(x)\displaystyle q_{3}(x) =\displaystyle= 57x13+150x12x2+120x1x22+27x23+1053,\displaystyle 57\,{x_{{1}}}^{3}+150\,{x_{{1}}}^{2}x_{{2}}+120\,x_{{1}}{x_{{2}}}^{2}+27\,{x_{{2}}}^{3}+1053,
q4(x)\displaystyle q_{4}(x) =\displaystyle= 16x13x2+44x12x22+24x1x23400.\displaystyle 16\,{x_{{1}}}^{3}x_{{2}}+44\,{x_{{1}}}^{2}{x_{{2}}}^{2}+24\,x_{{1}}{x_{{2}}}^{3}-400.

The polynomial system {q1(x)=0,q2(x)=0,q3(x)=0,q4(x)=0}\{q_{1}(x)=0,q_{2}(x)=0,q_{3}(x)=0,q_{4}(x)=0\} has a unique solution ξ=(1,2).\xi^{*}=(-1,-2). One can check that ξ\xi^{*} is the solution of (1) for the permutation π=[1,2,4,3]\pi=[1,2,4,3].

We can check that ξ\xi^{*} is the unique solution of the polynomial system {q1(x)=0,q2(x)=0,q3(x)=0}\{q_{1}(x)=0,q_{2}(x)=0,q_{3}(x)=0\}, which gives an example supporting the positive answer to the open problem (Open Problem).

The polynomial system {q1(x)=0,q2(x)=0}\{q_{1}(x)=0,q_{2}(x)=0\} has two solutions η=(45,115),ξ=(1,2).\eta^{*}=\left(-{\frac{4}{5}},-{\frac{11}{5}}\right),~{}~{}\xi^{*}=(-1,-2). One can check that η\eta^{*} is not a solution of (1).

The following theorem shows that one can obtain the unique and desired vector ξ\xi^{*} for a generic matrix AA^{*} and a permutation yy^{*} of a generic vector in the column space of AA^{*} by solving the polynomial system Qn+1Q_{n+1}.

Theorem 1.

The map

f:m×n×n\displaystyle f:\mathbb{C}^{m\times n}\times\mathbb{C}^{n} m×n×n+1\displaystyle\rightarrow\mathbb{C}^{m\times n}\times\mathbb{C}^{n+1} (5)
(A,ξ)\displaystyle(A^{*},\xi^{*}) (A,p1(Aξ),,pn+1(Aξ))\displaystyle\mapsto(A^{*},p_{1}(A^{*}\xi^{*}),\dots,p_{n+1}(A^{*}\xi^{*}))

is a birational morphism to the image f(m×n×n)m×n×n+1f(\mathbb{C}^{m\times n}\times\mathbb{C}^{n})\subseteq\mathbb{C}^{m\times n}\times\mathbb{C}^{n+1}.

Remark 1.

The fact that the morphism (5) is birational implies that there exists a dense open subset U2U_{2} of VV such that ff induces a bijection from U1:=f1(U2)U_{1}:=f^{-1}(U_{2}) to U2U_{2}, see [22, Chapter I, Corollary 4.5]. Therefore, if the sample point (A,y)(A^{*},y^{*}) satisfies the condition (A,p1(y),,pn+1(y))U2(A^{*},p_{1}(y^{*}),\dots,p_{n+1}(y^{*}))\in U_{2}, then there exists a unique vector ξ\xi^{*} such that (A,ξ)U1(A^{*},\xi^{*})\in U_{1} and

(A,p1(Aξ),,pn+1(Aξ))=(A,p1(y),,pn+1(y)).(A^{*},p_{1}(A^{*}\xi^{*}),\dots,p_{n+1}(A^{*}\xi^{*}))=(A^{*},p_{1}(y^{*}),\dots,p_{n+1}(y^{*})).

This answers the open problem above in the affirmative.

As an alternative to existing RANSAC [18], expectation maximization [2], branch & bound [52] or homotopy and Groebner bases [62] approaches for solving the unlabeled sensing problem in its generality (that is without special assumptions such as sparsity, as, e.g., considered in [56]), in this paper we focus on finding the unique solution of the polynomial system Qn+1Q_{n+1} efficiently by reducing it to a rank-one moment matrix completion problem, which can then be solved by many efficient algorithms.

Structure of the paper

In Section 2, we briefly review the main idea in the proof of Theorem 1 in [37] and its relationship with previous known results. In Section 3, we show that solving polynomial system Qn+1Q_{n+1} can be reduced to solving the rank-one moment matrix completion problem, which can be solved efficiently by solving a sequence of semidefinite programming (SDP) problems. In Section 4, we propose some symbolic and numerical algorithms based on Theorem 1 and LABEL:cor. We conducted eight numeric and symbolic experiments to test the efficiency and robustness of the proposed algorithms.

2. Proof of Theorem 1

For convenience, we will abbreviate the unlabeled sensing problem for case m=s>n>0m=s>n>0 by USP and call the data (A,y)(A^{*},y^{*}) in USP as the sample point of USP.

The morphism ff can be written as the composition of the following two morphisms

g:m×n×n\displaystyle g:\mathbb{C}^{m\times n}\times\mathbb{C}^{n} g(m×n×n)¯m×n×m\displaystyle\rightarrow\overline{g(\mathbb{C}^{m\times n}\times\mathbb{C}^{n})}\subseteq\mathbb{C}^{m\times n}\times\mathbb{C}^{m} (6)
(A,ξ)\displaystyle(A^{*},\xi^{*}) (A,p1(Aξ),,pm(Aξ))\displaystyle\mapsto(A^{*},p_{1}(A^{*}\xi^{*}),\dots,p_{m}(A^{*}\xi^{*}))
andγ:g(m×n×n)¯\displaystyle\text{and}~{}\gamma:\overline{g(\mathbb{C}^{m\times n}\times\mathbb{C}^{n})} f(m×n×n)¯m×n×m.\displaystyle\rightarrow\overline{f(\mathbb{C}^{m\times n}\times\mathbb{C}^{n})}\subseteq\mathbb{C}^{m\times n}\times\mathbb{C}^{m}. (7)
(A,p1(η),,pm(η))\displaystyle(A^{*},p_{1}(\eta^{*}),\dots,p_{m}(\eta^{*})) (A,p1(η),,pn+1(η))\displaystyle\mapsto(A^{*},p_{1}(\eta^{*}),\dots,p_{n+1}(\eta^{*})) (8)

Since the composition of two birational morphisms is still a birational morphism, the proof of Theorem 1 can be divided into two parts:

  • Show the morphism gg is birational, and

  • Show the morphism γ\gamma is birational.

When m=n+1m=n+1, the morphism gg equals the morphism ff in Theorem 1. We show first in Section 2.1 that the morphism gg is birational. Hence, it implies Theorem 1 for the case m=n+1m=n+1, i.e., one can obtain the unique and desired vector ξ\xi^{*} for a generic matrix AA^{*} and a permutation yy^{*} of a generic vector in the column space of AA^{*} by solving the polynomial system Qn+1Q_{n+1} for the case m=n+1m=n+1. Moreover, as gg is birational, it implies that for generic sample points (A,y)(A^{*},y^{*}) of the unlabeled sensing problem (USP), the polynomial system QmQ_{m} in (2) has a unique solution x=ξx=\xi^{*}, which is also the unique solution of USP.

Since mm can be much larger than nn in the application, we need to show that ff is birational for a general mn+1m\geqslant n+1. This implies that the smaller polynomial system Qn+1Q_{n+1} in (3) has a unique solution for generic (A,y)(A^{*},y^{*}). The key point of showing that ff is birational is to prove

[(A,x):(A,p1(Ax),,pn(Ax))]=n!,and\displaystyle[\mathbb{C}(A,x):\mathbb{C}(A,p_{1}(Ax),\dots,p_{n}(Ax))]=n!,~{}\text{and}
[(A,p1(Ax),,pn+1(Ax)):(A,p1(Ax),,pn(Ax))]=n!.\displaystyle[\mathbb{C}(A,p_{1}(Ax),\dots,p_{n+1}(Ax)):\mathbb{C}(A,p_{1}(Ax),\dots,p_{n}(Ax))]=n!.

The first degree is the Bézout number of the complete intersection given by the regular sequence p1(Ax),,pn(Ax)p_{1}(Ax),\dots,p_{n}(Ax), and the second one is exactly the degree of the minimal polynomial of pn+1(Ax)p_{n+1}(Ax) over (A,p1(Ax),,pn(Ax))\mathbb{C}(A,p_{1}(Ax),\dots,p_{n}(Ax)). In Section 2.2, we give a sketch proof showing that this polynomial is the resultant of the system {pk(Ax)zk:k=1,,n+1}\{p_{k}(Ax)-z_{k}:k=1,\dots,n+1\}.

2.1. The morphism gg is birational

In this subsection, we prove the following theorem, which implies that one can obtain the unique and desired vector ξ\xi^{*} for a generic matrix AA^{*} and a permutation yy^{*} of a generic vector in the column space of AA^{*} by solving the polynomial system QmQ_{m}.

Theorem 2.

For m>n>0m>n>0, the morphism

g:m×n×n\displaystyle g:\mathbb{C}^{m\times n}\times\mathbb{C}^{n} m×n×m\displaystyle\rightarrow\mathbb{C}^{m\times n}\times\mathbb{C}^{m}
(A,ξ)\displaystyle(A^{*},\xi^{*}) (A,p1(Aξ),,pm(Aξ))\displaystyle\mapsto(A^{*},p_{1}(A^{*}\xi^{*}),\dots,p_{m}(A^{*}\xi^{*}))

is birational onto the image g(m×n×n)g(\mathbb{C}^{m\times n}\times\mathbb{C}^{n}).

When m=n+1m=n+1, the morphism gg equals ff, i.e., Theorem 2 implies Theorem 1 in this special case.

We decompose gg into the composition of two dominant morphisms α\alpha and β\beta defined by (9) and (12) respectively. To show gg is birational, it suffices to show α\alpha and β\beta are birational.

Define the morphism α\alpha

α:m×n×n\displaystyle\alpha:\mathbb{C}^{m\times n}\times\mathbb{C}^{n} Dn+1(A|y)m×n×m,\displaystyle\rightarrow D_{n+1}(A|y)\subseteq\mathbb{C}^{m\times n}\times\mathbb{C}^{m}, (9)
(A,ξ)\displaystyle(A^{*},\xi^{*}) (A,Aξ)\displaystyle\mapsto(A^{*},A^{*}\xi^{*})

where Dn+1(A|y)D_{n+1}(A|y) is the zero locus of the all (n+1)(n+1)-minors of the matrix (A|y)(A|y), A=[aij]A=[a_{ij}] is a mm-by-nn matrix of mnmn variables and y=[yi]y=[y_{i}] is a column vector of mm variables, i[m],j[n]{i\in[m],j\in[n]}. Namely, Dn+1(A|y)D_{n+1}(A|y) is the determinantal variety consisting of mm-by-(n+1)(n+1) matrices of rank less than n+1n+1. The determinantal variety Dn+1(A|y)D_{n+1}(A|y) is irreducible [6, Proposition 1.1], i.e., it is not a union of any two proper closed subvarieties.

Lemma 1.

The morphism α\alpha is birational.

Proof.

The determinantal variety Dn+1(A|y)D_{n+1}(A|y) is irreducible [6, (1.1) Proposition], and the morphism α\alpha maps the non-empty open subset S0S_{0} of m×n×n\mathbb{C}^{m\times n}\times\mathbb{C}^{n} to the non-empty open subset S1=α(S0)S_{1}=\alpha(S_{0}) in Dn+1(A|y)D_{n+1}(A|y), where

S0\displaystyle S_{0} :={(A,ξ)m×n×n:rank(A)=n}and\displaystyle:=\{(A^{*},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\xi^{*}})\in\mathbb{C}^{m\times n}\times\mathbb{C}^{n}:\mathrm{rank}(A^{*})=n\}~{}\text{and} (10)
S1\displaystyle S_{1} :={(A,η)Dn+1(A|y):rank(A)=n}\displaystyle:=\{(A^{*},\eta^{*})\in D_{n+1}(A|y):\mathrm{rank}(A^{*})=n\} (11)

are respectively the set of the solutions to USP and the set of correctly sorted sample points in USP. Hence, we deduce that α\alpha is dominant. Furthermore, for any (A,η)S1(A^{*},\eta^{*})\in S_{1}, using Cramer’s rule in linear algebra, we see that the fiber α1(A,η)\alpha^{-1}(A^{*},\eta^{*}) is a singleton, hence the morphism α\alpha is restricted to a bijection from S0S_{0} to S1S_{1}. Therefore, α\alpha is birational.

The second dominant morphism, to be shown to be birational, is

β:Dn+1(A|y)\displaystyle\beta:D_{n+1}(A|y) β(Dn+1(A|y))¯m×n×m\displaystyle\rightarrow\overline{\beta(D_{n+1}(A|y))}\subseteq\mathbb{C}^{m\times n}\times\mathbb{C}^{m} (12)
(A,η)\displaystyle(A^{*},\eta^{*}) (A,p1(η),,pm(η))\displaystyle\mapsto(A^{*},p_{1}(\eta^{*}),\dots,p_{m}(\eta^{*}))

Note that β\beta is the restriction of the finite morphism

m×n×m\displaystyle\mathbb{C}^{m\times n}\times\mathbb{C}^{m} m×n×m,\displaystyle\rightarrow\mathbb{C}^{m\times n}\times\mathbb{C}^{m}, (13)
(A,η)\displaystyle(A^{*},\eta^{*}) (A,p1(η),,pm(η))\displaystyle\mapsto(A^{*},p_{1}(\eta^{*}),\dots,p_{m}(\eta^{*}))

hence β\beta is also a finite and closed morphism [22, Chapter II, Exercise 3.5], and

β(Dn+1(A|y))=β(Dn+1(A|y))¯=g(m×n×n)¯.\beta(D_{n+1}(A|y))=\overline{\beta(D_{n+1}(A|y))}=\overline{g(\mathbb{C}^{m\times n}\times\mathbb{C}^{n})}. (14)

To prove the following Lemma 2, we shall first clarify the actions of Σm\Sigma_{m} on both m×n×m\mathbb{C}^{m\times n}\times\mathbb{C}^{m} and [A,y]\mathbb{C}[A,y].

A permutation σΣm\sigma\in\Sigma_{m} acts on [A,y]\mathbb{C}[A,y] as a \mathbb{C}-algebra isomorphism given by σ(yi)=yσ(i)\sigma(y_{i})=y_{\sigma(i)} and fixing AA. Then for any f[A,y]f\in\mathbb{C}[A,y] and y=(y1,,ym)my^{*}=(y_{1}^{*},\dots,y_{m}^{*})\in\mathbb{C}^{m}, we have

(σ(f))(y)=f(yσ(1),,yσ(m)).(\sigma(f))(y^{*})=f\left(y_{\sigma(1)}^{*},\dots,y_{\sigma(m)}^{*}\right).

A permutation σΣm\sigma\in\Sigma_{m} acts on y=(y1,,ym)my^{*}=(y_{1}^{*},\dots,y_{m}^{*})\in\mathbb{C}^{m} by

σ(y1,,ym)=(yσ1(1),,yσ1(m)).\sigma(y_{1}^{*},\dots,y_{m}^{*})=\left(y_{\sigma^{-1}(1)}^{*},\dots,y_{\sigma^{-1}(m)}^{*}\right). (15)

Then (σ(f))(y)=f(σ1(y))(\sigma(f))(y^{*})=f\left(\sigma^{-1}(y^{*})\right). This action induces the fiberwise action of Σm\Sigma_{m} on m×n×m\mathbb{C}^{m\times n}\times\mathbb{C}^{m}, i.e., for (A,y)m×n×m(A^{*},y^{*})\in\mathbb{C}^{m\times n}\times\mathbb{C}^{m}, we define

σ(A,y)=(A,σ(y)).\sigma(A^{*},y^{*})=(A^{*},\sigma(y^{*})). (16)

For the variable vector y=[y1,,ym]y=[y_{1},\dots,y_{m}],

σ(y):=[yσ(1),,yσ(m)]\sigma(y):=[y_{\sigma(1)},\dots,y_{\sigma(m)}]

is also a vector of variables, and we denote by Dn+1(A|σ(y))D_{n+1}(A|\sigma(y)) the determinantal variety on which all the (n+1)(n+1)-minors of (A|σ(y))(A|\sigma(y)) vanish. Then we have

Dn+1(A|σ(y))=σ(Dn+1(A|y)).D_{n+1}(A|\sigma(y))=\sigma(D_{n+1}(A|y)).
Lemma 2.

The morphism β\beta is birational.

Proof.

To show that β\beta is birational, it suffices to prove the following three facts

  • The subset

    W1:=Dn+1(A|y)σΣm{1}Dn+1(A|σ(y))W_{1}:=D_{n+1}(A|y)\setminus\bigcup_{\sigma\in\Sigma_{m}\setminus\{1\}}{D_{n+1}(A|\sigma(y))} (17)

    is dense open in Dn+1(A|y)D_{n+1}(A|y).

    By definition, W1W_{1} is an open subset of Dn+1(A|y)D_{n+1}(A|y). According to [62, Proof of Theorem 1], for generic matrix Am×nA^{*}\in\mathbb{C}^{m\times n}, generic column vector ξn\xi^{*}\in\mathbb{C}^{n} and any mm-permutation σ1\sigma\neq 1, we have

    rank(A|σ(Aξ))=n+1.{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathrm{rank}(A^{*}|\sigma(A^{*}\xi^{*}))=n+1.}

    Hence W1W_{1} is nonempty and dense open in Dn+1(A|y)D_{n+1}(A|y).

  • β1(β(W1))=W1\beta^{-1}(\beta(W_{1}))=W_{1}, and the restriction β\beta^{\prime} of β\beta on W1W_{1} is a bijection from W1W_{1} to β(W1)\beta(W_{1}).

    We first notice that for any (A,η)(A^{*},\eta^{*}), all the entries in η\eta^{*} are different from each other. Using Vieta’s Theorem, we derive that for any (A,ζ)β1(β(W1))(A^{*},\zeta^{*})\in\beta^{-1}(\beta(W_{1})), there is a permutation σΣm\sigma\in\Sigma_{m} such that

    σ(ζ)=η.\sigma(\zeta^{*})=\eta^{*}.

    Thus σ=1\sigma=1 and (A,ζ)=(A,η)W1(A^{*},\zeta^{*})=(A^{*},\eta^{*})\in W_{1}. Moreover, we also deduce that β\beta^{\prime} is injective.

    It is clear that β\beta^{\prime} is surjective, so the bijectivity of β\beta^{\prime} follows.

  • β(W1)\beta(W_{1}) is dense open in β(Dn+1(A|y))\beta(D_{n+1}(A|y)).

    Vieta’s Theorem implies that β\beta is a finite morphism, hence β\beta is a closed morphism. Since W1W_{1} is open in Dn+1(A|y)D_{n+1}(A|y) and β1(β(W1))=W1\beta^{-1}(\beta(W_{1}))=W_{1}, we conclude that β(W1)\beta(W_{1}) is also open in β(Dn+1(A|y))\beta(D_{n+1}(A|y)).

Now Theorem 2 follows, because gg is the composition of the birational morphisms α\alpha and β\beta.

Remark 2.

In this remark, we construct explicitly the open subset WW of 𝒳\mathcal{X} such that if the sample point (A,y)(A^{*},y^{*}) of USP lies in WW, then there exists a unique mm-permutation π\pi and a unique vector ξn\xi^{*}\in\mathbb{C}^{n} such that Aξ=π(y)A^{*}\xi^{*}=\pi(y^{*}).

  • The set SS of sample points (A,y)(A^{*},y^{*}) in USP is

    S:=σΣmσ(S1)m×n×m,S:=\bigcup_{\sigma\in\Sigma_{m}}{\sigma(S_{1})}\subseteq\mathbb{C}^{m\times n}\times\mathbb{C}^{m}, (18)

    where S1S_{1} is defined in (11), and the action of Σm\Sigma_{m} on m×n×m\mathbb{C}^{m\times n}\times\mathbb{C}^{m} is defined by (15) and (16). SS is a dense open subset of the Σm\Sigma_{m}-equivariant closed subvariety

    𝒳:=σΣmσ(Dn+1(A|y))m×n×m,\mathcal{X}:=\bigcup_{\sigma\in\Sigma_{m}}{\sigma(D_{n+1}(A|y))}\subseteq\mathbb{C}^{m\times n}\times\mathbb{C}^{m}, (19)

    A key observation is that, the finite morphism in (13) parameterizes the Σm\Sigma_{m}-orbits of 𝒳\mathcal{X}.

  • Since W1W_{1} in (17) is nonempty, for any τΣm\tau\in\Sigma_{m} the set

    τ(W1)=τ(Dn+1(A|y))σΣm{τ}σ(Dn+1(A|y)))\tau(W_{1})=\tau(D_{n+1}(A|y))\setminus\bigcup_{\sigma\in\Sigma_{m}\setminus\{\tau\}}{\sigma(D_{n+1}(A|y))})

    is also nonempty and dense open in τ(Dn+1(A|y))\tau(D_{n+1}(A|y)). Therefore, the variety 𝒳\mathcal{X} has exactly |Σm|=m!\lvert\Sigma_{m}\rvert=m! irreducible components σ(Dn+1(A|y))\sigma(D_{n+1}(A|y)) for σΣm\sigma\in\Sigma_{m}.

  • Write W1W_{1} in (17), then the set

    W:=σΣmσ(W1),{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W:=\bigcup_{\sigma\in\Sigma_{m}}{\sigma(W_{1})},}

    is dense open in 𝒳\mathcal{X}, hence SWS\cap W\neq\varnothing. Since any point in WW lies in a unique irreducible component of 𝒳\mathcal{X}, we conclude that if the sample point (A,y)(A^{*},y^{*}) of USP lies in WW, then there exists a unique mm-permutation π\pi and a unique vector ξn\xi^{*}\in\mathbb{C}^{n} such that Aξ=π(y)A^{*}\xi^{*}=\pi(y^{*}). Therefore, the solution (A,ξ)(A^{*},\xi^{*}) to USP is unique for most sample points, and in this case (A,ξ)(A^{*},\xi^{*}) is exactly the solution of system QmQ_{m} in (2).

2.2. The morphism γ\gamma is birational

To show the morphism γ\gamma in (7) is birational, we need the projection

δ:f(m×n×n)¯\displaystyle\delta:\overline{f(\mathbb{C}^{m\times n}\times\mathbb{C}^{n})} m×n×n.\displaystyle\rightarrow\mathbb{C}^{m\times n}\times\mathbb{C}^{n}. (20)
(A,p1(η),,pn+1(η))\displaystyle(A^{*},p_{1}(\eta^{*}),\dots,p_{n+1}(\eta^{*})) (A,p1(η),,pn(η))\displaystyle\mapsto(A^{*},p_{1}(\eta^{*}),\dots,p_{n}(\eta^{*})) (21)

Notice that (A,p1(η),,pn+1(η))(A^{*},p_{1}(\eta^{*}),\dots,p_{n+1}(\eta^{*})) is a point in f(m×n×n)f(\mathbb{C}^{m\times n}\times\mathbb{C}^{n}). Hence, δ\delta is uniquely determined by (21).

Now we give the sketch of the proof that the morphism γ\gamma is birational.

  • For generic (A,w)m×n×n(A^{*},w^{*})\in\mathbb{C}^{m\times n}\times\mathbb{C}^{n}, the fiber (δγg)1(A,w)(\delta\gamma g)^{-1}(A^{*},w^{*}) consists of at most n!n! points [62, Theorem 2] . Since gg is birational, we deduce that the fiber (δγ)1(A,w)(\delta\gamma)^{-1}(A^{*},w^{*}) also generically consists of at most n!n! points. Hence, it suffices to show that δ\delta is also a generically finite morphism with generic covering degree n!n!.

  • Let A=[aij]m×n,t=[t1,,tn]A=[a_{ij}]_{m\times n},t=[t_{1},\dots,t_{n}] and z=[z1,,zn+1]z=[z_{1},\dots,z_{n+1}] be the matrix and vectors in variables aija_{ij}, ti,zjt_{i},z_{j} respectively. Let R(A,z)[A,z]R(A,z)\in\mathbb{Z}[A,z] and c(A)[A]c(A)\in\mathbb{Z}[A] be the resultants of the polynomial systems

    {pk(At)zk:k=1,,n+1},\displaystyle\{p_{k}(At)-z_{k}:k=1,\dots,n+1\}, (22)
    {pk(At):k=1,,n}\displaystyle\{p_{k}(At):k=1,\dots,n\} (23)

    eliminating the variables t1,,tnt_{1},\ldots,t_{n}, respectively. We shall show that

    1. (i)

      The leading term of R(A,z)R(A,z) in zn+1z_{n+1} is c(A)n+1zn+1n!c(A)^{n+1}z_{n+1}^{n!};

    2. (ii)

      R(A,p1(Ax),,pn+1(Ax))=0R(A,p_{1}(Ax),\dots,p_{n+1}(Ax))=0;

    3. (iii)

      R(A,z)R(A,z) is irreducible in (A,z1,,zn)[zn+1]\mathbb{C}(A,z_{1},\dots,z_{n})[z_{n+1}].

    • To show (i), we use Macaulay’s assertions on resultants [40, Section 8] to deduce that degzn+1(R(A,z))n!\deg_{z_{n+1}}(R(A,z))\leqslant n! and the coefficient of zn+1n!z_{n+1}^{n!} in R(A,z)R(A,z) is c(A)n+1c(A)^{n+1}. According to Lemma 4 in [62], p1(At),,pn(At)p_{1}(A^{*}t),\dots,p_{n}(A^{*}t) is a homogeneous regular sequence of [t]\mathbb{C}[t] for generic matrix Am×nA^{*}\in\mathbb{C}^{m\times n}, hence p1(At),,pn(At)p_{1}(At),\dots,p_{n}(At) is also a homogeneous (A)[t]\mathbb{C}(A)[t]-regular sequence. Thus, the system (23) in variables tt has no solution in the projective (n1)(n-1)-space over the algebraic closure of (A)\mathbb{C}(A). Therefore, [40, Section 10] asserts that c(A)c(A) is nonzero in [A]\mathbb{C}[A], whence the leading term of R(A,z)R(A,z) in zn+1z_{n+1} is c(A)n+1zn+1n!c(A)^{n+1}z_{n+1}^{n!}.

      Moreover, the fact that p1(At),,pn(At)p_{1}(At),\dots,p_{n}(At) is a homogeneous (A)[t]\mathbb{C}(A)[t]-regular sequence implies that the field extension degree

      [(A,x):(A,p1(Ax),,pn(Ax))]=n![\mathbb{C}(A,x):\mathbb{C}(A,p_{1}(Ax),\dots,p_{n}(Ax))]=n!

      is the product of the degrees of pk(At),k[n]p_{k}(At),k\in[n] in tt. This fact also implies that for generic (A,w)m×n×n(A^{*},w^{*})\in\mathbb{C}^{m\times n}\times\mathbb{C}^{n}, the fiber (δγ)1(A,w)(\delta\gamma)^{-1}(A^{*},w^{*}) consists of n!n! points.

    • To show (ii), we notice that the system (22) substituting zkz_{k} with pk(Ax)p_{k}(Ax) has a solution t=xt=x in the affine space over the algebraic closure of (A,p1(Ax),,pn+1(Ax))\mathbb{C}(A,p_{1}(Ax),\dots,p_{n+1}(Ax)). According to arguments in [40, Section 10], we know that

      R(A,p1(Ax),,pn+1(Ax))=0.R(A,p_{1}(Ax),\dots,p_{n+1}(Ax))=0.
    • To show (iii), let A¯\overline{A} be the first n+1n+1 rows of AA, and A>n+1A_{>n+1} be the last m(n+1)m-(n+1) rows of AA. Then

      1. (a)

        R(A¯,z)R(A,z)modA>n+1R\left(\overline{A},z\right)\equiv R(A,z)\mod A_{>n+1}, and c(A¯)0c\left(\overline{A}\right)\neq 0;

      2. (b)

        The leading term of R(A¯,z)R\left(\overline{A},z\right) in zn+1z_{n+1} is c(A¯)n+1zn+1n!c\left(\overline{A}\right)^{n+1}z_{n+1}^{n!};

      3. (c)

        R(A¯,p1(A¯x),,pn+1(A¯x))=0R\left(\overline{A},p_{1}\left(\overline{A}x\right),\dots,p_{n+1}\left(\overline{A}x\right)\right)=0;

      4. (d)

        [(A¯,x):(A¯,p1(A¯x),,pn(A¯x))]=n!\left[\mathbb{C}\left(\overline{A},x\right):\mathbb{C}\left(\overline{A},p_{1}\left(\overline{A}x\right),\dots,p_{n}\left(\overline{A}x\right)\right)\right]=n!.

      Note that the birational morphism gg for the case m=n+1m=n+1 implies that

      (A¯,p1(A¯x),,pn+1(A¯x))=(A¯,x),\mathbb{C}\left(\overline{A},p_{1}\left(\overline{A}x\right),\dots,p_{n+1}\left(\overline{A}x\right)\right)=\mathbb{C}\left(\overline{A},x\right),

      hence the polynomial R(A¯,p1(A¯x),,pn(A¯x),zn+1)R\left(\overline{A},p_{1}\left(\overline{A}x\right),\dots,p_{n}\left(\overline{A}x\right),z_{n+1}\right) in zn+1z_{n+1} is irreducible over (A¯,p1(A¯x),,pn(A¯x))\mathbb{C}\left(\overline{A},p_{1}\left(\overline{A}x\right),\dots,p_{n}\left(\overline{A}x\right)\right). Since p1(A¯x),,pn(A¯x)p_{1}\left(\overline{A}x\right),\dots,p_{n}\left(\overline{A}x\right) is (A¯)[x]\mathbb{C}\left(\overline{A}\right)[x]-regular, R(A¯,z)R\left(\overline{A},z\right) in zn+1z_{n+1} is also irreducible over (A¯,z1,,zn)\mathbb{C}\left(\overline{A},z_{1},\dots,z_{n}\right). Using Gauss’ Lemma, we see that any factorization R(A,z)=uvR(A,z)=uv in (A,z1,,zn)[zn+1]\mathbb{C}(A,z_{1},\dots,z_{n})[z_{n+1}] induces a factorization R(A,z)=u1v1R(A,z)=u_{1}v_{1} in [A,z]\mathbb{C}[A,z] satisfying

      degzn+1u1=degzn+1u,degzn+1v1=degzn+1v\deg_{z_{n+1}}u_{1}=\deg_{z_{n+1}}u,~{}\deg_{z_{n+1}}v_{1}=\deg_{z_{n+1}}v

      and the leading coefficients of u1u_{1} and v1v_{1} in zn+1z_{n+1} ly in [A]\mathbb{C}[A]. Modulo A>n+1A_{>n+1} and noticing that c(A¯)n+10c\left(\overline{A}\right)^{n+1}\neq 0 is the leading coefficient of R(A¯,z)R\left(\overline{A},z\right), we deduce that R(A¯,z)=u1¯v1¯R\left(\overline{A},z\right)=\overline{u_{1}}\overline{v_{1}}, degzn+1u1¯=degzn+1u1\deg_{z_{n+1}}\overline{u_{1}}=\deg_{z_{n+1}}u_{1} and degzn+1v1¯=degzn+1v1\deg_{z_{n+1}}\overline{v_{1}}=\deg_{z_{n+1}}v_{1}. Thus, R(A,z)R(A,z) is irreducible over (A,z1,,zn)\mathbb{C}(A,z_{1},\dots,z_{n}).

    Hence, the algebraic field extension degree

    [(A,p1(Ax),,pn+1(Ax)):(A,p1(Ax),,pn(Ax))]=n![\mathbb{C}(A,p_{1}(Ax),\dots,p_{n+1}(Ax)):\mathbb{C}(A,p_{1}(Ax),\dots,p_{n}(Ax))]=n!

    coincides with

    [(A,x):(A,p1(Ax),,pn(Ax))]=n!,[\mathbb{C}(A,x):\mathbb{C}(A,p_{1}(Ax),\dots,p_{n}(Ax))]=n!,

    which implies that δ\delta is a generically finite morphism with a generic covering degree n!n!, and

    (A,p1(Ax),,pn+1(Ax))=(A,x),\mathbb{C}(A,p_{1}(Ax),\dots,p_{n+1}(Ax))=\mathbb{C}(A,x),

    namely, γ\gamma is birational, and Theorem 1 follows. We refer to [37] for detailed proof of the above arguments.

Remark 3.

From the fact that γ\gamma is birational, we can also derive that for generic (A,w)m×n×n(A^{*},w^{*})\in\mathbb{C}^{m\times n}\times\mathbb{C}^{n}, the fiber (δγg)1(A,w)(\delta\gamma g)^{-1}(A^{*},w^{*}) exactly consists of n!n! points, among which there is, however, only one solution to USP. This result indicates that the bound n!n! in [62, Theorem 2] is optimal.

3. Rank-one moment matrix completion

There are many well-known symbolic or numeric algorithms and software packages for solving zero-dimensional polynomial systems in variables x=[x1,,xn]x=[x_{1},\ldots,x_{n}]

{g1(x)=0,,gm(x)=0},\{g_{1}(x)=0,\ldots,g_{m}(x)=0\}, (24)

e.g. [19, 4, 66, 3, 34, 54]. A desirable method for us is the semidefinite relaxation method proposed by Parrilo [50, 51], Lasserre [33] for solving a sequence of SDP problems with constant objective function 11

{min1s. t.𝗒0=1,Mt(𝗒)𝟢,Mtdj(gj𝗒)=0,j=1,,m,\left\{\begin{array}[]{cl}\displaystyle\min&1\\ \text{s.\ t.}&{\sf y}_{0}=1,\\ &M_{t}(\sf{y})\succeq 0,\\ &M_{t-d_{j}}(g_{j}{\sf y})=0,\quad{j=1,\ldots,m,}\\ \end{array}\right.\quad (25)

where dj:=deg(gj)/2d_{j}:=\lceil\deg(g_{j})/2\rceil, 𝗒=(𝗒α)αnn{\sf y}={({\sf y}_{\alpha})}_{\alpha\in\mathbb{N}^{n}}\in\mathbb{R}^{\mathbb{N}^{n}}, and Mt(𝗒)M_{t}({\sf y}) is the truncated moment matrix of order tt

Mt(𝗒):=(𝗒α+β)|α|,|β|tn×n,M_{t}({\sf y}):=({\sf y}_{\alpha+\beta})_{|\alpha|,|\beta|\leqslant t}\in\mathbb{R}^{\mathbb{N}^{n}\times\mathbb{N}^{n}},

|α|=iαi|\alpha|=\sum_{i}\alpha_{i}, |β|=iβi,|\beta|=\sum_{i}\beta_{i}, which is the principal submatrix of the full (infinite) moment matrix

M(𝗒)=(𝗒α+β)n×n,α,βn.M({\sf y})=({\sf y}_{\alpha+\beta})\in\mathbb{R}^{\mathbb{N}^{n}\times\mathbb{N}^{n}},~{}\alpha,\beta\in\mathbb{N}^{n}.

Suppose gj(x)=Σαngj,αxα[x]g_{j}(x)=\Sigma_{\alpha\in\mathbb{N}^{n}}g_{j,\alpha}x^{\alpha}\in\mathbb{R}[x]. If the (k,l)(k,l)-th entry of Mt(𝗒)M_{t}({\sf y}) is yβy_{\beta}, then the (k,l)(k,l)-th entry of the localizing matrix Mtdj(gj𝗒)M_{t-d_{j}}(g_{j}{\sf y}) with respect to 𝗒{\sf y} and gjg_{j} is defined by

Mtdj(gj𝗒)(k,l):=αgj,α𝗒α+β.\displaystyle M_{t-d_{j}}(g_{j}{\sf y})(k,l):=\sum_{\alpha}{g_{j,\alpha}{\sf y}_{\alpha+\beta}}.

Lasserre, Laurent, and Rostalski gave explicitly rank conditions that guarantee us to find all real solutions of the polynomial system (24) [30, 31, 32] by the semidefinite relaxation method. Let d=maxj=1,,mdj,d=\max_{j=1,\ldots,m}d_{j}, for tdt\geqslant d, define the convex set

Kt:={𝗒2tn|𝗒0=1,Mt(𝗒)0,Mtdj(gj𝗒)=0,j=1,,m}K_{t}^{\mathbb{R}}:=\left\{{\sf y}\in\mathbb{R}^{\mathbb{N}_{2t}^{n}}|{\sf y}_{0}=1,M_{t}({\sf y})\succeq 0,M_{t-d_{j}}(g_{j}{\sf y})=0,j=1,\ldots,m\right\} (26)
K:={𝗒n|𝗒0=1,M(𝗒)0,M(gj𝗒)=0,j=1,,m}K^{\mathbb{R}}:=\left\{{\sf y}\in\mathbb{R}^{\mathbb{N}^{n}}|{\sf y}_{0}=1,M({\sf y})\succeq 0,M(g_{j}{\sf y})=0,j=1,\ldots,m\right\} (27)

We introduce the following notations: Tn:={xα|αn}T^{n}:=\{x^{\alpha}|\alpha\in\mathbb{N}^{n}\}, Ttn:={xαTn||α|t}T^{n}_{t}:=\{x^{\alpha}\in T^{n}||\alpha|\leq t\}.

Proposition 1.

Let I=g1,,gmI=\langle g_{1},\dots,g_{m}\rangle be an ideal in [x],V(I)Ø\mathbb{R}[x],V_{\mathbb{R}}(I)\neq\O. If vnv\in\mathbb{R}^{n} is a real root of II, i.e. vV(I)v\in V_{\mathbb{R}}(I), then for tdt\geq d,

𝗒tδv=(vα)αT2tn{\sf y}^{\delta_{v}}_{t}=(v^{\alpha})_{\alpha\in T^{n}_{2t}}

is an optimal solution for the optimization problem:

{minrank(Mt(𝗒))s. t.𝗒0=1,Mt(𝗒)𝟢,Mtdj(gj𝗒)=0,j=1,,m,\left\{\begin{array}[]{cl}\displaystyle\min&\mathrm{rank}(M_{t}({\sf y}))\\ \text{s.\ t.}&{\sf y}_{0}=1,\\ &M_{t}(\sf{y})\succeq 0,\\ &M_{t-d_{j}}(g_{j}{\sf y})=0,\quad{j=1,\ldots,m,}\\ \end{array}\right.\quad (28)
Proof.

Since 𝗒tδv{\sf y}^{\delta_{v}}_{t} comes from a Dirac measure δv\delta_{v}, by definition (𝗒tδv)0=v0=1({\sf y}^{\delta_{v}}_{t})_{0}=v^{0}=1 and

Mt(𝗒tδv)=ζv,tζv,t0,M_{t}({\sf y}^{\delta_{v}}_{t})=\zeta_{v,t}\zeta_{v,t}^{\top}\succeq 0,

where ζv,t=(vα)αTtn\zeta_{v,t}=(v^{\alpha})_{\alpha\in T_{t}^{n}} is a column vector. For gj(x)=αTngj,αxαg_{j}(x)=\sum_{\alpha\in T^{n}}g_{j,\alpha}x^{\alpha}, define vec(gj)=(gj,α)αvec(g_{j})=(g_{j,\alpha})_{\alpha} Therefore vec(gj)Mt(𝗒tδv)vec(gj)=gj(v)2=0vec(g_{j})^{\top}M_{t}({\sf y}^{\delta_{v}}_{t})vec(g_{j})=g_{j}(v)^{2}=0, combining Mt(𝗒tδv)0M_{t}({\sf y}^{\delta_{v}}_{t})\succeq 0, it follows that

Mt(𝗒tδv)vec(gj)=0,vec(gj)kerMt(𝗒tδv),Mtdj(gj𝗒)=0.M_{t}({\sf y}^{\delta_{v}}_{t})vec(g_{j})=0,vec(g_{j})\in\ker M_{t}({\sf y}^{\delta_{v}}_{t}),M_{t-d_{j}}(g_{j}{\sf y})=0. (29)

We have shown that 𝗒tδvKt{\sf y}^{\delta_{v}}_{t}\in K_{t}^{\mathbb{R}} is a feasible solution of (28) and rankMt(𝗒tδv)=1\mathrm{rank}M_{t}({\sf y}^{\delta_{v}}_{t})=1. For 𝗒Kt,𝗒0=1rank(Mt(𝗒))1\forall{\sf y}\in K_{t}^{\mathbb{R}},{\sf y}_{0}=1\Rightarrow\mathrm{rank}(M_{t}({\sf y}))\geq 1, so the optimal value of (28) is 11, which can be achieved at 𝗒tδv{\sf y}^{\delta_{v}}_{t}. ∎

Proposition 2.

Let I=g1,,gmI=\langle g_{1},\dots,g_{m}\rangle be an ideal in [x],V(I)Ø\mathbb{R}[x],V_{\mathbb{R}}(I)\neq\O. For tdt\geq d, suppose 𝗒^tKt\hat{{\sf y}}_{t}\in K_{t}^{\mathbb{R}} is an optimal solution of the optimization problem (28). Then there vV(I)\exists v\in V_{\mathbb{R}}(I), s.t. 𝗒^t=𝗒tδv\hat{{\sf y}}_{t}={\sf y}^{\delta_{v}}_{t}.

Proof.

Because V(I)ØV_{\mathbb{R}}(I)\neq\O, choose some ξ0V(I)\xi_{0}\in V_{\mathbb{R}}(I) and construct ytδξ0y^{\delta_{\xi_{0}}}_{t}. According to Proposition 1, the optimal value of (28) is 11 and can be achieved at ytδξ0y^{\delta_{\xi_{0}}}_{t}. Since 𝗒^tKt\hat{{\sf y}}_{t}\in K_{t}^{\mathbb{R}} is an optimal solution of (28), rank(Mt(𝗒^t))=1\mathrm{rank}(M_{t}(\hat{{\sf y}}_{t}))=1. Notice that M0(𝗒^t)=(1)M_{0}(\hat{{\sf y}}_{t})=(1),

1=rank(M0(𝗒^t))rank(Ms(𝗒^t))rank(Mt(𝗒^t))=1,0st.1=\mathrm{rank}(M_{0}(\hat{{\sf y}}_{t}))\leq\mathrm{rank}(M_{s}(\hat{{\sf y}}_{t}))\leq\mathrm{rank}(M_{t}(\hat{{\sf y}}_{t}))=1,0\leq s\leq t.

Mt(𝗒^t)M_{t}(\hat{{\sf y}}_{t}) satisfies flat extension condition. Hence 𝗒^t\hat{{\sf y}}_{t} can be uniquely extended to 𝗒~K\tilde{{\sf y}}\in K^{\mathbb{R}}, s.t. rank(M(𝗒~))=rank(Mt(𝗒^t))=1\mathrm{rank}(M(\tilde{{\sf y}}))=\mathrm{rank}(M_{t}(\hat{{\sf y}}_{t}))=1. Then Theorem 3.3 in [30] shows that there exists vnv\in\mathbb{R}^{n}, 𝗒=𝗒δv=(vα)αTn{\sf y}={\sf y}^{\delta_{v}}=(v^{\alpha})_{\alpha\in T^{n}}, with ker(M(𝗒~))=I({v})\ker(M(\tilde{{\sf y}}))=I(\{v\}). Then M(gj𝗒~)=0gj𝗒~=M(y~)vec(gj)=0M(g_{j}\tilde{{\sf y}})=0\Rightarrow g_{j}\tilde{{\sf y}}=M(\tilde{y})vec(g_{j})=0, i.e. vec(gj)ker(M(𝗒~))vec(g_{j})\in\ker(M(\tilde{{\sf y}})), so gj(v)=0,j=1,mg_{j}(v)=0,j=1,\dots m, which implies that vV(I)v\in V_{\mathbb{R}}(I).

Combining the above results, we have proved that:

Theorem 3.

Suppose V(I)Ø,tdV_{\mathbb{R}}(I)\neq\O,t\geq d then the map

V(I)\displaystyle V_{\mathbb{R}}(I) Opt(I):={𝗒Kt|𝗒 is an optimal solution of (28)}\displaystyle\to Opt(I):=\{{\sf y}\in K_{t}^{\mathbb{R}}~{}|~{}{\sf y}\text{ is an optimal solution of (\ref{sdp*})}\}
v\displaystyle v 𝗒tδv\displaystyle\mapsto{\sf y}^{\delta_{v}}_{t}

gives an 1:11:1 correspondence between a real root in V(I)V_{\mathbb{R}}(I) and an optimal solution of (28)(which implies that the rank of the moment matrix is 11).

According to the above Theorem, to obtain the unique solution of II, one must find the truncated moment sequence 𝗒Kd{\sf y}\in K_{d}^{\mathbb{R}} with rank 11. Hence, finding the unique real solution of a polynomial system is reduced to solving the rank-one moment matrix completion problem. There are many well-known methods for solving low-rank matrix completion problems [10, 11, 21, 53, 12, 65, 25, 47, 14, 15, 23, 48, 49, 16].

The nuclear norm of a matrix is the sum of the singular values of the matrix. If the matrix is symmetric and positive semidefinite, then its nuclear norm is equal to the trace of the matrix. In [53, Theorem 2.2], Recht, Fazel, and Parrilo have shown that the nuclear norm is the best convex lower approximation of the rank function over the set of matrices with spectral norm less than or equal to 11. Therefore, the problem of finding a rank-one moment matrix can be relaxed to the following form:

{mintr(Mt(𝗒))s. t.𝗒0=1,Mt(𝗒)0,Mtdj(gj𝗒)=0,j=1,,m\left\{\begin{array}[]{cl}\displaystyle\min&\mathrm{tr}(M_{t}({\sf y}))\\ \text{s.\ t.}&{\sf y}_{0}=1,\\ &M_{t}({\sf y})\succeq 0,\\ &M_{t-d_{j}}(g_{j}{\sf y})=0,\quad{j=1,\ldots,m}\\ \end{array}\right. (30)

In [14, 15], Cossc and Demanet also show that when the rank-one matrix completion problem has a unique solution, the second-order semidefinite relaxation with minimization of the nuclear norm will be enough for finding the unique solution. In [39, 67, 68], we have shown that by solving the nuclear norm optimization problem (30), one can find some real solutions (not all solutions) for polynomial systems more efficiently.

According to Theorem 1, the polynomial system Qn+1Q_{n+1} has a unique real solution. Hence, in the next section, instead of minimizing the constant objective function 11 in (25), we solve the nuclear norm minimization problem (30) to find the unique real solution of Qn+1Q_{n+1}.

4. Experiments

4.1. Algorithm Design

For simplicity, we assume that our data Am×n,ynA^{*}\in\mathbb{R}^{m\times n},y^{*}\in\mathbb{R}^{n} are generic in the sense of Theorem 1 and [62, Theorem 2]. This section presents two symbolic algorithms and one symbolic-numeric algorithm to solve the unlabeled sensing problem (1).

A straightforward attempt to recover ξ\xi^{*} is to solve the polynomial system QnQ_{n} using Groebner basis [7, 8, 9], which is formalized as the following algorithm. The experiment 4.2.1 was conducted on Maple’s symbolic computation software, and Maple automatically chose the monomial order, usually the reverse-graded lexicographic order.

Input: Am×n,ynA^{*}\in\mathbb{Q}^{m\times n},y^{*}\in\mathbb{Q}^{n} such that πΣm\exists\pi\in\Sigma_{m}, the permutation group, and π(y)im(A)\pi(y^{*})\in\mathrm{im}(A^{*}), the image space of AA^{*}
Output: ξ\xi^{*} such that Aξ=π(y)A^{*}\xi^{*}=\pi(y^{*})
compute Groebner basis of the ideal In=(Qn)I_{n}=(Q_{n}), and extracted the roots symbolically
roots{ξn|qi(ξ)=0,i=1,,n}roots\leftarrow\{\xi\in\mathbb{Q}^{n}|q_{i}(\xi)=0,i=1,\dots,n\};
for ξroots\xi\in roots do
      if  coordinates of AξA^{*}\xi are a permutation of yy^{*} then
            return ξξ\xi^{*}\leftarrow\xi
       end if
      
end for
Algorithm 1 Groebner basis for solving square system QnQ_{n} (4)

The complexity of solving zero-dimensional polynomial systems in nn variables is known to be single exponential in nn [13, 28, 19], even for well-behaved cases [20]. The complexity of examining whether coordinates of AξA^{*}\xi are a permutation of yy^{*} is quasi-linear in mm using a quick sorting algorithm.

According to Theorem 1, the polynomial system Qn+1Q_{n+1} has only one solution. Experimental results in Section 4.2 show that it is far more efficient to compute the Groebner basis for polynomial systems with only one solution. This observation leads to the following algorithm 2.

Input: Am×n,ynA^{*}\in\mathbb{Q}^{m\times n},y^{*}\in\mathbb{Q}^{n} such that πΣm\exists\pi\in\Sigma_{m}, the permutation group, and π(y)im(A)\pi(y^{*})\in\mathrm{im}(A^{*}), the image space of AA^{*}
Output: ξ\xi^{*} such that Aξ=π(y)A^{*}\xi^{*}=\pi(y^{*})
compute Groebner basis of the ideal In+1=(Qn+1)I_{n+1}=(Q_{n+1}), and extracted the roots symbolically
roots{ξn|qi(ξ)=0,i=1,,n+1}roots\leftarrow\{\xi\in\mathbb{Q}^{n}|q_{i}(\xi)=0,i=1,\dots,n+1\};
roots={ξ}roots=\{\xi^{*}\};
return ξ\xi^{*}
Algorithm 2 Groebner basis for solving overdetermined system Qn+1Q_{n+1} (3)

Both algorithms above use symbolic computation. Given precise input and clean measurement, the solution one obtains using these two algorithms is exact. When the given data are corrupted by noise, the perturbed polynomial system Qn+1Q_{n+1} will have no solution as it is an overdetermined polynomial system. However, in real-world applications, all measurements are inevitably influenced by noise. Therefore, a robust and efficient symbolic-numerical algorithm is needed to deal with corrupted measurements yy^{*} one observes in the unlabeled sensing problem. In Section 3, we have shown that the unique solution of the polynomial system Qn+1Q_{n+1} can be obtained by solving the rank-one moment matrix completion problem (30). Below, we develop the rank-one moment matrix completion method for solving polynomial system Qn+1Q_{n+1} (3) numerically.

Input: Am×n,ynA^{*}\in\mathbb{R}^{m\times n},y^{*}\in\mathbb{R}^{n} such that πΣm,π(y)im(A)\exists\pi\in\Sigma_{m},\pi(y^{*})\in im(A^{*})
Output: ξ\xi^{*} such that Aξ=π(y)A^{*}\xi^{*}=\pi(y^{*})
solve SDP problem (30) for gi=qig_{i}=q_{i} to the relaxation order s=(n+1)/2s=\lceil(n+1)/2\rceil to get the moment matrix MsM_{s}
compute singular value decomposition of MsM_{s} [U,S,V]svd(Ms)[U,S,V]\leftarrow svd(M_{s});
select the first n+1 rows of UU
x1U(1:n+1,1)x_{1}\leftarrow U(1:n+1,1);
normalize x1x_{1} to recover ξsdp\xi_{sdp} ξsdp(ξ(k)x1(k+1,1)/x1(1,1))k=1,,n\xi_{sdp}\leftarrow(\xi(k)\leftarrow x_{1}(k+1,1)/x_{1}(1,1))_{k=1,\dots,n};
use EM method in [60] to refine ξsdp\xi_{sdp}
ξEMEM(ξsdp)\xi_{EM}\leftarrow EM(\xi_{sdp});
return ξξEM\xi^{*}\leftarrow\xi_{EM}
Algorithm 3 Matrix completion solving overdetermined system Qn+1Q_{n+1} (3)

Semidefinite programming problems can be solved using the interior-point method [5], whose complexity is bounded by a polynomial in the matrix size. There are well-developed software packages to solve SDP problems, for example, Sedumi [58] and SDPNAL+ [59].

The Expectation-Maximization (EM) method proposed in [62] is an algebraic trick to refine the solution extracted from a numerical polynomial system solver. By sorting AξsdpA^{*}\xi_{sdp} and yy^{*}, if the solution extracted is accurate enough, one expects that AξsdpA^{*}\xi_{sdp} and yy^{*} after sorting should match, hence recovers the permutation π\pi and reduces the problem to classical linear regression. The numerical precision of classical linear regression is determined by the condition number of AA^{*} and precision of yy^{*}, which, in practice, can be better controlled in numerical algorithms.

4.2. Symbolic and Numerical Experimental Results

The experiments were conducted on a desktop with an Intel(R) Core(TM) i9-10900X CPU @ 3.70GHz and 128 GB RAM. The matrix entries of AA^{*} were randomly generated from standard Gaussian distribution in these experiments, and entries ξ\xi^{*} were generated from a uniform distribution in (0,1)(0,1). Equations (qi=0)i=1,,n+1(q_{i}=0)_{i=1,\dots,n+1} were then constructed from AA^{*}, ξ\xi^{*}, and possible noise. SDP in algorithm 3 was solved using Sedumi [58]. Source codes are uploaded on GitHub at unlabeled-sensing-ISSAC-2024.

4.2.1. symbolic algorithm 1 and 2, clean input data in \mathbb{Q}, algorithm performance on different nn

Experiment 4.2.1 compares the running time of algorithm 1 (TA1/s) and algorithm 2 (TA2/s) with different n=3,4,8n=3,4\dots,8, and mm was set as 2n2n, clean measurement. Algorithm 1 uses the square system Qn={q1(x)=0,,qn(x)=0}Q_{n}=\{q_{1}(x)=0,\ldots,q_{n}(x)=0\}. Algorithm 2 uses the overdeterimined system Qn+1={q1(x)=0,,qn(x),qn+1(x)=0}Q_{n+1}=\{q_{1}(x)=0,\ldots,q_{n}(x),q_{n+1}(x)=0\}. All of the Groebner basis computations are done in Maple.

Table 1. cpu time for Algorithm 1 and Algorithm 2, with different nn, m=2nm=2n.
nn TA1/TA1/s TA2/TA2/s
3 0.05 0.08
4 1.95 0.24
5 153.05 0.56
6 3.30
7 200.39
8 72549.63

Table 1 shows that as the number of variables and the degree of the polynomial system increase, computing a Groebner basis for the polynomial system QnQ_{n} or Qn+1Q_{n+1} becomes more difficult. As Qn+1Q_{n+1} has only one solution, its Groebner basis in lexicographic order consists of nn linear polynomials x1ξ1,,xnξnx_{1}-\xi^{*}_{1},\ldots,x_{n}-\xi^{*}_{n}. On the other hand, the polynomial system QnQ_{n} has at most n!n! solutions; its Groebner basis consists of a set of polynomials with higher degrees and larger coefficients. Therefore, solving the overdetermined polynomial system Qn+1Q_{n+1} is much faster than solving the square polynomial system QnQ_{n}. When n6n\geqslant 6, algorithm 1 could not return the result within 24 hours and was terminated manually.

4.2.2. numerical algorithm 3, clean input data in \mathbb{R}, algorithm performance on different nn

The following experiments examine the performance of numerical algorithm 3. Since algorithm 3 computes the moment optimization problem numerically, the relative error of the solution extracted from the moment matrix is defined as errsdp:=ξξsdp2ξ2,err_{sdp}:=\frac{\lVert\xi^{*}-\xi_{sdp}\rVert_{2}}{\lVert\xi^{*}\rVert_{2}}, where ξsdp\xi_{sdp} represents the solution extracted by semidefinite programming part of algorithm 3.

Experiment 4.2.2 examines the running time of algorithm 3 (TA3/s) with different n=3,4,5,6n=3,4,5,6, and mm was set as 2n2n, tt records the relaxation order used in the semidefinite programming, ranksranks records the rank sequence of the truncated moment matrix MM, sizesize denotes the dimension of the truncated moment matrix for the given relaxation order. In this experiment, no noise was imposed on yy.

Table 2. CPU time, relative error of ξsdp\xi_{sdp}, rank sequence and moment matrix size for Algorithm 3 with different nn, m=2nm=2n.
nn tt TA3/TA3/s errsdperr_{sdp} ranksranks sizesize
33 22 0.200.20 2.73E-08 1,1,11,1,1 1010
44 33 1.881.88 6.53E-07 1,1,1,11,1,1,1 3535
55 33 5.585.58 1.62E-06 1,1,1,11,1,1,1 5656
66 44 290.95290.95 4.93E-06 1,1,1,1,11,1,1,1,1 210210

Table 2 shows that as the number of variables and the degree of polynomial rise, the time consumed to compute the moment matrix increases. The rank sequence of the moment matrix MM stabilizing at 11 indicates that the semidefinite programming reaches the flat extension criterion at the lowest relaxation order required, which is (n+1)/2\lceil(n+1)/2\rceil, therefore Theorem LABEL:cor guarantees that the solution extracted from the MM is exactly the unique solution of the polynomial system Qn+1Q_{n+1}. When using clean measurement yy, the errsdperr_{sdp} behaves well and mildly increases as nn becomes larger for n=3,4,5,6n=3,4,5,6.

4.2.3. numerical algorithm 3, corrupted measurement yy, algorithm performance on different mm

The following experiments included noise on measurement. For real solution ξ\xi^{*}, and clean measurement y:=Aξy^{*}:=A^{*}\xi^{*}, noise ycy_{c} was assumed to be a Gaussian random vector with expectation 0 and covariance matrix σ2I\sigma^{2}I, where II denotes the n×nn\times n identity matrix. The magnitude of noise in signal processing is measured using signal-noise-ratio SNRSNR in decibels. In our experiments, SNRSNR can be computed from the formula

SNR=10(log10(n3σ2)).SNR=10\left(\log_{10}\left(\frac{n}{3\sigma^{2}}\right)\right).

The equations generated using corrupted measurement y:=y+ycy:=y^{*}+y_{c} as

qic:=pi(Ax)pi(y)=0.q_{ic}:=p_{i}(A^{*}x)-p_{i}(y)=0.

We further report the relative error of the solution refined by EM algorithm [62] as errEM:=ξξEM2ξ2.err_{EM}:=\frac{\lVert\xi^{*}-\xi_{EM}\rVert_{2}}{\lVert\xi^{*}\rVert_{2}}. EM method is essentially a sorting procedure. It sorts yy according to its coordinates to ysorty_{sort}. For possible solutions (xi)i(x_{i})_{i}, compute AxiAx_{i} and sort its coordinations to (Axi)sort(Ax_{i})_{sort}. The most likely solution should be the nearest one to ysorty_{sort}. Hence, the permutation is determined accordingly. Once the permutation is found, solving USP is reduced to solving a linear system.

Experiment 4.2.3 examines the running time and relative error of algorithm 3 with respect to mm under the noise of SNR=60SNR=60dB, fixing n=4n=4.TEMT_{EM}, TsdpT_{sdp} record the CPU time of EM refinement and solving SDP problem in algorithm 3, TA3=TEM+TsdpTA3=T_{EM}+T_{sdp} is the total CPU time of algorithm 3.

Table 3. CPU time, and relative error of ξEM\xi_{EM} and ξsdp\xi_{sdp} for algorithm 3 with different m, observation y is corrupted by Gaussian random noise of SNR=60SNR=60dB, n=4n=4 fixed, 2020 trials median.
mm errEMerr_{EM} errsdperr_{sdp} TEM/T_{EM}/s Tsdp/T_{sdp}/s TA3/sTA3/s
500 0.010% 0.141% 0.003 2.477 2.450
1000 0.006% 0.102% 0.003 2.586 2.589
2000 0.005% 0.125% 0.016 2.602 2.617
5000 0.006% 0.250% 0.188 2.672 2.859
Remark 4.

Algorithm 3 may fail to find a good approximate solution using SDP solvers and EM refinement possibly cannot recover the right permutation from the solution extracted from ξsdp\xi_{sdp}, leading to the relative error of the solution returned by algorithm 3 diverging very far away from most of the trials. Here we manifest a typical outlier on the condition n=4,m=500,SNR=60n=4,m=500,SNR=60dB, the median of relative errors of the trials are errEM=0.016%err_{EM}=0.016\% and errsdp=0.396%err_{sdp}=0.396\%, the outlier appeared in our experiment is errEM=165.957%err_{EM}=165.957\% and errsdp=167.938%err_{sdp}=167.938\%. For small n=4n=4 and noise SNR=80SNR=80dB or SNR=100SNR=100dB, the proportion of the outliers is less than 5%5\%. The proportion of the outliers is between 5%5\% and 35%35\% for other nn and SNRSNR, increasing along with the number of variables and noise.

The outliers are mainly caused by noise interfering SDP algorithm. The rank sequence of SDP relaxation without noise under the assumption of a unique solution should stabilize at 11. However, in outliers, the numerical rank sequence increases beyond 1 under a certain numerical tolerance concerning the noise. Set n=4,m=50,SNR=60n=4,m=50,SNR=60dB and tolerance 0.010.01 for the numerical rank. In 5050 trials, 22 outliers were recorded; these two outliers have numerical rank sequence 1,3,5,51,3,5,5, and 1,2,3,31,2,3,3. Many factors may lead to the increase of the rank sequence. We observed that the scale of the coordinates of the real solution ξ\xi^{*} varied acutely in two outliers. The outlier with numerical rank sequence 1,3,5,51,3,5,5 corresponds to the real solution ξ=(0.753,0.081,0.326,0.879)t\xi^{*}=(0.753,0.081,0.326,0.879)^{t}; the outlier with rank sequence 1,2,3,41,2,3,4 corresponds to the real solution ξ=(0.414,0.515,0.004,0.624)t\xi^{*}=(0.414,0.515,0.004,0.624)^{t} .

Table 3 shows that the accuracy errsdperr_{sdp} and errEMerr_{EM} of algorithm 3 is stable with different mm. As mm increases, the time TA3TA3 consumed mildly rises because the EM refinement part of algorithm 3 needs to sort more entries. However, compared with the sorting procedure, the main computation in algorithm 3 occurs in solving the SDP problem, which is irrelevant to mm. Therefore, TA3TA3 increases rather slowly when mm becomes larger. Table 3 also reveals that for m=500m=500 or 10001000, the portion of time consumed on EM refinement is insignificant in the overall time consumed; hence, the overall time of algorithm 3 will be reported in further experiments, and the time of EM refinement will be omitted in tables below.

4.2.4. numerical algorithm 3, corrupted measurement yy, algorithm performance on different magnitude of noise SNRSNR, n=5,6n=5,6

Experiment 4.2.4 examines for each nn the running time and relative error of algorithm 3 with respect to the magnitude of noise, fixing m=500m=500.

Table 4. CPU time, and relative error of ξEM\xi_{EM} and ξsdp\xi_{sdp} for algorithm 3 with different magnitude of noise, n=5,m=500n=5,m=500 fixed, 2020 trials median.
SNR/SNR/dB TA3/TA3/s errsdperr_{sdp} errEMerr_{EM}
100 5.46 0.003% 0.003%
80 5.58 0.043% 0.021%
60 6.50 0.220% 0.010%
50 6.97 1.430% 0.058%
Table 5. CPU time, and relative error of ξEM\xi_{EM} and ξsdp\xi_{sdp} for algorithm 3 with different magnitude of noise, n=6,m=500n=6,m=500 fixed, 2020 trials median.
SNR/SNR/dB TA3/TA3/s errsdperr_{sdp} errEMerr_{EM}
100 185.83 0.002% 0.002%
80 184.59 0.023% 0.019%
60 253.38 0.546% 0.011%
50 220.52 1.224% 0.044%

Tables 4, 5 show that as the noise builds up, the errsdperr_{sdp} becomes larger, and EM refinement technique successfully improves the precision of the solution in all these experiments. For small noise SNR=100SNR=100dB or 8080dB, the errEMerr_{EM} is not sensitive to nn. This is not so strange because the precision of ξEM\xi_{EM} depends majorly on the sorting procedure. When the noise is relatively small, with high probability EM technique finds the right permutation π\pi. Once the permutation π\pi is deduced correctly, the precision of classical linear regression only depends on the condition number of AA and the relative noise ycAξ2\lVert\frac{y_{c}}{A^{*}\xi^{*}}\rVert_{2} imposed on the measurement. For general n=5,6n=5,6, the condition number of AA with random entries of standard normal distribution behaves well, and the magnitude of relative noise is less sensitive to nn. Hence, the relative error after EM refinement errEMerr_{EM} behaves well in small noise circumstances.

4.2.5. comparison between algorithm AIEM [62] and algorithm 3

Algorithm AIEM was introduced in [62]. AIEM uses homotopy method [43] to solve the polynomial system QnQ_{n} (4). Homotopy continuation method solves polynomial systems by transforming the original system into a continuous path from an easily solvable system to the desired one. It introduces an auxiliary parameter tt and constructs a family of homotopy equations with tt, which smoothly deforms from a known system to the target system to be solved. The homotopy solver returns at most n!n! roots, among which the real solution will be selected by EM sorting procedure. The relative error for ξhomo\xi_{homo} the solution from homotopy solver after selection by sorting is defined as: errhomo:=ξξhomo2ξ2.err_{homo}:=\frac{\lVert\xi^{*}-\xi_{homo}\rVert_{2}}{\lVert\xi^{*}\rVert_{2}}.

Table 6. CPU time and relative error for AIEM and algorithm 3 with different nn. m=500,SNR=80m=500,SNR=80dB, 20 trials median.
AIEM(homotopy) algorithm 3(SDP)
nn TAIEM/T_{AIEM}/s errhomoerr_{homo} errEMerr_{EM} TA3/TA3/s errsdperr_{sdp} errEMerr_{EM}
3 0.12 0.043% 0.012% 0.29 0.029% 0.016%
4 0.23 0.037% 0.016% 3.04 0.021% 0.019%
5 0.79 0.040% 0.014% 5.58 0.043% 0.021%
6 9.06 0.039% 0.019% 184.80 0.053% 0.011%

Table 6 shows that as nn increases, AIEM using the homotopy solver HOM4PS2 [36] is more efficient than algorithm 3 using the SDP solver Sedumi. The SDP solver Sedumi used in algorithm 3 implemented the traditional interior-point method in its inner iterations, which is numerically more stable but is relatively slow compared to modern SDP solvers such as SDPNAL+. In [62], the homotopy solver Bertini [3] was used to solve the polynomial system QnQ_{n}, for n=6n=6 AIEM needs 22432243 seconds to return the solution. Our experiments suggest that different software and implementations can influence the time and accuracy of AIEM and algorithm 3, which should be further examined and compared in the future.

4.2.6. Computational complexity analysis

The complexity of solving a zero-dimensional polynomial system using Groebner bases is single exponential dO(n)d^{O(n)} [29], where dd is the degree of polynomials and nn is the number of variables. In the context of unlabeled sensing and algorithm 1(d=nd=n), 2(d=n+1d=n+1), the overall complexity of symbolic algorithm should be nO(n)=O(exp(nlogn))n^{O(n)}=O(\exp(n\log n)).

For general input, n!\exists n! solutions for homotopy continuation algorithm, and EM procedure picks one among the n!n! solutions [62]. Fixing the precision of the result, the complexity of the homotopy continuation algorithm is proportional to the number of solutions, providing a lower bound Ω(n!)\Omega(n!).

For the SDP relaxation algorithm 3, the interior method costs O(NM)O(N\sqrt{M}), where NN is the cost of solving Newton system in iterations, MM characterizes the scale of semi-definite constrains [46]. Denote msm_{s} as the size of the moment matrix. In the context of unlabeled sensing and algorithm 3, N=O(ms6)N=O(m_{s}^{6}), M=O(msn)M=O(m_{s}n), and ms(n+dd)=(2n+1n)m_{s}\leq\binom{n+d}{d}=\binom{2n+1}{n}.

Fixing the precision of results, an upper bound of the complexity of algorithm 3 can be estimated as O((2n+1n)6.5n)O\left(\binom{2n+1}{n}^{6.5}\sqrt{n}\right). Comparing the complexity of homotopy continuation algorithm and algorithm 3:

cost(HC)cost(SDP)=Ω(n!(2n+1n)6.5n)Ω(exp(n)),\frac{cost(HC)}{cost(SDP)}=\Omega\left(\frac{n!}{\binom{2n+1}{n}^{6.5}\sqrt{n}}\right)\subseteq\Omega(\exp(n)), (31)

we know that in the perspective of computational complexity, algorithm 3 has an exponential speed-up over the homotopy continuation algorithm.

For the case n=5n=5, the homotopy continuation algorithm in Bertini took half an hour to find 120 solutions [62], while the SDP algorithm in Sedumi used only 5.585.58 seconds to find the unique solution. In our experiment, we used HOM4PS2 to solve the polynomial system, and it took only 0.790.79 seconds to find 120 solutions. Different implementations seem to influence the actual runtime of the algorithms. For low-rank matrix completion problems, there are also other practical algorithms; future research should test other implementations.

5. Conclusion

In this paper, we give a positive answer to the open problem Open Problem by showing the polynomial system Qn+1Q_{n+1} consisting of the n+1n+1 Newton polynomials has a unique solution. The main Theorem 1 is proved by considering the birationality of the algebraic morphisms defined among several algebraic varieties. Since Qn+1Q_{n+1} has a unique solution, the unlabeled sensing problem can be reduced to the classic rank-one moment matrix completion problem, and the unique solution can be recovered at the lowest relaxation order n+12\lceil\frac{n+1}{2}\rceil.

Enlightened by the main Theorem 1 , we develop symbolic and numerical algorithms 1, 2 and 3. Algorithm 1 and 2 are based on the Groebner basis computation, and algorithm 3 is based on semidefinite programming.

We test and compare the performance of these algorithms for different n,mn,m and SNRSNR. Our experiments show that for clean input without noise, algorithm 2 solving the overdetermined system Qn+1Q_{n+1} is faster than algorithm 1, 3 and the existing algorithm AIEM, as the Groebner basis of the ideal generated by polynomials in Qn+1Q_{n+1} has simple linear form for lexicographic monomial order and graded monomial order. For corrupted input data with noise, our experiment shows that algorithm 3 can return an approximate solution within several minutes. The rank sequence of the moment matrices in algorithm 3 stabilizes at 11 within the numeric tolerance when using clean input data. However, we also notice that due to the degree of the polynomial system Qn+1Q_{n+1} is relatively high for the moment relaxation method, algorithm 3 may fail to recover the unique solution, which should be improved in the future.

Acknowledgements

Lihong Zhi is supported by the National Key R&\&D Program of China (2023YFA1009401) and the National Natural Science Foundation of China (12071467). Manolis C. Tsakiris is supported by the National Key R&\&D Program of China (2023YFA1009402). Jingyu Lu would like to thank the help of Tianshi Yu on the experiments. The authors acknowledge the support of the Institut Henri Poincaré (UAR 839 CNRS-Sorbonne Université) and LabEx CARMIN (ANR-10-LABX-59-01).

References

  • [1] A. Abid and J. Zou. A stochastic expectation-maximization approach to shuffled linear regression. In 2018 56th Annual Allerton Conference on Communication, Control, and Computing (Allerton), pages 470–477, 2018.
  • [2] A. Abid and J. Zou. A stochastic expectation-maximization approach to shuffled linear regression. In Annual Allerton Conference on Communication, Control, and Computing, pages 470–477, 2018.
  • [3] D. J. Bates, J. D. Hauenstein, A. J. Sommese, and C. W. Wampler. Bertini: Software for numerical algebraic geometry. 2019.
  • [4] J. Berthomieu, C. Eder, and M. Safey El Din. Msolve: A library for solving polynomial systems. In Proceedings of the 2021 on International Symposium on Symbolic and Algebraic Computation, pages 51–58, 2021.
  • [5] S. P. Boyd and L. Vandenberghe. Convex optimization. Cambridge university press, 2004.
  • [6] W. Bruns and U. Vetter. Determinantal Rings, volume 1327 of Lecture Notes in Mathematics. Springer, Berlin, Heidelberg, 1988.
  • [7] B. Buchberger. Ein Algorithmus zum Auffinden der Basiselemente des Restklassenringes nach einem nulldimensionalen Polynomideal. PhD thesis, Innsbruck, 1965.
  • [8] B. Buchberger. An algorithmic criterion for the solvability of an algebraic system of equations. Aequationes Mathematicae, 4:374–383, 1970.
  • [9] B. Buchberger. A criterion for detecting unnecessary reductions in the construction of gröbner-bases. In Symbolic and Algebraic Computation, pages 3–21, Berlin, Heidelberg, 1979. Springer Berlin Heidelberg.
  • [10] E. Candes and B. Recht. Exact matrix completion via convex optimization. Communications of the ACM, 55(6):111–119, 2012.
  • [11] E. J. Candes and Y. Plan. Matrix completion with noise. Proceedings of the IEEE, 98(6):925–936, 2010.
  • [12] E. J. Candes and B. Recht. Exact low-rank matrix completion via convex optimization. In 2008 46th Annual Allerton Conference on Communication, Control, and Computing, pages 806–812, 2008.
  • [13] L. Caniglia, A. Galligo, and J. Heintz. Some new effectivity bounds in computational geometry. In International Conference on Applied Algebra, Algebraic Algorithms, and Error-Correcting Codes, pages 131–151. Springer, 1988.
  • [14] A. Cosse and L. Demanet. Rank-one matrix completion is solved by the sum-of-squares relaxation of order two. In 2015 IEEE 6th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), pages 9–12, 2015.
  • [15] A. Cosse and L. Demanet. Stable rank-one matrix completion is solved by the level 2 lasserre relaxation. Foundations of Computational Mathematics, 21:891–940, 2021.
  • [16] A. P. da Silva, P. Comon, and A. L. de Almeida. A finite algorithm to compute rank-1 tensor approximations. IEEE Signal Processing Letters, 23(7):959–963, 2016.
  • [17] P. David, D. DeMenthon, R. Duraiswami, and H. Samet. Softposit: Simultaneous pose and correspondence determination. International Journal of Computer Vision, 59:259–284, 2004.
  • [18] G. Elhami, A. Scholefield, B. B. Haro, and M. Vetterli. Unlabeled sensing: Reconstruction algorithm and theoretical guarantees. In IEEE International Conference on Acoustics, Speech and Signal Processing, pages 4566–4570, Mar. 2017.
  • [19] J.-C. Faugére. A new efficient algorithm for computing gröbner bases (f4). Journal of Pure and Applied Algebra, 139(1):61–88, 1999.
  • [20] J.-C. Faugère, M. S. El Din, and T. Verron. On the complexity of computing gröbner bases for weighted homogeneous systems. Journal of Symbolic Computation, 76:107–141, 2016.
  • [21] M. Fazel, E. Candes, B. Recht, and P. Parrilo. Compressed sensing and robust recovery of low rank matrices. In 2008 42nd Asilomar Conference on Signals, Systems and Computers, pages 1043–1047, 2008.
  • [22] R. Hartshorne. Algebraic Geometry, volume 52 of Graduate Texts in Mathematics. Springer, New York, NY, 1 edition, 1977.
  • [23] D. Henrion, S. Naldi, and M. S. E. Din. Real root finding for low rank linear matrices. Applicable Algebra in Engineering, Communication and Computing, 31(2):101–133, 2020.
  • [24] X. Huang and A. Madan. Cap3: A dna sequence assembly program. Genome Research, 9(9):868–877, 1999.
  • [25] P. Jain, P. Netrapalli, and S. Sanghavi. Low-rank matrix completion using alternating minimization. In Proceedings of the forty-fifth annual ACM symposium on Theory of computing, pages 665–674, 2013.
  • [26] R. Ji, Y. Liang, L. Xu, and W. Zhang. A concave optimization-based approach for joint multi-target track initialization. IEEE Access, 7:108551–108560, 2019.
  • [27] L. Keller, M. Jafari Siavoshani, C. Fragouli, K. Argyraki, and S. Diggavi. Identity aware sensor networks. In IEEE INFOCOM 2009, pages 2177–2185, 2009.
  • [28] Y. N. Lakshman. A single exponential bound on the complexity of computing gröbner bases of zero dimensional ideals. In Effective methods in algebraic geometry, pages 227–234. Springer, 1991.
  • [29] Y. N. Lakshman and D. Lazard. On the complexity of zero-dimensional algebraic systems. In Effective methods in algebraic geometry, pages 217–225. Springer, 1991.
  • [30] J. Lasserre, M. Laurent, and P. Rostalski. Semidefinite characterization and computation of zero-dimensional real radical ideals. Foundations of Computational Mathematics, 8:607–647, 2008.
  • [31] J. Lasserre, M. Laurent, and P. Rostalski. A prolongation-projection algorithm for computing the finite real variety of an ideal. Theoretical Computer Science, 410(27-29):2685–2700, 2009.
  • [32] J. Lasserre, M. Laurent, and P. Rostalski. A unified approach to computing real and complex zeros of zero-dimensional ideals. In Emerging applications of algebraic geometry, volume 149 of IMA Vol. Math. Appl., pages 125–155. Springer, New York, 2009.
  • [33] J. B. Lasserre. Global optimization with polynomials and the problem of moments. SIAM Journal on optimization, 11(3):796–817, 2001.
  • [34] T.-L. Lee, T.-Y. Li, and C.-H. Tsai. Hom4ps-2.0: a software package for solving polynomial systems by the polyhedral homotopy continuation method. Computing, 83(2-3):109–133, 2008.
  • [35] F. Li, K. Fujiwara, F. Okura, and Y. Matsushita. Shuffled linear regression with outliers in both covariates and responses. International Journal of Computer Vision, 131:732–751, 2023.
  • [36] T.-Y. Li. Numerical solution of multivariate polynomial systems by homotopy continuation methods. Acta numerica, 6:399–436, 1997.
  • [37] H. Liang, J. Lu, M. C. Tsakiris, and L. Zhi. A field-theoretic approach to unlabeled sensing, 2023.
  • [38] R. Ma, T. T. Cai, and H. Li. Optimal estimation of bacterial growth rates based on a permuted monotone matrix. Biometrika, 108(3):693–708, 10 2020.
  • [39] Y. Ma and L. Zhi. Computing real solutions of polynomial systems via low-rank moment matrix completion. In Proceedings of the 37th international symposium on Symbolic and algebraic computation, pages 249–256, 2012.
  • [40] F. Macaulay. The Algebraic Theory of Modular Systems, volume 19 of Cambridge Tracts in Mathematics and Mathematical Physics. Cambridge, 1916.
  • [41] M. Marques, M. Stošić, and J. Costeira. Subspace matching: Unique solution to point matching with geometric constraints. In 2009 IEEE 12th International Conference on Computer Vision, pages 1288–1294, 2009.
  • [42] H. Melánová, B. Sturmfels, and R. Winter. Recovery from power sums. Experimental Mathematics, pages 1–10, 2022.
  • [43] A. Morgan and A. Sommese. Computing all solutions to polynomial systems using homotopy continuation. Applied Mathematics and Computation, 24(2):115–138, 1987.
  • [44] A. Narayanan and V. Shmatikov. Robust de-anonymization of large sparse datasets. In 2008 IEEE Symposium on Security and Privacy (sp 2008), pages 111–125, 2008.
  • [45] A. Nejatbakhsh and E. Varol. Neuron matching in c. elegans with robust approximate linear regression without correspondence. 2021 IEEE Winter Conference on Applications of Computer Vision (WACV), pages 2836–2845, 2021.
  • [46] A. Nemirovski. Interior point polynomial time methods in convex programming. Lecture notes, 42(16):3215–3224, 2004.
  • [47] L. T. Nguyen, J. Kim, and B. Shim. Low-rank matrix completion: A contemporary survey. IEEE Access, 7:94215–94237, 2019.
  • [48] J. Nie and L. Wang. Semidefinite relaxations for best rank-1 tensor approximations. SIAM Journal on Matrix Analysis and Applications, 35(3):1155–1179, 2014.
  • [49] J. Nie, L. Wang, and Z. Zheng. Low rank tensor decompositions and approximations. Journal of the Operations Research Society of China, pages 1–27, 2023.
  • [50] P. A. Parrilo. Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization. California Institute of Technology, 2000.
  • [51] P. A. Parrilo. Semidefinite programming relaxations for semialgebraic problems. Mathematical programming, 96:293–320, 2003.
  • [52] L. Peng and M. C. Tsakiris. Linear regression without correspondences via concave minimization. IEEE Signal Processing Letters, 27:1580–1584, 2020.
  • [53] B. Recht, M. Fazel, and P. A. Parrilo. Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization. SIAM review, 52(3):471–501, 2010.
  • [54] G. Reid and L. Zhi. Solving polynomial systems via symbolic-numeric reduction to geometric involutive form. Journal of Symbolic Computation, 44(3):280–291, 2009.
  • [55] C. Rose and I. S. Mian. Signaling with identical tokens: Upper bounds with energy constraints. In 2014 IEEE International Symposium on Information Theory, pages 1817–1821, 2014.
  • [56] M. Slawski and E. Ben-David. Linear regression with sparsely permuted data. Electronic Journal of Statistics, 13(1):1–36, 2019.
  • [57] X. Song, H. Choi, and Y. Shi. Permuted linear model for header-free communication via symmetric polynomials. In 2018 IEEE International Symposium on Information Theory (ISIT), pages 661–665, 2018.
  • [58] J. F. Sturm. Using sedumi 1.02, a MATLAB toolbox for optimization over symmetric cones. Optimization Methods and Software, 11(1-4):625–653, 1999.
  • [59] D. Sun, K.-C. Toh, Y. Yuan, and X.-Y. Zhao. Sdpnal+: A matlab software for semidefinite programming with bound constraints (version 1.0), 2019.
  • [60] M. C. Tsakiris. Matrix recovery from permutations: An algebraic geometry approach. In 2023 IEEE International Symposium on Information Theory (ISIT), pages 2511–2516, 2023.
  • [61] M. C. Tsakiris and L. Peng. Homomorphic sensing. In K. Chaudhuri and R. Salakhutdinov, editors, Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 6335–6344. PMLR, 09–15 Jun 2019.
  • [62] M. C. Tsakiris, L. Peng, A. Conca, L. Kneip, Y. Shi, and H. Choi. An algebraic-geometric approach for linear regression without correspondences. IEEE Transactions on Information Theory, 66(8):5130–5144, 2020.
  • [63] J. Unnikrishnan, S. Haghighatshoar, and M. Vetterli. Unlabeled sensing: Solving a linear system with unordered measurements. In 2015 53rd Annual Allerton Conference on Communication, Control, and Computing (Allerton), pages 786–793, 2015.
  • [64] J. Unnikrishnan, S. Haghighatshoar, and M. Vetterli. Unlabeled sensing with random linear measurements. IEEE Transactions on Information Theory, 64(5):3237–3253, 2018.
  • [65] B. Vandereycken. Low-rank matrix completion by riemannian optimization. SIAM Journal on Optimization, 23(2):1214–1236, 2013.
  • [66] J. Verschelde. Algorithm 795: Phcpack: A general-purpose solver for polynomial systems by homotopy continuation. ACM Transactions on Mathematical Software (TOMS), 25(2):251–276, 1999.
  • [67] Z. Yang, H. Zhao, and L. Zhi. Verifyrealroots: A matlab package for computing verified real solutions of polynomials systems of equations and inequalities. Journal of Systems Science and Complexity, 36(2):866–883, 2023.
  • [68] Z. Yang, L. Zhi, and Y. Zhu. Verified error bounds for real solutions of positive-dimensional polynomial systems. In Proceedings of the 38th International Symposium on Symbolic and Algebraic Computation, pages 371–378, 2013.