Sum-of-square-of-rational-function based representations of positive semidefinite polynomial matrices

TL;DR

This paper establishes sum-of-square-of-rational-function representations for positive semidefinite polynomial matrices on various sets, providing a numerical method and demonstrating its effectiveness through tests.

Contribution

It introduces a novel method for representing positive semidefinite polynomial matrices using sosrf-based forms on specific sets, with a two-stage diagonalization and low-rank approximation approach.

Findings

Method successfully computes sosrf-representations for various polynomial matrices.

Numerical tests validate the approach's accuracy and efficiency.

Applicable to matrices positive on real line, intervals, and strips.

Abstract

The paper proves sum-of-square-of-rational-function based representations (shortly, sosrf-based representations) of polynomial matrices that are positive semidefinite on some special sets: $\mathbb{R}^n;$ $\mathbb{R}$ and its intervals $[a,b]$, $[0,\infty)$; and the strips $[a,b] \times \mathbb{R} \subset \mathbb{R}^2.$ A method for numerically computing such representations is also presented. The methodology is divided into two stages: (S1) diagonalizing the initial polynomial matrix based on the Schm\"{u}dgen's procedure \cite{Schmudgen09}; (S2) for each diagonal element of the resulting matrix, find its low rank sosrf-representation satisfying the Artin's theorem solving the Hilbert's 17th problem. Some numerical tests and illustrations with \textsf{OCTAVE} are also presented for each type of polynomial matrices.

Equations

$$\mathcal{D}=\{x\in {\mathbb{R}}^n|\mathbf{F}(x)⪰0\}=\{x\in {\mathbb{R}}^n|u^{\mathrm{T}}\mathbf{F}(x)u\ge 0,\forall u\in {\mathbb{R}}^n\}.$$

$$\mathbf{F}=\sum _{i=1}^r{\mathbf{A}}_i^{\mathrm{T}}{\mathbf{A}}_i,$$

$$\mathbf{F}=\sum _{i=1}^r\left(\frac{1}{b}{\mathbf{A}}_i^{\mathrm{T}}\right)\left(\frac{1}{b}{\mathbf{A}}_i\right).$$

$$f^M(x,y)=1+x^2{y}^4+x^4{y}^2-3x^2{y}^2.$$

$$p^2{f}^M(x,y)=(x^2-y^2{)}^2+[xy(p-2){]}^2+[x^{2}y(p-2){]}^2+[x{y}^2(p-2){]}^2$$

$$f(x)=p^2(x)+x{q}^2(x),$$

$$f\ge 0\text{mbox}on{\mathbb{R}}^n$$

$$⟺\exists b\in \mathbb{R}[x],G\in {\mathbb{S}}^e\mathbb{R}:b(x{)}^{2}f(x)=\pi (x{)}^{\mathrm{T}}G\pi (x),$$

$$\text{mbox}minimize\left\{\text{mbox}\mathrm{rank}(X)|X\in {\mathbb{R}}^{m\times p},\varphi (X)=u\right\},$$

$$X_{\pm }=\left[\begin{array}{cc}\alpha & 0\\ \pm {\beta }^{\mathrm{T}}& \alpha I\end{array}\right],\tilde{F}=\left[\begin{array}{cc}{\alpha }^3& 0\\ 0& \alpha (\alpha C-{\beta }^{\mathrm{T}}\beta )\end{array}\right]\in R^{m\times m},$$

$$F=\left[\begin{array}{cc}\alpha & \beta \\ {\beta }^{\mathrm{T}}& C\end{array}\right].$$

$$X_{+}{X}_{-}=X_{-}{X}_{+}={\alpha }^{2}I,{\alpha }^{4}F=X_{+}\tilde{F}{X}_{+}^{\mathrm{T}},\tilde{F}=X_{-}F{X}_{-}^{\mathrm{T}}.$$

$${\mathbf{X}}_{+}{\mathbf{X}}_{-}={\mathbf{X}}_{-}{\mathbf{X}}_{+}=b{I}_m,b^2\mathbf{F}={\mathbf{X}}_{+}\text{mbox}\mathrm{diag}(d_1,\dots ,d_m){\mathbf{X}}_{+}^{\mathrm{T}}$$

$$\mathbf{F}=\left[\begin{array}[]{cc}\alpha&\beta\\ \beta^{\mathrm{T}}&\mathbf{C}\end{array}\right],\text{ where }\mathbf{C}\in\mathbb{S}^{m-1}\mathbb{R}[t],\quad\beta\in\mathbb{R}[t]^{1\times(m-1)}.$$

$$\alpha_{0}=f_{11}^{(0)},\quad\mathbf{X}_{0\pm}=\left[\begin{array}[]{cc}\alpha_{0}&0\\ \pm\beta_{0}^{\mathrm{T}}&\alpha_{0}I\end{array}\right],\quad\widetilde{\mathbf{F}}_{0}=\left[\begin{array}[]{cc}\alpha_{0}^{3}&0\\ 0&\mathbf{B}_{0}\end{array}\right],\enskip\mathbf{D}^{(0)}=\widetilde{\mathbf{F}}_{0},$$

$$\alpha_{i-1},\enskip\mathbf{X}_{\pm}^{(i-1)},\enskip\mathbf{D}^{(i-1)}=\left[\begin{array}[]{cc}\mbox{\rm diag}(d_{0}^{(i-1)},\ldots,d_{i-1}^{(i-1)})&0\\ 0&\mathbf{B}_{i-1}\end{array}\right],$$

$$\alpha_{i}=f_{11}^{(i)},\enskip\mathbf{X}_{i\pm}=\left[\begin{array}[]{cc}\alpha_{i}&0\\ \pm\beta_{i}^{\mathrm{T}}&\alpha_{i}I\end{array}\right],\enskip\widetilde{\mathbf{F}}_{i}=\left[\begin{array}[]{cc}\alpha_{i}^{3}&0\\ 0&\mathbf{B}_{i}\end{array}\right]\in\mathbb{S}^{m-i}\mathbb{R}[t]$$

$${\mathbf{X}}_{i+}{\mathbf{X}}_{i-}={\mathbf{X}}_{i-}{\mathbf{X}}_{i+}={\alpha }_i^{2}I,{\alpha }_i^4{\mathbf{F}}_i={\mathbf{X}}_{i+}{\tilde{\mathbf{F}}}_i{\mathbf{X}}_{i+}^{\mathrm{T}},{\tilde{\mathbf{F}}}_i={\mathbf{X}}_{i-}{\mathbf{F}}_i{\mathbf{X}}_{i-}^{\mathrm{T}}.$$

$$b^{(i)},\quad\mathbf{X}_{\pm}^{(i)},\quad\mathbf{D}^{(i)}=\left[\begin{array}[]{cc}\mbox{\rm diag}(d_{0}^{(i)},\ldots,d_{i}^{(i)})&0\\ 0&\mathbf{B}_{i}\end{array}\right]\in\mathbb{S}^{m}\mathbb{R}[t],$$

$$b^{(0)}={\alpha }_0^2,{\mathbf{X}}_{\pm }^{(0)}={\mathbf{X}}_{0\pm }\text{ and }{\mathbf{D}}^{(0)}={\tilde{\mathbf{F}}}_0\in {\mathbb{S}}^m\mathbb{R}[t],d_0^{(0)}={\alpha }_0^3.$$

$${\mathbf{X}}_{+}^{(0)}{\mathbf{X}}_{+}^{(0)}={\mathbf{X}}_{-}^{(0)}{\mathbf{X}}_{+}^{(0)}={\alpha }_0^2{I}_m,{\alpha }_0^4\mathbf{F}={\mathbf{X}}_{+}^{(0)}{\mathbf{D}}^{(0)}[{\mathbf{X}}_{+}^{(0)}{]}^{\mathrm{T}},{\mathbf{D}}^{(0)}={\mathbf{X}}_{-}^{(0)}{\mathbf{F}}_0[{\mathbf{X}}_{-}^{(0)}{]}^{\mathrm{T}}.$$

$${\mathbf{X}}_{1+}{\mathbf{X}}_{1-}={\mathbf{X}}_{1-}{\mathbf{X}}_{1+}={\alpha }_1^2,{\alpha }_1^4{\mathbf{F}}_1={\mathbf{X}}_{1+}{\tilde{\mathbf{F}}}_1{\mathbf{X}}_{1+}^{\mathrm{T}},{\tilde{\mathbf{F}}}_1={\mathbf{X}}_{1-}{\mathbf{F}}_1{\mathbf{X}}_{1-}^{\mathrm{T}}.$$

$${\alpha }_1^4({\alpha }_0^4\mathbf{F})$$

$$\displaystyle=\mathbf{X}_{0+}\left[\begin{array}[]{cc}\alpha_{1}^{4}\alpha_{0}^{3}&0\\ 0&\mathbf{X}_{1+}\widetilde{\mathbf{F}}_{1}\mathbf{X}_{1+}^{\mathrm{T}}\end{array}\right]\mathbf{X}_{0+}^{\mathrm{T}}$$

$$\displaystyle=\mathbf{X}_{0+}\left[\begin{array}[]{cc}\alpha_{1}&0\\ 0&\mathbf{X}_{1+}\end{array}\right]\left[\begin{array}[]{cc}\alpha_{1}^{2}\alpha_{0}^{3}&0\\ 0&\widetilde{\mathbf{F}}_{1}\end{array}\right]\left[\begin{array}[]{cc}\alpha_{1}&0\\ 0&\mathbf{X}_{1+}^{\mathrm{T}}\end{array}\right]\mathbf{X}_{0+}^{\mathrm{T}}$$

$$\displaystyle\left[\begin{array}[]{cc}\alpha_{1}^{2}\alpha_{0}^{3}&0\\ 0&\widetilde{\mathbf{F}}_{1}\end{array}\right]$$

$$\displaystyle=\left[\begin{array}[]{cc}\alpha_{1}&0\\ 0&\mathbf{X}_{1-}\end{array}\right]\left[\begin{array}[]{cc}\alpha_{0}^{3}&0\\ 0&\mathbf{F}_{1}\end{array}\right]\left[\begin{array}[]{cc}\alpha_{1}&0\\ 0&\mathbf{X}_{1-}^{\mathrm{T}}\end{array}\right]$$

$$\displaystyle=\left[\begin{array}[]{cc}\alpha_{1}&0\\ 0&\mathbf{X}_{1-}\end{array}\right]\mathbf{X}_{0-}\mathbf{F}_{0}\mathbf{X}_{0-}^{\mathrm{T}}\left[\begin{array}[]{cc}\alpha_{1}&0\\ 0&\mathbf{X}_{1-}^{\mathrm{T}}\end{array}\right].$$

$${\mathbf{X}}_{+}^{(1)}$$

$${\mathbf{X}}_{-}^{(1)}$$

$${\mathbf{D}}^{(1)}$$

$${\mathbf{X}}_{+}^{(1)}{\mathbf{X}}_{-}^{(1)}={\mathbf{X}}_{-}^{(1)}{\mathbf{X}}_{+}^{(1)}=b^{(1)}{I}_m,[b^{(1)}{]}^2\mathbf{F}={\mathbf{X}}_{+}^{(1)}{\mathbf{D}}^{(1)}[{\mathbf{X}}^{(1)}{]}_{+}^{\mathrm{T}},{\mathbf{D}}^{(1)}={\mathbf{X}}_{-}^{(1)}\mathbf{F}[{\mathbf{X}}_{-}^{(1)}{]}^{\mathrm{T}}.$$

$$b^{(i)}$$

$${\mathbf{X}}_{+}^{(i)}$$

Full text

Keywords: nonnegative polynomial, positive semidefinite polynomial matrix, sum of Hermitian squares of polynomial matrices, Schmüdgen’s diagonalization

1 Introduction and preliminary

Nonnegative polynomials, i.e., polynomials with scalar-coefficients, appears in a variety of mathematical problems and applications. The problem of approximating such polynomials as sums of squares of polynomials have been studied more deeply in both theoretical and practical points of view, see eg., [13] and references therein. In our opinion, this idea can be generalized to matrix polynomials, i.e., polynomials with scalar-matrix coefficients. A matrix polynomial can also be written as a polynomial matrix, i.e., a matrix with polynomial entries. However, the terminology “nonnegative” for scalar polynomials must be generalized as “positive semidefinite”, which is usually used for (scalar or polynomial) matrices. Because of this reason, we prefer “polynomial matrices” to “matrix polynomial” in this paper.

There have been a variety of researches focusing on the generalization of Positivstellensätze to psd-polynomial matrices, i.e., $m\geq 2,$ and its applications (see in, eg., [8] and references therein). However, from the practical point of view, most applications lead to problems over univariate/bivariate polynomial matrices. This motivates us to focus on these two types of polynomial matrices.

Let $\mathbb{R}[x]$ be the ring of $n$-real variable polynomials with real coefficients. Let $\mathbb{R}^{m\times m}$ and $\mathbb{S}^{m}\mathbb{R}$ denote the set of all $m\times m$ matrices with real entries and its subset of symmetric matrices, respectively. By $.^{\mathrm{T}}$ we denote the transpose of matrices. For $A\in\mathbb{S}^{m}\mathbb{R},$ by $A\succeq 0$ we mean $A$ is positive semidefinite, i.e., $u^{\mathrm{T}}Au\geq 0$ for all $u\in\mathbb{R}^{n}.$ We denote $\mathbb{S}_{+}^{m}\mathbb{R}$ the set of all positive semidefinite matrices in $\mathbb{S}^{m}\mathbb{R}.$ Moreover, for two matrices $A$ and $B$ in $\mathbb{S}^{m}\mathbb{R},$ we write $A\succeq B$ if $A-B\in\mathbb{S}_{+}^{m}\mathbb{R}.$ For a subset $\mathcal{D}\subseteq\mathbb{R}^{n},$ and a polynomial matrix $\mathbf{F}\in\mathbb{S}^{m}\mathbb{R}[x],$ we say $\mathbf{F}$ to be positive semidefinite on $\mathcal{D}$ if

[TABLE]

We will say “psd” instead of “positive semidefinite”. Someone calls a psd polynomial matrix on $\mathcal{D}=\mathbb{R}^{n}$ the global psd one. Throughout this paper, we fix

• •

$m,n$ as above, i.e., $m$ is the size of a polynomial matrix and $n$ is the number of variables in the matrix;

• •

$d:$ the maximum of degrees of $m^{2}$ polynomials in the matrix.

A symmetric polynomial matrix $\mathbf{F}$ is called a “sum of squares of polynomial matrices” (resp., a *sum of squares of rational functions matrices111 A rational function is a ratio of polynomials. *) if it is a finitely many sum of the form

[TABLE]

where the matrices $\mathbf{A}_{i}$ has polynomial entries (resp., rational-function entries) defined on the whole space $\mathbb{R}^{n}.$ We shortly say a sum of squares matrix a sos-matrix while the other a sosrf-matrix. It is clear that if $b^{2}\mathbf{F},$ $b\in\mathbb{R}[x],$ is sos then $\mathbf{F}$ is sosrf since

[TABLE]

It can be seen that any sos-matrix or sosrf-matrix is psd on $\mathbb{R}^{n}.$ The converse direction is always true for sosrf-matrices since Proposition 4, but not for sos-matrices. The observation of this issue is easy to see for scalar polynomials including a famous counter example of the Motzkin polynomial

[TABLE]

It cannot be expressed as a sum of squares of polynomials in $\mathbb{R}[x,y]$ but

[TABLE]

with $p(x,y)=x^{2}+y^{2}.$

When $m=1,$ polynomial matrix concept here coincides with the scalar polynomial as usual, and “psd-polynomial matrices” now means “nonnegative polynomials” as well. There is a number of well-known Positivstellensätze for nonnegative/positive scalar polynomials [7, 22], leading to representations of scalar polynomials nonnegative on several subsets of $\mathbb{R}$ using sos-representations of other polynomials. For example, it is well known in literature that, see also in eg., [2], any univariate real polynomial nonnegative on $\mathbb{R}$ can be expressed as a sum of two squares of real polynomials. This is called a sos-representation of the initial polynomial. And, a polynomial $f$ nonnegative on $[0,\infty)$ can be written as

[TABLE]

for some real polynomials $p,q.$ We will call the later representation the “sos-based representation” of $f.$

In [9], the authors propose an algorithm for decomposing a nonnegative polynomial on $\mathbb{R}^{n}$ as a sum of squares of rational functions relying on the idea of Reznick [18], saying that the common denominator is a power of the sum of squares of coordinate functions. In this paper, we deal with the problem of finding a sos-based and/or sosrf-based representation of a polynomial matrix which is positive semidefinite over one of the sets: $\mathbb{R}^{n},$ $[a,b],$ $[0,+\infty)$ and strips in $\mathbb{R}^{2}.$ The idea is to combine Schmüdgen’s procedure of diagonalizing polynomial matrices and Levenberg-Marquardt algorithm [15] finding sosrf-representations of the diagonal polynomials of the resulting diagonal matrices. Our method for finding sosrf-based representations for psd-polynomial matrices is hence divided into two stages:

(S1) diagonalizing the initial polynomial matrix $\mathbf{F}$ suggested by Proposition 1 below so that (3) holds true. Note that the initial polynomial matrix is positive semidefinite if and only if the resulting diagonal polynomials are all nonnegative;

(S2) finding low-rank sos-representations of the resulting diagonal polynomials by applying the algorithm proposed in [9]. The representation for the initial polynomial matrix we need is then obtain by substituting to the relation in the first stage.

As mentioned earlier, the second stage is to find sosrf-representations of scalar polynomials. It is well-known that any sum of squares polynomial can be determined by its Gram matrix. With the help of Artin’s theorem answering to the Hilbert’s 17th problem, one notes

[TABLE]

where $e={n+\deg(b^{2}f)\choose n},$ $\pi(x)$ is the vector of monomials $x^{\alpha}:=x_{1}^{\alpha_{1}}\ldots x_{n}^{\alpha_{n}}$ with degree not exceed $\deg(b^{2}f),$ and $G$ is real symmetric and positive semidefinite, which is called a Gram matrix of $b^{2}f.$ By comparing the coefficients, we see that each coefficient on the left side is linearly dependent on the entries of $G.$

In addition, a Gram matrix of a sos-polynomial can be found to be low rank. This suggests us to propose and solve the matrix rank minimization problem

[TABLE]

where $\phi:\mathbb{R}^{m\times p}\rightarrow\mathbb{R}^{l}$ is a differentiable map, that can be applied in our second stage. The rank function is non-convex, even the feasibility region is smooth. The idea of this stage is to check whether the problem (1) has a numerical solution of rank $r,$ step by step, for $r=1,2,\ldots$ At each step, with respect to each $r,$ the Levenberg-Marquardt method (for short, LM-method) [15] solving least square problems is applied to the function defined by $\phi(X)=u.$ The differentiability of $\phi$ guarantees the existence of the Jacobian matrix of $\phi$ in Levenberg-Marquardt iterations.

The problem of finding a sos-representation of a scalar sos-polynomial is thus equivalent to which of solving a linear system finding a polynomial’s Gram matrix with low rank. This is reasonable to apply the model (1) with respect to a suitable linear map $\phi.$ Furthermore, any polynomial’s Gram matrix must be symmetric, and the existence of a Gram matrix is equivalent to which of its Cholesky factor222 A positive semidefinite matrix $X$ can always be factorized as $X=YY^{\mathrm{T}},$ where $Y$ is called a Cholesky factor of $X.$ , where both are the same rank. In this sense, $\phi$ will be a quadratic function in the Cholesky factor. Our second stage, as seen in Section 3, deals properly with finding such Cholesky factors. The polynomial’s Gram matrix is then derived, and so is the polynomial. This is done by applying the LM-method.

The LM-method has been generalized for complex setting in [20, 21], which is called the generalized LM-method, or gLM-method. It solves the least squares problems with real-valued function in complex variables. The corresponding Matlab toolbox is called COT [21]. The authors there convert the complex setting to the real one by using 1-1 linear transformations that map any complex $n$-tuples $z\in\mathbb{C}^{n},$ to $(\mbox{\rm Re}(z),\mbox{\rm Im}(z))\in\mathbb{R}^{2n},$ and/or to $(z,\bar{z})$ with the help of Wirtinger calculus and the complex Taylor series expansions. Even the gLM-method aims at solving least square problems with real-valued functions in complex variables, it is also valid for real-valued functions in real variables. In our work we recode the gLM-method in OCTAVE that makes more convenience for the readers. Unlike as “lsqnonlin” function in Matlab, OCTAVE codes, as those of COT as well, additionally need the Jacobian/gradient of the function as an input. We hence give detailed formulas for Jacobian matrices with respect to particular cases that gLM-method needs. The numerical tests in this paper are all implemented in GNU OCTAVE 4.4.2333Available at https://octave.sourceforge.io/ and Ubuntu 14.04, on a desktop Intel(R) Core(TM) i3-3220 CPU@3.30GHz and 4GB RAM.

This paper is organized as follows. Section 2 gives a clear explaination to Schmüdgen diagonalization, which is considered as the first stage of our method, and presents the corresponding algorithm. Some numerical tests are also illustrated. The second stage is based on a rank matrix minimization model devoted in Section 3. After proposing an algorithm solving such a model, we apply to solve our second stage. We give some numerical examples as well. In section 4 we prove sosrf-based representations for polynomial matrices that are psd on the sets $\mathbb{R}^{n},$ ($n\geq 1$), some intervals in $\mathbb{R},$ and some strips in $\mathbb{R}^{2}.$ The conclusion of the paper is devoted in the last section.

2 Diagonalization of polynomial matrices

2.1 Schmüdgen’s procedure

We start this section by recalling a very important matrix decomposition by Schmüdgen [19] which is a key part for numerically computing sum of Hermitian squares of polynomial matrices below. For a symmetric matrix $F$ with entries in a commutative unitary ring $R,$ set

[TABLE]

where $\alpha,$ $\beta$ and $C$ are parts of $F$ in the following partition

[TABLE]

Then the following relations hold

[TABLE]

Schmüdgen’s procedure (2) is applied to diagonalize a polynomial matrix as in the following proposition, viewed as an extension of Artin’s theorem solving Hilbert’s 17th problem.

Proposition 1**.**

[19, Corollary 10]* Let $\mathbf{F}\in\mathbb{S}^{m}\mathbb{R}[x].$ Then there exist polynomial matrices $\mathbf{X}_{+},\mathbf{X}_{-}\in\mathbb{R}[x]^{m\times m}$ and polynomials $b,d_{j}\in\mathbb{R}[x],$ $j=1,\ldots,m,$ such that*

[TABLE]

and $X_{-}\mathbf{F}X_{-}^{\mathrm{T}}=\mbox{\rm diag}(d_{1},\ldots,d_{m}).$

2.2 Algorithmetic implementation

Even though Schmüdgen’s proof (of Proposition 1) suggests us an algorithm to find $b,\mathbf{X}_{\pm}$ and $\mathbf{D},$ it was not numerically clarified in his seminal paper. In this section, we explain in more detailed how Schmüdgen’s procedure (2) is applied to obtain $b,\mathbf{X}_{\pm}$ and $\mathbf{D}$ as in Proposition 1, step by step such that one can implement the algorithm in any scientific programming language.

Given a polynomial matrix $\mathbf{F}\in\mathbb{S}^{m}\mathbb{R}[t]$ that is partitioned as

[TABLE]

With the help of the procedure in (2), our algorithm below starts with $\mathbf{F}=[f_{ij}^{(0)}]\in\mathbb{S}^{m}\mathbb{R}[t]$ and

[TABLE]

where $\mathbf{B}_{0}=\alpha_{0}^{2}\mathbf{C}_{0}-\alpha_{0}\beta_{0}^{\mathrm{T}}\beta_{0}\in\mathbb{S}^{m-1}\mathbb{R}[t].$ The next iteration will apply procedure (2) to the matrix $\mathbf{B}_{0}.$

We would emphasize that, afterwards, the matrices with subscripts $\mathbf{X}_{i+},\mathbf{X}_{i-},$ $\mathbf{F}_{i}\in\mathbb{S}^{m-i}\mathbb{R}[t]$ have sizes that are changed through the different iterations; while the others, indexed by superscripts such as $b^{(i)},\mathbf{X}_{\pm}^{(i)},$ $\mathbf{D}^{(i)},$ will be updated from $\mathbf{X}_{i+},\mathbf{X}_{i-},\mathbf{F}_{i}$ as in Proposition 2. The matrices with superscripts must have the same size as $\mathbf{F}$ at every iteration. The algorithm proceeds by induction (on $m$) until $\mathbf{B}_{i}\in\mathbb{S}^{1}\mathbb{R}[t]$ or $\mathbf{B}_{i}$ is identically zero, i.e., $\mathbf{D}^{(i)}$ is diagonal and the outputs are $\mathbf{D}=\mathbf{D}^{(i)},$ $\mathbf{X}_{\pm}=\mathbf{X}_{\pm}^{(i)},b=b^{(i)}.$

From the computational point of view, the matrix $\mathbf{D}$ should be treated as its diagonal vector. We thus only need to compute its main diagonal element by element. In other words, at Iteration $i,$ for instant, we only need to compute $d_{i}$ (exactly equals to $\alpha_{i}^{3}$) and update $d_{0},d_{1},\ldots,d_{i-1},$ but do not pay attention for $d_{r},i<r\leq m.$

Below, we will see in our algorithm that the data at iteration $i$ ($i\geq 2$) is updated from that at iteration $i-1.$ More precisely, suppose we are given

[TABLE]

where $\mathbf{B}_{i-1}=\alpha_{i-1}^{2}\mathbf{C}_{i-1}-\alpha_{i-1}\beta_{i-1}^{\mathrm{T}}\beta_{i-1}\in\mathbb{S}^{m-i+1}\mathbb{R}[t].$ Then by applying procedure (2) to $\mathbf{F}_{i}=\widetilde{\mathbf{B}}_{i-1}=[f_{rs}^{(i)}]$ we compute

[TABLE]

that satisfy the relation (3):

[TABLE]

One then updates

[TABLE]

where $d_{k}^{(i)}$’s are determined as in Proposition 2 below.

To unify the notations and to make the readers more convenient, we implement two first iterations as follows.

Iteration 0. Define $\alpha_{0},\mathbf{X}_{0\pm},$ $\mathbf{B}_{0}\in\mathbb{S}^{m-1}\mathbb{R}[t]$ and $\widetilde{\mathbf{F}}_{0}\in\mathbb{S}^{m}\mathbb{R}[t]$ from $\mathbf{F}$ as in (5), then put

[TABLE]

It is clear that (3) holds:

[TABLE]

Iteration 1. Set $\mathbf{F}_{1}=\mathbf{B}_{0}$ and determine $\alpha_{1},\mathbf{X}_{1\pm},$ $\widetilde{\mathbf{F}}_{1}=\begin{bmatrix}\alpha_{1}^{3}&0\\ 0&\mathbf{B}_{1}\end{bmatrix}\in\mathbb{S}^{m-1}\mathbb{R}[t]$ by applying procedure (2) to $\mathbf{F}_{1}.$ These satisfy relation (3):

[TABLE]

This, combining with (10), implies that

[TABLE]

and

[TABLE]

We hence obtain that $b^{(1)}=\alpha_{1}^{2}\alpha_{0}^{2}$ and

[TABLE]

Provided by (20) and (37) , $\mathbf{X}_{\pm}^{(1)}$ and $\mathbf{D}^{(1)}$ also satisfy the relation (3):

[TABLE]

The diagonal of $\mathbf{D}^{(1)},$ as mentioned earlier, should be updated only $d_{1}^{(1)}=d_{1}^{(0)}\alpha_{1}^{2}$ such that $\mathbf{D}^{(1)}$ has the form (6).

The two early iterations show us the numerical formulas to compute $b,\mathbf{X}_{\pm}$ and $\mathbf{D}$ in the next iterations such that they satisfy (3). The mathematical explanation in more detail, compared with (6), (7) and (9), is showed in the following proposition.

Proposition 2**.**

With notation as in (6), (7) and (9). Suppose the data at Iteration $i-1$ with $i\geq 1$ is already given: $b^{(i-1)},\mathbf{X}_{\pm}^{(i-1)}$ and $\mathbf{D}^{(i-1)}.$ Then for all $i=2,\ldots,m-1,$ we have the update rules

[TABLE]

These matrices furthermore satisfy the relation described in (3):

[TABLE]

Proof.

Since the section of Iteration 1, the proposition is true with $i=1.$ Now suppose that the proposition is true for some $i\geq 1,$ namely,

[TABLE]

where $b^{(i)}=b^{(i-1)}\alpha_{i}^{2}=\alpha_{1}^{2}\alpha_{2}^{2}\ldots\alpha_{i}^{2},$

[TABLE]

and $\mathbf{D}^{(i)}$ is defined as in (2). That is

[TABLE]

We will show that (46) also holds for $i+1.$ By inductive hypothesis, with noting that ${b^{(i+1)}}=\alpha_{i+1}^{2}{b^{(i)}},$ one has

[TABLE]

This yields the form of $\mathbf{D}^{(i+1)}$ and $\mathbf{X}_{+}^{(i+1)}.$

For the last relation, $\mathbf{X}^{(i+1)}_{-}\mathbf{F}{\mathbf{X}^{(i+1)}_{-}}^{\mathrm{T}}=\mathbf{D}^{(i+1)},$ we write

[TABLE]

But $\mathbf{D}^{(i)}=\mathbf{X}_{-}^{(i)}\mathbf{F}\mathbf{X}_{-}^{(i)}$ provided by (46), we can conclude

[TABLE]

Thus

[TABLE]

and we are done. ∎

Remark 1**.**

We would note that $\mathbf{D}^{(i)}$ in Proposition 2 can also expressed via $\mathbf{D}^{(i-1)}$ as

[TABLE]

if $\mathbf{F}_{i}$ is invertible.

Proposition 2 gives us formulas updating the outputs $b,\mathbf{X}_{\pm}$ and $\mathbf{D}$ in each step of the algorithm. The algorithm supposes that $\alpha_{i}\neq 0$ at all iterations $i.$ Otherwise, according to [19], one can find an orthogonal matrix $\mathbf{T}$ such that the $(1,1)$-th entry of $\mathbf{T}\mathbf{F}_{i}\mathbf{T}^{\mathrm{T}}$ is nonzero. A more clear formula can also be found in [4, Proof of Proposition 5].

Algorithm 1. Dialgonalizing a polynomial matrix.

Input:    $0\neq\mathbf{A}\in\mathbb{S}^{m}\mathbb{R}[t],m\geq 2$

Output: polynomials $b,d_{j}\in\mathbb{R}[t],j=\overline{1,m},$ and matrices $\mathbf{X}_{+},\mathbf{X}_{-}\in\mathbb{R}[t]^{m\times m}$ $\bullet$ Assign: $\alpha_{0},$ $\mathbf{X}_{\pm}=\mathbf{X}_{0,\pm},$ $\widetilde{\mathbf{F}}_{0}=\widetilde{\mathbf{D}}^{(0)}$ and $\mathbf{B}_{0}$ as in (5).

$\bullet$ At Iteration $i\geq 1:$

1. Step 1:

Set $\mathbf{F}_{i}=\mathbf{B}_{i-1}.$ 2. Step 2:

Determine $\alpha_{i},\mathbf{X}_{i-},\mathbf{X}_{i+},\widetilde{\mathbf{F}}_{i}$ from $\mathbf{F}_{i}$ as in (7) which satisfy (8). 3. Step 3:

Update $b^{(i)},\mathbf{X}^{(i)}_{\pm},\mathbf{D}^{(i)}$ as in (2). 4. Step 4:

If $\mathbf{B}_{i}\in\mathbb{S}^{1}\mathbb{R}[t]$ or $\mathbf{B}_{i}$ is zero then stop.

Else set $i=i+1$ then go to Step 1.

In the case Algorithm 1 stops at Iteration $i$ with $0\neq\mathbf{B}_{i}\in\mathbb{S}^{1}\mathbb{R}[t],$ it defines the $(i+1)$-th diagonal element in $\mathbf{D}^{(i+1)}$ as $\mathbf{B}_{i}.$ In this case we set this element as $\alpha_{i+1}.$ It is clear that $\mathbf{F}$ is positive semidefinite on $\mathbb{R}^{n}$ if and only if the diagonal elements $d_{j}$’s are nonnegative on $\mathbb{R}^{n}.$ From Algorithm 1, $d_{j}$’s can be replaced by more simple polynomials as follows.

Corollary 3**.**

The polynomial matrix $\mathbf{F}$ in Proposition 1 is positive semidefinite on $\mathbb{R}^{n}$ if and only if $\alpha_{0},\alpha_{1},\ldots,\alpha_{i+1},$ $i=0,1,\ldots,m-1,$ are nonnegative on $\mathbb{R}^{n}.$ In other words,

[TABLE]

Another corollary of Algorithm 1 is as follows. Cimprič [4, Proposition 5], see also in [8], proves that for a set $\mathfrak{F}$ of polynomial matrices in $\mathbb{S}^{m}\mathbb{R}[x],$ one can find a set of scalar polynomials $\mathfrak{F}^{\prime}\subset\mathbb{R}[x]$ such that

[TABLE]

Even though his proof is suggested by Schmüdgen’s procedure, it is only inductive on matrices of smaller size than $m,$ so that the resulting set of scalar polynomials is not much clear. With the help of Corollary 3, such a set of scalar polynomials can be verified as

[TABLE]

2.3 Numerical illustrations

In OCTAVE, a polynomial will be stored as a vector of its coefficients. For example, the polynomial $f(t)=t^{2}$ will be saved as the vector $f=[0\enskip 0\enskip 1].$ The syntax “conv” of vectors was used to compute the multiplication of polynomials. A polynomial matrix will be saved as a cell, and the multiplication of matrices will be done as the cell multiplication.

2.3.1 Univariate polynomials

Table 1 shows some numerical tests for Schmüdgen’s decomposition of univariate polynomial matrices.

We give here in detailed explaination for the case $(m,d)=(3,2).$ That is, $\mathbf{F}(x)\in\mathbb{S}^{3}\mathbb{R}[t]$ with entries are univariate polynomials of degree 2, and their coefficient vectors are randomly chosen in $(0,1)^{4}.$ OCTAVE produces the $(i,j)$-entries of $\mathbf{F}$ in form $\mathbf{F}[i,j]$ as follows

[TABLE]

Algorithm 1 gives us the polynomial $b$ of degree 16, $\mathbf{X}_{p}$’s and $\mathbf{X}_{m}$’s entries are of degree 8, and diagonal elements’ degrees are all 18.

We next give a particular example that is also returned later.

Example 1**.**

We start with an example of $2\times 2$ univariate-polynomial matrix

[TABLE]

Applying the algorithm we then obtain

[TABLE]

In addition, it is not hard to see $\mathbf{F}$ is positive semidefinite on $\mathbb{R}$ and

[TABLE]

which illustrates (47).

2.3.2 Multivariate polynomials

To explain the tests in this situation, we need some more notations. Fix three positive numbers $m,n$ and $d.$ Let

[TABLE]

The cardinalities of $\Omega$ and $\Gamma$ are well-known and are determined by $|\Omega|=\binom{n+d}{n}$ and $|\Gamma|=\binom{n+2d}{n},$ respectively.

Our experiments in this situation admits the lexicographic order of monomials:

[TABLE]

For example, for $m=1,$ $n=2,$ $d=6$ then

[TABLE]

The polynomial’s coefficient vector is expressed with respect to the “lex” order as in $\Omega.$ In particular, the Motzkin polynomial mentioned earlier will be written as $f^{M}(x,y)=1-3x^{2}y^{2}+x^{2}y^{4}+x^{4}y^{2},$ and its coefficient vector will be stored in OCTAVE as

[TABLE]

It is worth mentioning that the lengths of multivariate polynomials’ coefficient vectors rapidly increase as $n$ and $d$ are large. This makes the size of the present problem rapidly large as well.

Example 2**.**

Consider the matrix, given in [19], $\mathbf{F}(x,y)=\begin{bmatrix}1+x^{4}y^{2}&xy\\ xy&1+x^{2}y^{4}\end{bmatrix}.$ Algorithm 1 produces $b=(1+x^{4}y^{2})^{2}$ and

[TABLE]

where the coefficient vectors are “full-length”. For example, for $b=(1+x^{4}y^{2})^{2}=1+2x^{4}y^{2}+x^{8}y^{4}$ above has degree 12, and OCTAVE produces its coefficient vector of length $|\Omega_{2,12}|=91$ with most zero entries except for the corresponding coefficients of $x^{0}y^{0},x^{4}y^{2},x^{8}y^{4}$ are $1,2,1,$ respectively.

3 Rank matrix minimization model (1)

3.1 Algorithm

As mentioned earlier, the problem (1) is solved by finding a matrix $X\in\mathbb{R}^{m\times p}$ step by step for $\mbox{\rm rank}(X)=1,2,\ldots,\min\{m,p\}$ such that $\phi(X)=u.$ The search of $X$ believes in the gLM-method [20, 21].

The Levenberg-Marquardt method [15] is a famous method that is used to solve the least squares problems. For convenience to the readers, we summarize this method as follows. A “real data” least square problem minimizes a real function $f(x),$ $x\in\mathbb{R}^{n},$ where

[TABLE]

and $F_{j}$’s are continuously differentiable. According to [15], starting with an initial point $x_{0}\in\mathbb{R}^{n},$ the method finds a sequence of points $\{x_{k}\}\subset\mathbb{R}^{n}$ that converges to a minimizer of $f.$ At each step $k,$ the next point $x_{k+1}$ is determined by applying the first-order approximation of $f$ over an appropriately closed “hyperellipsoid” with center $x_{k}.$ More precisely, one first approximates $f(x_{k}+p)$ by its first order approximation

[TABLE]

Then, one minimizes the approximating function $F(x_{k})+\mbox{\rm{Jac}}F(x_{k})p$ over a hyperellipsoid $E_{k}=\{p=(p_{1},\ldots,p_{n}):\ \sum_{i=1}^{n}(d_{i}^{(k)})^{2}p_{i}^{2}\leq\Delta_{k}\}$ to search a descent direction $p_{*}$ (a minimizer of the approximating function). Either the step bound $\Delta_{k}$ or the scaling parameters $d_{i}^{(k)}$ will be updated after each step $k,$ depending upon the difference between $f(x_{k}+p_{*})$ and $f(x_{k}).$ The next point $x_{k+1}$ will then be decided to be $x_{k}+p_{*}$ or exactly equal to $x_{k}.$ Roughly speaking, this method combines the advantages of the two well-known optimization methods: Gradient-descent and Gauss-Newton ones. This setting has already been implemented in several technology environments, for example in Matlab by calling the “lsqnonlin” or “fsolve” function.

In this situation, we consider the least squares problem with respect to the function $F:\mathbb{R}^{m\times p}\rightarrow\mathbb{R}^{l},$ for appropriate integer number $p,$ whose coordinate functions are defined as

[TABLE]

As mentioned in (50), the Jacobian matrix of the input function $F$ is needed:

[TABLE]

where $\dfrac{\partial F_{k}}{\partial x_{ij}}$ denotes the partial derivative of $F_{k}$ with respect to $x_{ij}.$ It is mentioned that $\dfrac{\partial F_{k}}{\partial x_{ij}}=\dfrac{\partial\phi_{k}}{\partial x_{ij}}.$

Two tolerances (small enough) below are need at each Levenberg-Marquardt iteration:

1. $\bullet$

the step tolerance $\hat{\tau}:$ at the LM step $k$th, we rely on that a numerical solution $X_{[k]}$ to system (51) exists if

[TABLE] 2. $\bullet$

the residual tolerance $\tau:$ we rely on that the least squares problem has a numerical solution if $\|F(X)\|_{2}<\tau.$

Algorithm 2. Solving the problem (1).

 

Input:

Vector $u=(u_{1},\ldots,u_{l})\in\mathbb{R}^{l}$ and function $\phi.$

Output:

a solution $X\in\mathbb{R}^{m\times p}$ to (1) .

Step 1.

Set $r=1.$

Step 2.

Solve the system $\phi(X)=u,$ by applying the gLM-method to the function defined as in (51). Let $X_{[k]}$ be a numerical solution determined by (52) with respect to the tolerance $\hat{\tau}$ at some LM-iteration $k.$ Compute $F(X_{[k]}).$

Step 3.

If $\|F(X_{[k]})\|_{2}<\tau$ then $X=X_{[k]}$ is a numerical solution and stop.

Else,

  3.1. Set $r=r+1.$

  3.2. Go to Step 2.

3.2 Affine rank minimization problem over psd-scalar matrices

3.2.1 Reformulation

As mentioned earlier, the coefficients of a sum of squares scalar-polynomial linearly depend on its Gram matrix’s entries. Finding a low rank Gram matrix of a sos scalar-polynomial hence leads to an affine rank minimization problem over positive semidefinite matrices:

[TABLE]

where $\ell:\mathbb{S}^{m}\mathbb{R}\rightarrow\mathbb{R}^{l}$ is a linear map and $b\in\mathbb{R}^{l}$ is given. It is clear that the above problem is not in form (1). It is well-known that any psd-matrix $X$ can be defined by its Cholesky factor $Y\in\mathbb{S}^{m}\mathbb{R}:$ $X=YY^{\mathrm{T}}.$ Therefore, by expressing the linear map $\ell$ via $l$ matrices $A_{i}\in\mathbb{S}^{m}\mathbb{R}:$

[TABLE]

we can define an objective function $F:\mathbb{S}^{m}\mathbb{R}\rightarrow\mathbb{R}$ that can apply LM-method [15, 20, 21] as:

[TABLE]

Moreover, if $X\in\mathbb{S}_{+}^{m}\mathbb{R}$ is needed of rank $r$ then so does $Y\in\mathbb{R}^{m\times r}.$ The problem (53) can be cast as

[TABLE]

which has the form (1). In this situation, the Jacobian matrix of $F$ can be directly computed as follows

[TABLE]

where

[TABLE]

Here “$\mathrm{vec}(X)$” means the column vector that is built by stacking the X’s columns as usual.

3.2.2 Numerical tests

It is very well-known from the theoretical point of view that the rank matrix minimization problem, say RM-problem, is NP-hard. One reason might be the nonconvexity of the rank function, so that there has not been any method directly solve this one in literature. A good way for solving the RM-problem over positive semidefinite matrices is to relax the rank function into nuclear norm function (see, eg., in [17, 14]). The resulting problemis then a semidefinite program [23] and can be efficiently solved by SDP solvers. Very recently, Huang and Wolkowicz [6] have proposed a method, that combines a facial reduction and low rank structure of semidefinite embedding, to solve a matrix completion problem. Another method for solving the NNM-problem over positive semidefinite matrices was proposed in [11] by using modified fixed point continuation method. Unfortunately, resulting matrices given by SDP solver are usually full rank.

The experiments in this section are also implemented in OCTAVE, where some codes are translated from which of COT [21] in Matlab. The input matrices $A_{1},\ldots,A_{k}$ and $b$ have entries randomly chosen in (0,1). Below we see that the resulting matrices given by gLM-method are low-rank.

4 Sosrf-based representations for some psd-polynomial matrices

4.1 Sosrf-based representations

In this section we present and prove sosrf-based and/or sos-based decompositions of a psd-polynomial matrix defined on the sets: $\mathbb{R}$ and its special subset, $\mathbb{R}^{n},$ and strips $[a,b]\times\mathbb{R}.$ The proofs are based on the Schmüdgen’s representation described in Proposition 1 and Artin’s Theorem [1] solving Hilbert’s 17th problem stating the existence of sosrf-representations of nonnegative scalar polynomials on $\mathbb{R}^{n}.$

The following result is a direct corollary of Proposition 1, its proof was presented in [19] but we would rewritten here that makes the readers more easily understand our proofs below.

Proposition 4**.**

[19, Proposition 11]* For each symmetric (polynomial) matrix $\mathbf{F}\in\mathbb{S}^{m}(\mathbb{R}[x]),$ the following statements are equivalent:*

1. (i)

$\mathbf{F}\succeq 0$* (on $\mathbb{R}^{n}$).* 2. (ii)

$b^{2}\mathbf{F}$* are sum of Hermitian squares for some polynomial $b,$ that is*

[TABLE]

where $\mathbf{A}_{i}$ are polynomial matrices.

Proof.

It is not hard to see $\mathbf{F}\succeq 0$ on $\mathbb{R}^{n}$ if $p^{2}\mathbf{F}=\sum_{i=1}^{r}\mathbf{A}_{i}^{\mathrm{T}}\mathbf{A}_{i}.$

Conversely, we consider two cases of $\mathbf{F}$ as follows:

(a) if $\mathbf{F}$ is diagonal and $\mathbf{F}=\mathbf{D}:=\mbox{\rm diag}(d_{1},\ldots,d_{m})\succeq 0$ on $\mathbb{R}^{n}.$ Then, denote $e_{i}$ the $i$th unit vector in $\mathbb{R}^{n},$ we have

[TABLE]

It thus follows from Artin’s Theorem solving Hilbert’s 17th problem that

[TABLE]

for some polynomials $p_{i},g_{ij}\in\mathbb{R}[x].$ Set $p=p_{1}\ldots p_{r},$ $h_{ij}=g_{ij}\prod\limits_{j\neq i}{p_{j}}$ and $A_{j}=\mbox{\rm diag}(h_{1j},\ldots,h_{rj})$ we have $p^{2}\mathbf{D}=\sum_{i=1}^{k}\mathbf{A}_{i}^{\mathrm{T}}\mathbf{A}_{i}.$

(b) if $\mathbf{F}$ is not diagonal then from (4), $b^{2}\mathbf{F}=\mathbf{X}_{+}\mathbf{D}\mathbf{X}_{+}^{\mathrm{T}}$ we are done by applying the representation of $\mathbf{D}$ above. ∎

The following result shows that the number of sosrf-terms in Proposition 4 can be restricted to 2 if polynomial matrices defined on $\mathbb{R}.$

Proposition 5**.**

Let $\mathbf{F}\in\mathbb{S}^{m}\mathbb{R}[x],$ $x\in\mathbb{R}.$ Then $\mathbf{F}\succeq 0$ on $\mathbb{R}$ if and only if there exists a polynomial $b\neq 0$ and two polynomial matrices $\mathbf{U},\mathbf{V}\in\mathbb{R}[x]^{m\times m}$ such that

[TABLE]

Proof.

We first prove the proposition for $\mathbf{F}$ is diagonal,

[TABLE]

where $d_{i}$’s are non-negative polynomials. One notes that each univariate polynomial which is non-negative on the real line can be expressed as a sum of at most two squares. So $d_{i}=a_{i}^{2}+b_{i}^{2},\ \forall i=1,\ldots,m$ where $a_{i},b_{i}\in\mathbb{R}[x].$ We then have

[TABLE]

The opposite direction is not hard to obtain by directly computing the summation of two diagonal matrices $\mathbf{A}^{\mathrm{T}}\mathbf{A}+\mathbf{B}^{\mathrm{T}}\mathbf{B}.$

Now, suppose $\mathbf{F}\in\mathbb{S}^{m}\mathbb{R}[x]$ is positive semidefinite on $\mathbb{R}^{n}.$ It follows from Proposition 1, one has

[TABLE]

for some $b\in\mathbb{R}[x],$ $\mathbf{X}_{-},\mathbf{X}_{+}\in\mathbb{R}[x]^{m\times m}$ and $\mathbf{D}$ is diagonal. Applying the above proof to the matrix $\mathbf{D}$ one obtains

[TABLE]

where $\mathbf{A},\mathbf{B}$ are defined as in (56), and we are done.

Conversely, if $b^{2}\mathbf{F}=\mathbf{U}^{\mathrm{T}}\mathbf{U}+\mathbf{V}^{\mathrm{T}}\mathbf{V}$ then for each $x\in\mathbb{R}$ and for all $y\in\mathbb{R}^{m}$ one has

[TABLE]

This yields $y^{\mathrm{T}}\mathbf{F}(x)y\geq 0$ for all $y\in\mathbb{R}^{m}$ and hence $\mathbf{F}(x)\succeq 0$ for all $x\in\mathbb{R}.$ ∎

On the interval $[0,+\infty)$ we have the following.

Proposition 6**.**

Let $\mathbf{F}\in\mathbb{S}^{m}\mathbb{R}[x],$ $x\in\mathbb{R}.$ Then $\mathbf{F}\succeq 0$ on $[0,+\infty)$ if and only if there exists a (scalar-coefficient) polynomial $b\neq 0$ and two polynomial matrices $\mathbf{U},\mathbf{V}\in\mathbb{R}[x]^{m\times m}$ such that

[TABLE]

Proof.

Analogously to the proof of the previous proposition, we prove the case of diagonal matrices. Suppose $\mathbf{D}=\mbox{\rm diag}(d_{1},\ldots,d_{m})\succeq 0,$ that is $d_{i}(x)\geq 0$ for all $x\in[0,+\infty).$ It is provided by the work in [2] that $d_{i}=a_{i}^{2}+b_{i}^{2}x,$ $a_{i},b_{i}\in\mathbb{R}[x],$ for all $\forall i=1,\ldots,m.$ So

[TABLE]

The “only if” part is easy to see by directly computing the multiplication of diagonal matrices.

We now consider an arbitrary $\mathbf{F}\succeq 0$ on $[0,+\infty).$ Proposition 1 leads us to the following decomposition

[TABLE]

for some $b\in\mathbb{R}[x],$ $\mathbf{X}_{-},\mathbf{X}_{+}\in\mathbb{R}[x]^{m\times m}$ and $\mathbf{D}$ is diagonal. Using the representation of $\mathbf{D}$ above we obtain that

[TABLE]

For the converse direction, take $x\geq 0$ arbitrary, then for every $y\in\mathbb{R}^{m}$ one has

[TABLE]

This implies $y^{\mathrm{T}}\mathbf{F}(x)y\geq 0$ for all $y\in\mathbb{R}^{m}.$ Thus $\mathbf{F}(x)\succeq 0$ on $[0,+\infty).$ ∎

We now consider the polynomial matrices positive semidefinite on $[a,b].$

Proposition 7**.**

Let $\mathbf{D}\in\mathbb{R}[x]^{m\times m}$ be a diagonal univariate-polynomial matrix. Then $\mathbf{D}\succeq 0$ on $[a,b]$ if and only if there are diagonal (polynomial) matrices $\mathbf{A}_{i},$ $i=1,\ldots,4$ such that

[TABLE]

Similarly to the two previous cases we are provided the following result. Its proof utilises the same techniques as in the four propositions above.

Proposition 8**.**

Let $\mathbf{F}\in\mathbb{S}^{m}\mathbb{R}[x].$ Then $\mathbf{F}\succeq 0$ on $[a,b]$ if and only if there exists a (scalar-coefficient) polynomial $b\neq 0$ and for polynomial matrices $\mathbf{U}_{1},\mathbf{U}_{2},\mathbf{U}_{3},\mathbf{U}_{4}\in\mathbb{R}[x]^{m\times m}$ such that

[TABLE]

Proof.

It is sufficient to prove for diagonal matrix $\mathbf{F}.$ It is well known that (see, eg., [2, 10]) a scalar polynomial $p$ of degree $d$ is nonnegative on $[a,b]$ if and only if there are two polynomials $p_{1},p_{2}$ such that

[TABLE]

By decomposing $\mathbf{D}=\mathbf{D}_{1}+\mathbf{D}_{2}$ where $\mathbf{D}_{1}$’s diagonal entries (resp., $\mathbf{D}_{2}$’s) are those of $\mathbf{D}$ with odd (resp., even) degree as the same position in $\mathbf{D},$ and are zeros otherwise. Then applying the above representation of polynomial nonnegative on $[a,b]$ to diagonal entries of $\mathbf{D}_{1}$ and $\mathbf{D}_{2},$ one can find $\mathbf{A}_{1},\mathbf{A}_{2}$ with respect to $\mathbf{D}_{1}$ and $\mathbf{A}_{3},\mathbf{A}_{4}$ with respect to $\mathbf{D}_{2}.$

The converse part is derived analogously to the two previous cases. ∎

On the strips $[a,b]\times\mathbb{R}$ we have the following.

Proposition 9**.**

Let $\mathbf{D}\in\mathbb{R}[x]^{m\times m}$ be a diagonal matrix defined on the strip $[0,1]\times\mathbb{R}.$ Then $\mathbf{D}\succeq 0$ on the strip $[0,1]\times\mathbb{R}$ if and only if there exist diagonal (polynomial) matrices $\mathbf{A}_{t},\mathbf{B}_{t},$ $t=1,\ldots,k$ such that

[TABLE]

Proof.

Suppose $\mathbf{D}=\mbox{\rm diag}(d_{1},\ldots,d_{m})\succeq 0$ on the strip $[0,1]\times\mathbb{R}.$ One can find in [12] that

[TABLE]

where $\sigma_{i},\tau_{i}\in\mathbb{R}[x_{1},x_{2}]$ are sums of squares:

[TABLE]

We then take $\mathbf{A}_{t}=\mbox{\rm diag}(g_{1t},\ldots,g_{kt})$ and $\mathbf{B}_{t}=\mbox{\rm diag}(h_{1t},\ldots,h_{kt})$ so that the desired representation of $\mathbf{D}$ is obtained.

The converse direction is done in the similar way to the previous proofs. ∎

Proposition 10**.**

Given $\mathbf{F}\in\mathbb{S}^{m}(\mathbb{R}[x_{1},x_{2}]).$ Then $\mathbf{F}\succeq 0$ on $[0,1]\times\mathbb{R}$ if and only if there exists a (scalar-coefficient) polynomial $b\neq 0$ and polynomial matrices $\mathbf{U}_{t},\mathbf{V}_{t}\in\mathbb{R}[x_{1},x_{2}]^{m\times m},$ $t=1,\ldots,k$ such that

[TABLE]

Proof.

This proposition is a consequence of Propositions 1 and 9. ∎

4.2 Numerical tests

Given a polynomial matrix $\mathbf{F}\in\mathbb{S}^{m}\mathbb{R}[x].$ The idea to find sos-based/sosrf-based representations of $\mathbf{F}$ is

• •

firstly apply Algorithm 1 to $\mathbf{F}$ to obtain a diagonal polynomial matrix $\mathbf{D}=\mbox{\rm diag}(d_{1},\ldots,d_{m})$ satisfying (4);

• •

next, find sosrf-representations of $d_{1},\ldots,d_{m},$ with assumption that the denominators in sosrf-representations of $d_{i}$’s are all exponents of sums of coordinates functions, by applying the algorithm in [9].

In the second stage of our computation, we will express polynomial matrices as matrix polynomials. The implementation needs to take advantage of a monomial order for them. In the rest of this paper, we prefer the lexicographic order of monomials earlier mentioned. We suppose that a $n$-variable matrix polynomial $\mathbf{P}(x)$ of degree $d,$ is always expressed under lexicographic order:

[TABLE]

With $\mathbf{P}(x)$ written as above, let furthermore $\mathcal{P}$ be the corresponding column vector of matrix coefficients $P_{\alpha},$ and $\Pi(x):=\Pi_{n,d}(x)$ be the column vector of matrices $x^{\alpha}I_{m}$ with respect to “$\leqq_{lex}$”. Then $\mathbf{P}(x)$ can be written as (see also in, eg., [16, 5])

[TABLE]

For example,

1. i)

$m=n=1,$ one has $\Omega(1,d)=\{0,1,\cdots,d\}$ and

[TABLE] 2. ii)

$m=n=d=2,$ $\Omega=\{(0,0),(0,1),(0,2),(1,0),(1,1),(2,0)\}$ and

[TABLE]

One can check that

[TABLE]

Example 3**.**

Reconsider the matrix in Example 2: $\mathbf{F}(x,y)=\begin{bmatrix}1+x^{4}y^{2}&xy\\ xy&1+x^{2}y^{4}\end{bmatrix}.$ It can be verified that $\mathbf{F}\succeq 0$ on $\mathbb{R}^{2}$ since $\mathbf{F}_{11}=1+x^{4}y^{2}>0,$ $\mathbf{F}_{22}=1+x^{2}y^{4}>0$ and

[TABLE]

Relying on the algorithm proposed in [9], we wish to obtain two diagonal polynomials of $\mathbf{D}$ from Example 2 are sums of squares of rational functions. It is clear that the 1st-one is a sum of squares. The 2nd-one is $(1+x^{4}y^{2})\det\mathbf{F}.$ We see that both $\det\mathbf{F}$ and $(1+x^{4}y^{2})\det\mathbf{F}$ consist of the term $-x^{2}y^{2}$ and they do not contain the terms

[TABLE]

But since

[TABLE]

both $\det\mathbf{F}$ and $(1+x^{4}y^{2})\det\mathbf{F}$ are not sums of squares of polynomial matrices.

Applying the algorithm proposed in [9] we obtain an approximation444 In fact we can directly compute that

$(1+x^{2}+y^{2})\det\mathbf{F}(x,y)=\frac{7}{8}+y^{2}+x^{2}+(\sqrt{2}x^{2}y^{2}-1)^{2}+x^{2}y^{6}+x^{6}(y^{2}+y^{6}+y^{8})+x^{8}y^{6}.$

of $(1+x^{2}+y^{2})\det\mathbf{F}(x,y)$ as a sum of 4 squares. A decomposition of $\mathbf{F}$ as in Proposition 5 is thus derived. It is worth mentioning that the algorithm in [9] is implemented in Matlab, we have converted these Matlab codes into OCTAVE for our present situation.

5 Conclusion

We have clearly expressed the Schmüdgen’s diagonalization for polynomial matrices and have then implemented in OCTAVE. We have then proposed an algorithm for representing a global positive semidefinite polynomial matrix as a sum of Hermitian squares of polynomial matrices. Our method combines the Schmüdgen’s diagonalization and an algorithm that approximates the diagonal elements of the resulting diagonal matrix as sums of squares, due to Reznick’s idea. This have been applied to some psd polynomial matrices on several sets such as $\mathbb{R}$ and its intervals, and strips $[a,b]\times\mathbb{R}.$ Corresponding numerical illustrations with OCTAVE have been presented in the paper. This might lead to other numerical diagonalizations of other psd polynomial matrices in the future.

Acknowledgement

The first author would like to thank Vietnam Institute for Advanced Study in Mathematics (VIASM) for warm hospitality during his visit.