Regularization Methods for Sum of Squares Relaxations in Large Scale Polynomial Optimization \author{Jiawang Nie\footnote{Department of Mathematics, University of California, 9500 Gilman Drive, La Jolla, CA 92093. Email: njw@math.ucsd.edu. The research was partially supported by NSF grants DMS-0757212, DMS-0844775 and Hellman Foundation Fellowship.}} \date{September 18, 2009}
We study how to solve sum of squares (SOS) and Lasserre's relaxations for large scale polynomial optimization. When interior-point type methods are used, typically only small or moderately large problems could be solved. This paper proposes the regularization type methods which would solve significantly larger problems. We first describe these methods for general conic semidefinite optimization, and then apply them to solve large scale polynomial optimization. Their efficiency is demonstrated by extensive numerical computations. In particular, a general dense quartic polynomial optimization with $100$ variables would be solved on a regular computer, which is almost impossible by applying prior existing SOS solvers.
0909.3551.tex
\documentclass[11pt]{article}

%\usepackage{showtags}
%\usepackage{showkeys}

\usepackage[margin=1.2in]{geometry}

%\input{../../LATEXHOME/njwlatex.tex}

\usepackage{amsmath}
\usepackage{amsthm}
\usepackage{graphicx}
%\usepackage{colortbl}
\usepackage{color}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amscd}
\usepackage{url}




\newcommand{\re}{\mathbb{R}}
\newcommand{\cpx}{\mathbb{C}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\F}{\mathbb{F}}
\newcommand{\N}{\mathbb{N}}
\newcommand{\Q}{\mathbb{Q}}
\renewcommand{\P}{\mathbb{P}}
\newcommand{\qq}{\mathbb{Q}}
\newcommand{\curL}{\cal L}
\newcommand{\curS}{\cal S}
\newcommand{\bdS}{\mathcal{S}}
\newcommand{\diag}{\mbox{diag}}
\newcommand{\dom}{\mathbf{dom}}
\newcommand{\lmda}{\lambda}
\newcommand{\half}{\frac{1}{2}}
\newcommand{\lmd}{\lambda}
\newcommand{\ntg}{\notag}
\newcommand{\nn}{\nonumber}
\newcommand{\eps}{\epsilon}
\newcommand{\veps}{\varepsilon}
\newcommand{\vareps}{\varepsilon}
\newcommand{\dt}{\delta}
\newcommand{\Dt}{\Delta}
\def\af{\alpha}
\def\bt{\beta}
\def\gm{\gamma}
\def\rank{\mbox{rank}}
\def\vec{\mathbf{vec}}
\def\mat{\mathbf{mat}}
\newcommand{\sgm}{\sigma}
\newcommand{\sig}{\sigma}
\newcommand{\Sig}{\Sigma}
%\newcommand{\proof}{\noindent {\it Proof \quad  }}
%\newcommand{\eproof}{{\hfill $\square$ }}
\newcommand{\eproof}{{\ $\square$ }}
\newcommand{\reff}[1]{(\ref{#1})}
\newcommand{\pt}{\partial}
\newcommand{\prm}{\prime}
\newcommand{\mc}[1]{\mathcal{#1}}
\newcommand{\mf}[1]{\mathfrak{#1}}
\newcommand{\td}[1]{\tilde{#1}}
\newcommand{\supp}[1]{\mbox{supp}(#1)}
\newcommand{\adj}[1]{\mathbf{adj}(#1)}




\newcommand{\bdes}{\begin{description}}
\newcommand{\edes}{\end{description}}

\newcommand{\bal}{\begin{align}}
\newcommand{\eal}{\end{align}}


\newcommand{\bnum}{\begin{enumerate}}
\newcommand{\enum}{\end{enumerate}}

\newcommand{\bit}{\begin{itemize}}
\newcommand{\eit}{\end{itemize}}

\newcommand{\bea}{\begin{eqnarray}}
\newcommand{\eea}{\end{eqnarray}}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}

\newcommand{\baray}{\begin{array}}
\newcommand{\earay}{\end{array}}

\newcommand{\bsry}{\begin{subarray}}
\newcommand{\esry}{\end{subarray}}

\newcommand{\bca}{\begin{cases}}
\newcommand{\eca}{\end{cases}}

\newcommand{\bcen}{\begin{center}}
\newcommand{\ecen}{\end{center}}


\newcommand{\bbm}{\begin{bmatrix}}
\newcommand{\ebm}{\end{bmatrix}}


\newcommand{\bmx}{\begin{matrix}}
\newcommand{\emx}{\end{matrix}}


\newcommand{\bpm}{\begin{pmatrix}}
\newcommand{\epm}{\end{pmatrix}}


\newcommand{\btab}{\begin{tabular}}
\newcommand{\etab}{\end{tabular}}



\newtheorem{theorem}{Theorem}[section]
\newtheorem{pro}[theorem]{Proposition}
\newtheorem{prop}[theorem]{Proposition}
\newtheorem{lem}[theorem]{Lemma}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{cor}[theorem]{Corollary}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{ass}[theorem]{Assumption}
\newtheorem{defi}[theorem]{Definition}
\theoremstyle{definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{exm}[theorem]{Example}
\newtheorem{alg}[theorem]{Algorithm}
\newtheorem{algorithm}[theorem]{Algorithm}
\newtheorem{remark}[theorem]{Remark}

% set up the counters in sections
\renewcommand{\theequation}{\thesection.\arabic{equation}}
\setcounter{equation}{0}
\setcounter{subsection}{0}

%
%\renewcommand{\subsection}[1]{
%    \stepcounter{subsection}
%    \settowidth{\hangindent}{\bf\thesubsection.~}
%    \hangafter=1
%    \bigskip\bigskip\noindent
%    {\bf\hbox{\thesubsection.~}#1}\par
%    \nobreak
%    \medskip
%}
%


\begin{document}

\title{Regularization Methods for Sum of Squares Relaxations in
Large Scale Polynomial Optimization
\author{Jiawang Nie\footnote{Department of Mathematics,
University of California, 9500 Gilman Drive, La Jolla, CA 92093.
Email: njw@math.ucsd.edu. The research was partially supported by NSF grants
DMS-0757212, DMS-0844775 and Hellman Foundation Fellowship.}}
\date{September 18, 2009}
}

\maketitle

\begin{abstract}
We study how to solve sum of squares (SOS) and Lasserre's relaxations
for large scale polynomial optimization.
When interior-point type methods are used,
typically only small or moderately large problems could be solved.
This paper proposes the regularization type methods
which would solve significantly larger problems.
We first describe these methods for general conic semidefinite optimization,
and then apply them to solve large scale polynomial optimization.
Their efficiency is demonstrated by extensive numerical computations.
In particular, a general dense quartic polynomial optimization with $100$ variables
would be solved on a regular computer,
which is almost impossible by applying prior existing SOS solvers.
\end{abstract}


\noindent
{\bf Key words} \, polynomial optimization, regularization methods,
semidefinite programming, sum of squares, Lasserre's relaxation

\bigskip
\noindent
{\bf AMS subject classification} \, 65K05, 90C22



\section{Introduction}

Consider the polynomial optimization
\be \label{pop:S}
\min_{x\in \re^n} \quad f(x) \quad s.t. \quad x \in S
\ee
where $f(x)$ is a multivariate polynomial and $S \subseteq \re^n$ is a semialgebraic set
(defined by a boolean combination of polynomial equalities or inequalities).
Recently there has been much work on solving \reff{pop:S}
by various sum of squares (SOS) or Lasserre's relaxations.
The basic idea is approximating nonnegative polynomials by SOS polynomials.
One of their attracting properties is
they are equivalent to certain semidefinite programming (SDP) problems.
Thus they would be solved by SDP packages
(like SDPT3 \cite{sdpt3}, SeDuMi \cite{sedumi}, SDPA \cite{sdpa}).
Numerical experiments show they give
good approximations for \reff{pop:S}, which is NP-hard.
For instance, in minimizing a general class of polynomials,
SOS relaxations work very well in finding their global minimum \cite{PS01}.
Typically SOS and Lasserre's relaxations are very successful in solving \reff{pop:S}.
We refer to \cite{JL05,Las01,Lau09,NDS,PS01,ParMP,Sweg05} for the work in this area.
However, their applications are very limited in solving big problems.
For instance, to minimize a quartic polynomial,
it is very difficult to solve its SOS relaxation on a regular computer
when it has more than $20$ variables.
%When interior points methods are used,
%large scale SOS relaxations require huge memory for storing big matrices
%and cost a huge amount of CPU time
%for solving large scale polynomial optimization.
So far, SOS relaxations can only be solved for small or moderately large problems,
which severely affects their practical applications.
The motivation of this paper is proposing cheaper numerical methods
for solving big polynomial optimization.


A standard semidefinite programming problem is
\be \label{sdp:primal}
\baray{rl}
\underset{X\in\mc{S}^N}{\min} &  C \bullet X \\
s.t. & \mc{A}(X) = b,\, X  \succeq 0.
\earay
\ee
Here $\mc{S}^N$ denotes the space of $N\times N$ real symmetric matrices,
$X\succeq 0$ (resp. $X\succ 0$) means $X$ is positive semidefinite (resp. definite),
and $\bullet$ denotes the standard Frobenius inner product.
Denote by $\mc{S}_+^N$ the cone of
positive semidefinite matrices in $\mc{S}^N$.
The $C \in \mc{S}^N$ and $b\in \re^m$ are constant,
and $\mc{A}: \, \mc{S}^N \rightarrow \re^m$ is a linear operator.
The dual of \reff{sdp:primal} is
\be  \label{sdp:dual}
\baray{rl}
\max &  b^Ty  \\
s.t. & \mc{A}^*(y) + Z = C, \, Z \succeq 0.
\earay
\ee
Here $\mc{A}^*$ is the adjoint of $\mc{A}$.
%When both \reff{sdp:primal} and \reff{sdp:dual} satisfy Slater's condition
%(existence of strictly feasible interior points),
An $X$ is optimal for \reff{sdp:primal} and $(y,Z)$ is optimal for \reff{sdp:dual}
if the triple $(X,y,Z)$ satisfies the optimality condition
\be \label{sdp:optcond}
\left.
\baray{rl}
\mc{A}(X) & = b \\
A^*(y) + Z & = C \\
 X, Z \succeq 0, \,\, XZ & = 0
\earay \right\}.
\ee
There is much work on solving SDP by interior point type methods.
Most of them generate a sequence $\{(X_k,y_k,Z_k)\}$ converging to an optimal solution.
%near the central path
%(replace $XZ=0$ by $XZ=\mu I_N$ for $\mu>0$).
%\be
%\left.
%\baray{rl}
%\mc{A}(X) & = b \\
%A^*(y) + Z & = C \\
%X \bullet Z & = \mu I \\
%X, Z & \succeq 0
%\earay \right\}.
%\ee
%As $\mu \to 0$, the sequence $\{(X_k,y_k,Z_k)\}$
%will approach an optimal triple $(X,y,Z)$.
At each step, a search direction
$(\Dt X, \Dt y, \Dt Z)$ need to be computed.
To compute $\Dt y$, typically an $m\times m$ linear system need to be solved.
To compute $\Dt X$ and $\Dt Z$, two linear matrix equations
need to be solved. Typically the cost for computing $\Dt y$ is $\mc{O}(m^3)$.
When $m = \mc{O}(N)$, the cost for computing $\Dt y$ is $\mc{O}(N^3)$.
In this case, solving SDP is not very expensive
if $N$ is not big (like less than $1,000$).
However, when $m = \mc{O}(N^2)$, the cost for computing $\Dt y$ would be $\mc{O}(N^6)$,
which is very expensive even for moderately large $N$ (like $500$).
In this case, the cost for computing $\Dt y$ is very big.
It requires the storage of a matrix of dimension $m\times m$
and $\mc{O}(m^3)$ arithmetic operations.
We refer to \cite{Tod01,VB96,WSV} for numerical methods solving SDP.



Unfortunately, SDP problems arising from polynomial optimization
belong to the latter case that $m=\mc{O}(N^2)$,
which is why the SDP packages based on interior point methods
have difficulty in solving big SOS relaxations (like degree $4$ with $100$ variables).
Now let us illustrate this.
Let $p(x)$ be a polynomial of degree $2d$.
Then $p(x)$ is SOS if and only if there exists a matrix $X \in \mc{S}_+^N$ \cite{ParMP}
such that $p(x) = [x]_d^T X [x]_d$,
where $[x]_d$ denotes the column vector of all monomials up to degree $d$ in
graded alphabetical ordering. Note the length of $[x]_d$ is $N=\binom{n+d}{d}$.
If we write $p(x)$ as
\[
p(x) = \sum_{\af\in\N^n: |\af| \leq 2d} p_\af x_1^{\af_1}\cdots x_n^{\af_n},
\]
then $p(x)$ being SOS is equivalent to the existence of $X$ satisfying
\be  \label{sos:sizes}
\baray{rcl}
A_\af \bullet X &=& p_\af \quad \forall \,\af\in\N^n: |\af| \leq 2d, \\
X &\succeq& 0.
\earay
\ee
Here $A_\af$ are certain constant symmetric matrices.
The number of equalities is $m = \binom{n+2d}{2d}$.
For any fixed $d$, $m=\mc{O}(N^2)=\mc{O}(n^{2d})$.
The size of SDP \reff{sos:sizes} is huge for moderately large $n$ and $d$.
Table~\ref{tab:(N,m)} lists the size of SDP \reff{sos:sizes} for some typical values of $(n,2d)$.
\begin{table}
\bcen
{\scriptsize
\btab[htb]{|c|c|c|c|c|c|} \hline
n= &  10 & 20 & 30 & 40 & 50  \\
2d = 4&  (66, 1001) & (231, 10626) & (496, 46376) & (861, 135751) & (1326, 316251) \\ \hline
n= &  60 & 70 & 80 & 90 & 100  \\
2d=4 & (1891, 635376) & (2556, 1150626) & (3321, 1929501) & (4186, 3049501) & (5151, 4598126) \\ \hline
n= &  10 & 15 & 20 & 25 & 30  \\
2d = 6& (286, 8008) & (816, 54264) & (1771, 230230) & (3276, 736281) & (5456, 1947792) \\ \hline
n= &  5 & 10 & 15 & 20 & 25  \\
2d = 8& (126, 1287) & (1001, 43758) & (3876, 490314) & (10626,3108105) & (23751,13884156) \\ \hline
n= \, &  5 & 8 & 9 & 10 & 15  \\
2d = 10& (252, 3003) & (1287,43758) & (2002, 92378) & ( 3003, 184756) & (15504,3268760) \\ \hline
\etab}
\caption{A list of sizes of SDP \reff{sos:sizes}.
In each pair (N,m), $N$ is the length of matrix and
$m$ is the number of equality constraints.}
\label{tab:(N,m)}
\ecen
\end{table}
As we can see, when interior point methods are applied
to solve \reff{sos:sizes},
at each step we need to solve a linear system of $m$ variables and two matrix equations.
To compute $\Delta y$, we need to store an $m\times m$ matrix
and do $\mc{O}(n^{6d})$ arithmetic operations.
This is very expensive for even moderately large $n$ and $d$,
and hence severely limits the applicability of SOS relaxations.
On a regular computer, to solve a quartic polynomial optimization,
it is quite difficult to apply interior point methods
when there are more than $20$ variables.

Recently there has been much work on designing efficient
numerical methods on solving large scale SDPs.
%where the number of equality constraints $m$ is significantly bigger than
%the matrix length $N$.
The so-called {\it regularization methods} are specially designed
for solving SDPs whose number of equality constraints $m$ is significantly bigger than
the matrix length $N$.
We refer to \cite{MPRW,PRW06,ZST08} for the work in this area.
Their numerical experiments show that
these methods are practical and efficient in solving large scale SDPs.
Currently, these methods are only designed for standard SDPs of
form \reff{sdp:primal} or its dual \reff{sdp:dual}.
However, the SDPs arising from
polynomial optimization are usually in the form of
conic semidefinite optimization, that is, they have block structures.
Therefore, we first propose the regularization methods
for solving general conic semidefinite optimization,
and then apply them to solve polynomial optimization.
Our numerical experiments show that
these new methods could solve significantly bigger polynomial optimization problems efficiently.



This paper is organized as follows.
Section~\ref{sec:regalg} describes the regularization methods
for solving general conic semidefinite optimization.
Section~\ref{sec:unopt} discusses how to solve unconstrained polynomial optimization.
Section~\ref{sec:hmgopt} discusses how to solve homogeneous polynomial optimization.
Section~\ref{sec:conopt} discusses how to solve constrained polynomial optimization.
Section~\ref{sec:cncl-dis} concludes this paper with some discussions.



\bigskip
\noindent
{\bf Notations.} \,
The symbol $\N$ (resp., $\re$) denotes the set of nonnegative integers (resp., real numbers).
For any $t\in \re$, $\lceil t\rceil$ (resp., $\lfloor t\rfloor$)
denotes the smallest integer not smaller
(resp., the largest integer not bigger) than $t$.
For $x \in \re^n$, $x_i$ denotes the $i$-th component of $x$,
that is, $x=(x_1,\ldots,x_n)$.
The $\mathbb{S}^{n-1}$ denotes the $n-1$ dimensional
unit sphere $\{x\in\re^n: x_1^2+\cdots+x_n^2=1\}$.
For $\af \in \N^n$, denote $|\af| = \af_1 + \cdots + \af_n$.
The symbol $\N_{\leq k}$ denotes the set $\{\af\in\N^n: |\af| \leq k \}$,
and $\N_{k}$ denotes $\{\af\in\N^n: |\af| = k \}$.
For each $i$, $e_i$ denotes the $i$-th standard unit vector.
The $\mathbf{1}$ denotes the vector of all ones.
For $x \in \re^n$ and $\af \in \N^n$,
$x^\af$ denotes $x_1^{\af_1}\cdots x_n^{\af_n}$.
The $[x]_d$ denotes the vector of all monomials
having degrees at most $d$ in graded alphabetical ordering, that is,
\[
[x]_d^T = [\, 1 \quad  x_1 \quad \cdots \quad x_n \quad x_1^2 \quad
x_1x_2 \quad \cdots \cdots
\quad x_1^d \quad x_1^{d-1}x_2 \quad \cdots \cdots \quad x_n^d \,],
\]
and $[x^d]$ denotes the homogeneous part of $[x]_d$ having degree $d$, that is,
\[
[x^d]^T = [\, x_1^d \quad x_1^{d-1}x_2 \quad \cdots \cdots \quad x_n^d \,].
\]
%The $\re[x]$ denotes the ring of polynomials in $(x_1,\ldots,x_n)$
%whose coefficients are real numbers,
%$\re[x]_k$ denotes the subspace of homogeneous polynomials whose degrees are equal to $k$;
%$\re[x]_{\leq k}$ denotes the subspace of polynomials
%whose degrees are at most $k$.
For a polynomial $p(x)$, $\supp{p}$ denotes the support of $p(x)$, i.e.,
the set of $\af\in\N^n$ such that the monomial $x^\af$ appears in $p(x)$.
For a finite set $S$, $|S|$ denotes its cardinality.
For a matrix $A$, $A^T$ denotes its transpose.
The $I_N$ denotes the $N\times N$ identity matrix.
%For a symmetric matrix $X$, $\lmd_{max}(X)$ and $\lmd_{min}(X)$
%denote the maximum and minimum eigenvalues of $X$ respectively.
For any vector $u\in \re^N$, $\| u \|_2 = \sqrt{u^Tu}$ denotes the standard Euclidean norm.





\section{Regularization methods for conic semidefinite optimization}
\setcounter{equation}{0}
\label{sec:regalg}

This section describes the regularization methods
for solving conic semidefinite optimization problems having block diagonal structures.
They are natural generalizations of the regularizatioin methods introduced in
\cite{MPRW,PRW06,ZST08} for solving standard SDP problems.



\subsection{Conic semidefinite optimization}

Let $\mc{K}$ be a cross product of semidefnite cones
\[
\mc{K} = \mc{S}^{N_1}_+ \times \cdots \times \mc{S}^{N_\ell}_+.
\]
It belongs to the space
$
\mc{M} = \mc{S}^{N_1} \times \cdots \times \mc{S}^{N_\ell}.
$
Each $X \in \mc{M}$ is a tuple  $X=(X_1,\ldots,X_\ell)$
with every $X_i\in \mc{S}^{N_i}$.
So $X$ could be thought of as a symmetric block diagonal matrix.
Thus $X\in\mc{K}$ if and only if its every block $X_i \succeq 0$.
The notation $X \succeq_{\mc{K}} 0$ (resp. $X \succ_{\mc{K}} 0$)
means every block of $X$ is positive semidefinite (resp. definite).
For $X=(X_1,\ldots,X_\ell) \in \mc{M}$ and $Y=(Y_1,\ldots,Y_\ell) \in \mc{M}$,
define their inner product as
$
X \bullet Y =  X_1 \bullet Y_1 + \cdots + X_\ell \bullet Y_\ell.
$
Denote by $\| \cdot \|$ the norm in $\mc{M}$ induced by this inner product.
%This inner product equips $\mc{M}$ with the norm
%\[
%\|X\| =  \Big(Trace (X_1^2) + \cdots + Trace(X_\ell^2) \Big)^{1/2}.
%\]
Note $\mc{K}$ is a self-dual cone, that is,
\[
\mc{K}^*= \{ Y\in \mc{M}: Y \bullet X \geq 0 \quad \forall X \in \mc{K} \} = \mc{K}.
\]


For a symmetric matrix $W$,
denote by $(W)_+$ (resp. $(W)_{-}$) the projection of $W$
into the positive (resp. negative) semidefinite cone,
that is, if $W$ has spectral decomposition
\[
W = \sum_{\lmd_i>0} \lmd_i u_iu_i^T + \sum_{\lmd_i< 0} \lmd_i u_iu_i^T,
\]
then $(W)_+$ and $(W)_{-}$ are defined as
\[
(W)_+ = \sum_{\lmd_i>0} \lmd_i u_iu_i^T, \quad
(W)_{-} = \sum_{\lmd_i< 0} \lmd_i u_iu_i^T.
\]
For $X=(X_1,\ldots,X_\ell) \in \mc{M}$,
its projections into $\mc{K}$ and $-\mc{K}$ are given by
\[
(X)_{\mc{K}} =((X_1)_+,\ldots,(X_\ell)_+), \quad
(X)_{-\mc{K}} =((X_1)_{-},\ldots,(X_\ell)_{-}).
\]
The projection function $(\cdot)_{\mc{K}}$ is Lipschitz continuous,
because for any $X,Y \in \mc{K} $
\[
\|X_{\mc{K}}-Y_{\mc{K}}\|^2 = \|X-Y\|^2 +2(X_{\mc{K}}\bullet Y_{-\mc{K}} + Y_{\mc{K}}\bullet X_{-\mc{K}})
- \|X_{-\mc{K}}-Y_{-\mc{K}}\|^2
\leq  \|X-Y\|^2.
\]
The projection norm function $\|(\cdot)_\mc{K}\|$ is convex.
To see this, note
\[
\|(X)_{\mc{K}}\| = \min_{F \in \mc{K}}  \| X+F\|.
\]
Then, for any two $X,Y \in \mc{M}$, it holds
\begin{align*}
\left\|\left(\frac{X+Y}{2}\right)_{\mc{K}}\right\|
&\leq \left\| \frac{X+Y}{2} - \frac{(X)_{-\mc{K}}+(Y)_{-\mc{K}}}{2} \right\|
= \left\| \frac{X-(X)_{-\mc{K}}}{2} + \frac{Y - (Y)_{-\mc{K}}}{2} \right\|  \\
& \leq   \left\| \frac{X-(X)_{-\mc{K}}}{2} \right\| +
 \left\| \frac{Y - (Y)_{-\mc{K}}}{2} \right\|
 = \half \|(X)_{\mc{K}}\| + \half \|(Y)_{\mc{K}}\|.
\end{align*}


A general conic semidefinite optimization problem is
\be \label{Ksdp:prm}
\baray{rl}
\min &  C \bullet X \\
s.t. & \mc{A}(X) = b, \,  X  \in \mc{K}.
\earay
\ee
Here $C \in \mc{M}$, $b\in\re^m$,
and $\mc{A}: \mc{M} \to \re^m$ is a linear operator.
The dual of \reff{Ksdp:prm} is
\be \label{Ksdp:dual}
\baray{rl}
\max &  b^Ty \\
s.t. & \mc{A}^*(y) + Z =C, \, Z \in \mc{K}.
\earay
\ee
If $\mc{A}(X)$ is a summation like
\[
\mc{A}(X) = \mc{A}_1(X_1) + \cdots + \mc{A}_\ell(X_\ell),
\]
where every $\mc{A}_i: \mc{S}^{N_i} \to \re^m$
is a linear operator,
then its adjoint $\mc{A}^*(y)$ has the form
\[
\mc{A}^*(y) = (\mc{A}_1^*(y), \ldots, \mc{A}_\ell^*(y)).
\]



\subsection{Regularization methods}


There are two typical regularizations for a standard SDP:
{\it Moreau-Yosida regularization} for the primal \reff{sdp:primal}
and {\it Augmented Lagrangian regularization} for the dual \reff{sdp:dual}.
They would be naturally generalized to
conic semidefinite optimization \reff{Ksdp:prm} and its dual \reff{Ksdp:dual}.
The Moreau-Yosida regularization for \reff{Ksdp:prm} is
\be \label{sdp:MYreg}
\baray{rl}
\underset{X,Y \in \mc{M} }{\rm min} &  C \bullet X + \frac{1}{2\sigma} \|X-Y\|^2 \\
s.t. & \mc{A}(X) = b,\, X  \in \mc{K}.
\earay
\ee
Obviously \reff{sdp:MYreg} is equivalent to \reff{Ksdp:prm},
because for each fixed feasible $X \in \mc{M}$
the optimal $Y\in\mc{M}$ in \reff{sdp:MYreg} is equal to $X$.
The Augmented Lagrangian regularization for \reff{Ksdp:dual} is
\be \label{sdp:augLag:1st}
\baray{rl}
\underset{y\in\re^m, Z\in \mc{M}}{\rm max} &
b^Ty - (Z+\mc{A}^*(y)-C)\bullet Y - \frac{\sigma}{2} \|Z+\mc{A}^*(y) -C\|^2 \\
s.t. &  Z \in \mc{K}.
\earay
\ee
When $\mc{K}=\mc{S}_+^N$ is a single semidefinite cone,
it was shown (see Section~2 of \cite{MPRW}) that
for every fixed $Y$, \reff{sdp:augLag:1st} is the dual of the following optimization
\[
\min_{X\in\mc{M} } \quad  C \bullet X + \frac{1}{2\sigma} \|X-Y\|^2 -y^T( \mc{A}(X)-b) - Z \bullet X.
\]
By fixing $y\in\re^m$ and optimizing over $Z\succeq 0$,
Malick, Povh, Rendl, and Wiegele \cite{MPRW} further showed that
\reff{sdp:augLag:1st}
can be reduced to
\be  \label{sdp:augLag}
\baray{rl}
\underset{y\in\re^m}{\rm max} &
 b^Ty - \frac{\sigma}{2} \|(\mc{A}^*(y)-C+Y/\sigma)_{\mc{K}}\|^2
+ \frac{1}{2\sigma} \|Y\|^2.
\earay
\ee


\begin{prop} %[Proposition~2.6, \cite{MPRW}]
\label{optcond:mprw}
Assume there exists $\bar{X}\succ_{\mc{K}} 0$ such that $\mc{A}(\bar{X})=b$.
The following statements are equivalent:
\bit
\item [(i)] $(X,y,Z)$ is optimal for \reff{sdp:MYreg}-\reff{sdp:augLag:1st}.

\item [(ii)] $(X,y,Z)$ satisfies
\[
\left\{\baray{l}
\mc{A}(X) = b, \, X = \sigma (Z+\mc{A}^*(y)-C)+Y, \\
X\in \mc{K},\, Z\in\mc{K}, \, X\bullet Z = 0.
\earay \right.
\]

\item [(iii)] $(X,y,Z)$ satisfies
\be \label{bpm:optcond}
\left\{\baray{l}
X = \sigma(Y/\sigma+\mc{A}^*(y)-C)_{\mc{K}}, \\
Z = -(Y/\sigma+\mc{A}^*(y)-C)_{-\mc{K}}, \\
\mc{A}\mc{A}^* y + \mc{A}(Z-C) = (b-\mc{A}(Y))/\sigma.
\earay \right.
\ee
\eit
\end{prop}

\noindent
The proof of Proposition~2.6 in \cite{MPRW} would be naturally generalized
from a single semidefinite cone $\mc{S}_+^N$ to the cone $\mc{K}$
of block diagonal positive semidefinite matrices.
So we omit the proof of Proposition~\ref{optcond:mprw}.



In \reff{bpm:optcond}, $\mc{A}\mc{A}^*$ is a symmetric matrix defined in the way that
\[
(\mc{A}\mc{A}^* y)^Tz = \mc{A}^*(y) \bullet \mc{A}^*(z) \quad
\forall \, y,z\in \re^m.
\]
The optimality condition~\reff{bpm:optcond}
leads to the following algorithm for solving \reff{Ksdp:prm}:

\begin{alg} \label{BPM:Ksdp}
Choose parameters $\sigma>0$,
$\eps>0$, $Y\in \mc{S}^N$,  set $Z=0$. Do the loop: \\
{\it While} ($\|b-\mc{A}(X)\| \geq \eps$ or $\|C-Z-\mc{A}^*(y)\| \geq \eps$) \\
\indent \quad  Solve $\mc{A}\mc{A}^*y = \mc{A}(C-Z) + (b-\mc{A}(Y))/\sigma$ for $y$.\\
\indent \quad  Compute the projections
\[
X = \sigma(Y/\sigma+\mc{A}^*(y)-C)_{\mc{K}}, \quad
Z = -(Y/\sigma+\mc{A}^*(y)-C)_{-\mc{K}}.
\]
\indent \quad  Set $Y:=X$ and update $\sigma$.
\end{alg}
When $\mc{K}=\mc{S}_+^N$ is a single semidefinite cone,
Algorithm~\ref{BPM:Ksdp} reduces to the famous
{\it Boundary Point Method (BPM)}, which was originally
proposed by Povh, Rendl, and Wiegele \cite{PRW06}
(also see \cite{MPRW}) for solving big SDPs.
This method was proven surprisingly efficient in some applications,
as illustrated in \cite{MPRW,PRW06}.
Algorithm~\ref{BPM:Ksdp} can be thought of as a natural generalization
of BPM to general conic semidefinite optimization \reff{Ksdp:prm}.
\iffalse
% the following is dropped since we do not use Algorithm \ref{BPM:Ksdp} much in the paper
%
An attractive property of Algorithm \ref{BPM:Ksdp} is that at every iteration
we only need to solve a linear system of type $\mc{A}\mc{A}^* y = w$
and compute two projections. Before the starting of iterations,
we can compute the Cholesky factorization of $\mc{A}\mc{A}^*$
and save it for future use.
Then solving $\mc{A}\mc{A}^*y = w$ can be done by solving two triangular systems,
whose cost is $\mc{O}(m^2)$ instead of $\mc{O}(m^3)$.
This saves a lot of work when $m$ is huge.
Of course, there is a one-time $\mc{O}(m^3)$ cost
for computing Cholesky factorization. For unconstrained polynomial optimization,
the resulting $\mc{A}$ satisfies certain orthogonality properties
($\mc{A}\mc{A}^*$ is diagonal),
which makes solving $\mc{A}\mc{A}^* y = w$ is trivial.
For homogeneous polynomial optimization,
there is a neat formula for the inverse $(\mc{A}\mc{A}^*)^{-1}$,
which makes solving $\mc{A}\mc{A}^* y = w$ is also trivial.
This is a big advantage for solving SDPs arising from  polynomial optimization.
\fi
It is a specification of
the following general regularization algorithm:

\begin{alg} \label{regalg:MPRW}
Choose $Y\in \mc{S}^N$, $y\in \re^m$, and
$\eps^{in},\eps^{out}$,
set $Z=0$ and do the loop \\
\indent  Repeat until $\|Z+\mc{A}^*(y)-C\| \leq \eps^{out}$ ({\it outer loop}): \\
\indent \qquad  Repeat until $\|b-\mc{A}(X)\| \leq \eps^{in}$ ({\it inner loop}): \\
\indent \qquad \qquad  Compute the projections \\
\indent \qquad \qquad  \qquad
$X = \sigma(Y/\sigma+\mc{A}^*(y)-C)_{ \mc{K} },\,
Z = -(Y/\sigma+\mc{A}^*(y)-C)_{-\mc{K}}. $\\
\indent \qquad \qquad Set $g=b-\mc{A}(X)$.\\  \qquad
\indent \qquad \qquad
Update $y\leftarrow y + \tau W g$ with appropriate $\tau$ and $W$. \\
\indent \qquad end   \\
\indent \qquad  Set $Y:=X$ and update $\sigma$. \\
\indent end
\end{alg}

Algorithm~\ref{regalg:MPRW} is also a natural generalization of
the regularization framework proposed by Malick, Povh, Rendl, and Wiegele \cite[Algorithm~4.3]{MPRW}
for solving large scale SDP when $\mc{K}=\mc{S}_+^N$ is a single semidefinite cone.
In that case, if $W$ is chosen to be $(\mc{A}\mc{A}^*)^{-1}$
and $\tau = 1/\sigma$, Algorithm~\ref{regalg:MPRW}
reduces to Algorithm~\ref{BPM:Ksdp}, as shown in \cite[Prop.~4.4]{MPRW}.
Its convergence was discussed in \cite[Theorem~4.4]{MPRW}.


For solving large scale SOS relaxations,
Algorithm~\ref{BPM:Ksdp} might converge fast at the beginning,
but generally has slow convergence when we get close to optimal solutions.
This is because it is basically a gradient type method.
To get faster convergence,
we need better solvers for the inner loop of Algorithm~\ref{regalg:MPRW}.
Typically, Newton type methods have good properties
like local superlinear or quadratic convergence.
This leads to another kind of regularization methods
called {\it Newton-CG Augmented Lagrangian},
which was proposed by Zhao, Sun and Toh \cite{ZST08}.
It gives surprisingly good numerical performance in solving big standard SDP problems.
In the rest of this section, we generalize this important method
to solve general conic semidefinite optimization \reff{Ksdp:prm}.

Denote by $\varphi_{\sigma}(Y,y)$ the objective in \reff{sdp:augLag}.
%Define function
%\be \label{def:varphi}
%\varphi_{\sigma}(Y,y):=b^Ty - \frac{\sigma}{2} \|(\mc{A}^*(y)-C+Y/\sigma)_{\mc{K}}\|^2
%+ \frac{1}{2\sigma} \|Y\|^2.
%\ee
When $\mc{K}=\mc{S}_+^N$,
$\varphi_{\sigma}(Y,y)$ is differentiable \cite[Proposition~3.2]{MPRW} and
\[
\nabla_y \varphi_{\sigma}(Y,y) =
b - \sigma \mc{A}\Big((\mc{A}^*(y)-C+Y/\sigma)_{\mc{K}}\Big).
\]
The above is also true when $\mc{K}$ is a cross product of semidefnite cones.
The inner loop of Algorithm~\ref{regalg:MPRW} is
essentially solving the maximization problem
\be  \label{max:varphi}
\max_{y\in\re^m}  \quad \varphi_{\sigma}(Y_k,y).
\ee
Since $\varphi_{\sigma}$ is concave in $y$,
a point $\hat{y}$ is a maximizer of \reff{max:varphi}
if and only if
\[
\nabla_y \varphi_{\sigma}(Y_k,\hat{y})=0.
\]
The function $\varphi_{\sigma}(Y,y)$ is not twice differentiable,
so the standard Newton's method is not applicable.
However, the function $\varphi_{\sigma}(Y,y)$ is semismooth,
and semismooth Newton's method could be applied
to get local superliner or quadratic convergence,
as pointed out by Zhao, Sun and Toh \cite[Section~3.2]{ZST08}.
So we need to compute the generalized Hessian of $\varphi_\sig$.
If
\[
\mc{A}^*(y) = ( \mc{A}_1^*(y), \ldots, \mc{A}_\ell^*(y) ), \quad
C =( C_1, \ldots, C_\ell),
\]
then for $Y=(Y_1,\ldots,Y_\ell)\in\mc{M}$
\[
\mc{A}^*(y)-C+Y/\sigma =
\Big(\mc{A}_1^*(y)-C_1+Y_1/\sigma, \, \ldots, \, \mc{A}_\ell^*(y)-C_\ell+Y_\ell/\sigma\Big).
\]
So we get
$\varphi_{\sigma}(Y,y) = b^Ty + \frac{1}{2\sigma} \|Y\|^2
 - \sum_{k=1}^\ell \varphi_{\sigma}^{(k)}(Y,y)$, where
\be
\varphi_{\sigma}^{(k)}(Y,y) =
\frac{\sigma}{2} \|(\mc{A}_k^*(y)-C_i+Y_k/\sigma)_+\|^2.
\ee
This implies the formula
\[
\nabla_y ^2\varphi_{\sigma}(Y,y) =
 - \sum_{k=1}^\ell  \nabla_y^2 \varphi_{\sigma}^{(k)}(Y,y).
\]
A numerical method for evaluating $\nabla_y^2 \varphi_{\sigma}^{(k)}(Y,y)$
was given in \cite[Section~3.2]{ZST08}.
This can be done by computing eigenvalue decompositions. Let
\[
\mc{A}_k^*(y)-C_k+Y_k/\sigma = Q_k \Lambda^{(k)} Q_k^T,
\]
where $Q_k\in \re^{N_k\times N_k}$ is orthogonal and
$\Lambda^{(k)} = \diag(\lmd_1^{(k)}, \ldots,\lmd_{N_k}^{(k)})$ is diagonal with
$
\lmd_1^{(k)} \geq \lmd_2^{(k)} \geq \cdots \geq \lmd_N^{(k)}.
$
Denote index sets
\[
\Gamma_+^{(k)} = \{ 1 \leq i \leq N_k :  \lmd_i^{(k)} \geq 0 \}, \quad
\Gamma_-^{(k)} = \{ 1 \leq i \leq N_k :  \lmd_i^{(k)} < 0 \}.
\]
Then define a symmetric matrix $\Omega^{(k)} \in \re^{N_k\times N_k}$ as follows:
\[
\Omega^{(k)}(i,j) = \bca
1  & \text{ if }  i,j \in \Gamma_+^{(k)}, \\
\frac{\lmd_i}{\lmd_i - \lmd_j}  & \text{ if }
 i \in \Gamma_+^{(k)}, \, j\in \Gamma_-^{(k)}, \\
\frac{\lmd_j}{\lmd_j - \lmd_i}  & \text{ if }
 i \in \Gamma_-^{(k)}, \, j\in \Gamma_+^{(k)}, \\
0 & \text{ otherwise.}
\eca
\]
It was shown \cite[Section~3.2]{ZST08} that
\[
\nabla_y^2 \varphi_{\sigma}^{(k)}(Y,y) \cdot z =
\sigma  \mc{A}_k \Big( Q_k\big(\Omega^{(k)} \circ
(Q_k^T \mc{A}_k^*(z) Q_k) \big)Q_k^T \Big),
\]
where $\circ$ denotes Hadamard product.
Now we define a matrix $\mc{L}_{Y,y}^{\sigma} \in \re^{m\times m}$ satisfying
\[
\mc{L}_{Y,y}^{\sigma} \cdot z = \sigma \sum_{k=1}^\ell
\mc{A}_k \Big( Q_k\big(\Omega^{(k)} \circ (Q_k^T \mc{A}_k^*(z) Q_k) \big)Q_k^T \Big).
\]
It was shown \cite[Prop.~3.2]{ZST08} that we always have $\mc{L}_{Y,y}^{\sigma} \succeq 0$.
Furthermore, $\mc{L}_{Y,y}^{\sigma} \succ 0$ if some nondegeneracy condition holds.
In either case, an approximate semismooth Newton direction $d_{new}$
for \reff{max:varphi}
can be determined by the linear system
\be  \label{linsys:Newton}
\Big(\mc{L}_{Y,y}^{\sigma} + \eps \cdot I_N \Big)  d_{new} = \nabla_y  \varphi_\sigma.
\ee
Here $\eps>0$ is a tiny number ensuring the positive definiteness of
the above linear system.


When $m$ is huge, it is usually not practical to solve \reff{linsys:Newton}
by direct methods like Cholesky factorization.
To avoid this difficulty, Zhao, Sun and Toh \cite{ZST08} proposed applying
conjugate gradient (CG) iterations to solve \reff{linsys:Newton} approximately.
To fasten the convergence of CG, a good preconditioner for \reff{linsys:Newton}
is $(\mc{A}\mc{A}^*)^{-1}$, as mentioned in \cite{ZST08}.
In SOS relaxations for unconstrained and homogeneous polynomial optimization,
there is a nice property that the resulting $(\mc{A}\mc{A}^*)^{-1}$ has explicit neat formulae,
as will be shown in Sections~\ref{sec:unopt} and \ref{sec:hmgopt}.
So we use $(\mc{A}\mc{A}^*)^{-1}$ as a preconditioner in this situation.
However, in Lasserre's relaxation for constrained polynomial optimization,
typically there is no explicit expression for $(\mc{A}\mc{A}^*)^{-1}$,
and usually computing $(\mc{A}\mc{A}^*)^{-1}$ is quite difficult.
In this case, we just use $\big(\diag (\mc{A}\mc{A}^*) \big)^{-1}$
as a preconditioner, as suggested in \cite{ZST08}.
In either case, we set a cap $K$ on the maximum number of CG iterations
to insure solving  \reff{linsys:Newton} do not take too much time.
A standard preconditioned CG algorithm (like Algorithm~6.12 in \cite{demmel})
can be applied to solve \reff{linsys:Newton}.
The numerical experiments in \cite{ZST08} show that
this method is surprisingly powerful in
solving some big standard SDPs when $\mc{K}=\mc{S}_+^N$ is a single semidefinite cone.


Now we present the Newton-CG Augmented Lagrangian regularization method
for solving general conic semidefinite optimization \reff{Ksdp:prm}-\reff{Ksdp:dual}.

\begin{alg} \label{alg:Newton-CG}
Choose $\eps^{in},\eps^{out} \in (0,1)$,
$\eps >0$, $\dt\in(0,1)$, $\rho>1$, $\sigma_{\max},K\in \N$. \\
\indent Initialization: Choose $X_0,Z_0\in \mc{S}^N$,
$y_0\in \re^m$, $\sig_0$ and set $k=0$. \\
\indent While ($\|Z_k+\mc{A}^*(y)_k-C\| \leq \eps^{out}$) ({\large \it outer loop}): \\
\indent  \indent  Set $Y_k:=X_k$. \\
\indent  \indent  Set $j=0$ and $y_{k,j} = y_k$. \\
\indent \indent While ($\|\nabla_y \varphi_{\sigma_k}(Y_k,y_{k,j})\| \leq \eps^{in}$)
({\large \it inner loop}): \\
\indent  \indent \indent Set $g_j = \nabla_y \varphi_{\sigma_k}(Y_k,y_{k,j})$. \\
\indent  \indent \indent Compute $d_{new}$ from \reff{linsys:Newton}
by applying preconditioned CG at most $K$ steps. \\
\indent  \indent \indent  Find the smallest integer $\af>0$ such that
\indent  \indent \indent  \indent \indent \indent
\be  \label{cond:linsch}
\varphi_{\sigma_k}(Y_k, y_{k,j}+ \dt^\af \cdot d_{new} ) \geq
\varphi_{\sigma_k}(Y_k, y_{k,j}) + \dt^\af \cdot g_j^T d_{new}.
\ee
\indent  \indent \indent  Set $y_{k,j+1} := y_{k,j} + \dt^\af \cdot d_{new}$. \\
\indent \indent  \indent  Compute the projections into the cone $\mc{K}$: \\
\indent \indent  \indent  \indent
$X_k = \sigma(Y_k/\sigma+\mc{A}^*(y_{k,j+1})-C)_{\mc{K}},\,
Z_k = -(Y_k/\sigma+\mc{A}^*(y_{k,j+1})-C)_{-\mc{K}}$.\\
\indent \indent  \indent  Set $j:=j+1$. \\
\indent \indent  end ({\large \it inner loop}) \\
\indent  \indent  Set $y_{k+1} = y_{k,j}$ \\
\indent  \indent  If $\sigma_k \leq \sigma_{\max}$, set $\sigma_{k+1} = \rho \sigma_k$. \\
\indent  \indent  Set $k:=k+1$. \\
\indent end ({\large \it outer loop})
\end{alg}


When $\mc{K}=\mc{S}_+^N$ is a single semidefinite cone,
the convergence of Algorithm~\ref{alg:Newton-CG}
has been discussed in \cite[Theorems~3.5, 4.1,4.2]{ZST08}.
These results could be generalized naturally
to the cone $\mc{K}$ of block diagonal positive semidefinite matrices.
%To get the convergence of the outer loop,
%usually one of the following three criteria
%(cf. Section~4 of \cite{Roc76MOR})
%\[ \qquad \qquad
%\baray{rclr}
%\varphi_k(y_{k+1}) - \inf \varphi_k &\leq & \eps_k^2/2\sigma_k,
%\quad \eps_k\geq 0, \sum_{k=0}^\infty \eps_k < \infty & \qquad (A^{\prm\prm}) \\
%\varphi_k(y_{k+1}) - \inf \varphi_k &\leq & (\dt_k^2/2\sigma_k) \|Y_{k+1}-Y_k \|^2,
%\quad \dt_k \geq 0, \sum_{k=0}^\infty \dt_k < \infty &  \qquad  (B_1^{\prm\prm}) \\
%\| \nabla \varphi_k(y_{k+1}) \| &\leq &
%(\dt_k^\prm/2\sigma_k) \|Y_{k+1}-Y_k \|^2, \quad 0\leq \dt_k^\prm \to 0
%&  \qquad  (B_2^{\prm\prm})
%\earay
%\]
%should be satisfied.
The specifics about the convergence are beyond the scope of this paper.
We refer to \cite{MPRW,Roc76CO,Roc76MOR,ZST08}.



Algorithm~\ref{BPM:Ksdp} and \ref{alg:Newton-CG}
are basically specifications of Algorithm~\ref{regalg:MPRW}.
Generally Algorithm~\ref{BPM:Ksdp} has slower convergence than Algorithm~\ref{alg:Newton-CG} does,
because it is a gradient type method while the latter is a Newton type method.
They have similar performance when computing eigenvalue decompositions is not expensive.
However, when eigenvalue decompositions are expensive to compute,
Algorithm~\ref{BPM:Ksdp} usually takes more time
because it requires computing more eigenvalue decompositions.
In large scale polynomial optimization, the matrices often have big dimension,
and hence Algorithm~\ref{alg:Newton-CG} is more suitable.
So in later sections, the numerical experiments are mainly based on using Algorithm~\ref{alg:Newton-CG}.





\section{Unconstrained polynomial optimization}
\setcounter{equation}{0}
\label{sec:unopt}

Consider unconstrained optimization problem
\be \label{uc:minf(x)}
f_{min}^{uc}:=
\min_{x\in \re^n} \quad f(x)
\ee
where $f(x)$ is a multivariate polynomial of degree $2d$.
Here $f_{min}^{uc}$ denotes the global minimum of $f(x)$ over $\re^n$.
Problem \reff{uc:minf(x)} is NP-hard when $\deg(f)\geq 4$, as shown in \cite{Nes00}.
%Many difficult problems can be formulated in the form \reff{uc:minf(x)}.
%For instance, checking co-positivity of symmetric matrices, which is NP-hard,
%is equivalent to checking nonnegativity of certain quartic polynomials.
There has much recent work on solving \reff{uc:minf(x)} via SOS relaxations,
like \cite{JL05,Las01,Las06b,NDS,ParMP,PS01}.
The standard SOS relaxation for \reff{uc:minf(x)} is
\be  \label{opt:stdsos}
\baray{rl}
f_{sos}^{uc}:=
\max & \gamma \\
s.t. & f(x) - \gamma \mbox{ is SOS}.
\earay
\ee
Obviously the above optimal value $f_{sos}^{uc}$ satisfies
the relation $f_{sos}^{uc} \leq f_{min}^{uc}$.
Though it is possible that $f_{sos}^{uc} < f_{min}^{uc}$,
it was observed in \cite{PS01,ParMP} that
\reff{uc:minf(x)} works very well in practice.
%It almost always solves \reff{uc:minf(x)} successfully
%when $f(x)$ is not constructed artificially.
Typically \reff{opt:stdsos} gives a reasonably well
approximation for \reff{uc:minf(x)}, which is good enough in many applications.
SOS relaxations are usually very expensive to solve
when interior-point methods are used.
In this section, we show how to apply regularization methods described in Section~\ref{sec:regalg}
to solve \reff{opt:stdsos}.
We begin with the SDP formulation of \reff{opt:stdsos}.



Note $f(x)-\gamma$ is SOS if and only if there exists $X\in\mc{S}^{\binom{n+d}{d}}$
satisfying
\[
f(x)-\gamma = [x]_d^TX[x]_d = X \bullet ([x]_d[x]_d^T), \quad X \succeq 0.
\]
Here $\binom{n+d}{d}$ is the total number of monomials in $x$ of degrees up to $d$.
Define 0/1 constant symmetric matrices $C$ and $A_\af$ in the way that
\be \label{def:ucpop-A}
[x]_d[x]_d^T = C + \sum_{\af\in\N^n: 0< |\af| \leq 2d} A_\af x^\af.
\ee
Denote $\mathbb{U}_{2d}^n=\{\af\in \N^n: 0 < |\af| \leq 2d\}$. If
$
f(x) = f_0+\sum_{\af\in \mathbb{U}_{2d}^n } f_\af x^\af,
$
then $\gamma$ is feasible for \reff{opt:stdsos}
if and only if $X$ satisfies
\[
\baray{rcl}
C \bullet X + \gamma  & = & f_0, \\
A_\af \bullet X & = & f_\af \quad \forall \, \af\in\mathbb{U}_{2d}^n, \\
X & \succeq & 0.
\earay
\]
Let $b=(f_\af)_{\af\in\mathbb{U}_{2d}^n}$ be a vector of dimension
$m =\binom{n+2d}{2d}-1$. Define $\mc{A}$ and $\mc{K}$ as
\[
\mc{A}(X) = (A_\af \bullet X)_{\af\in\mathbb{U}_{2d}^n}, \quad
\mc{K} = \mc{S}_+^{\binom{n+d}{d}}.
\]
Then SOS relaxation \reff{opt:stdsos} is equivalent to
\be  \label{sdp:stdsos}
\baray{rl}
f_{sdp}^{uc} :=
\min &  C \bullet X \\
s.t. & \mc{A}(X) = b, \, X  \in \mc{K}.
\earay
\ee
The dual of the above is
\be  \label{sdp:momsos}
\baray{rl}
\max &  b^Ty  \\
s.t. & \mc{A}^*(y) + Z =C, Z \in \mc{K}.
\earay
\ee
Here the adjoint $\mc{A}^*(y) = \sum_{\af \in \mathbb{U}_{2d}^n } y_\af A_\af$.



\begin{theorem} \label{thm-bd:uc}
For notations defined above, it holds
$f_{sos}^{uc} = -f_{sdp}^{uc}+f_0 \leq f_{min}^{uc}$.
If SOS relaxation \reff{opt:stdsos} is exact,
then $f_{min}^{uc} = -f_{sdp}^{uc}+f_0$.
\end{theorem}
\begin{proof}
As we have seen in the above,
optimization problems \reff{opt:stdsos} and \reff{sdp:stdsos}
are equivalent. Their objectives satisfy the relation
\[
C \bullet X = -\gamma + f_0.
\]
Maximizing $\gamma$ in \reff{opt:stdsos} is equivalent to
minimizing $C\bullet X$ in \reff{sdp:stdsos}.
The feasible sets of \reff{opt:stdsos} and \reff{sdp:stdsos}
are the same.
Hence we have $f_{sos}^{uc} = -f_{sdp}^{uc}+f_0$.
Furthermore, $f_{sos}^{uc}\leq f_{min}^{uc}$ is obvious
because \reff{opt:stdsos} is a relaxation of \reff{uc:minf(x)}.

If SOS relaxation \reff{opt:stdsos} is exact,
i.e., $f_{sos}^{uc} = f_{min}^{uc}$, then obviously
$f_{min}^{uc} = -f_{sdp}^{uc}+f_0$.
\end{proof}



Now we turn to finding $(\mc{A}\mc{A}^*)^{-1}$.
Sine $A_\af$ defined in \reff{def:ucpop-A} are orthogonal to each other,
the matrix $\mc{A}\mc{A}^*$ is diagonal.
%$
%\mc{A}\mc{A}^* = \diag( \mathbf{1}^T A_\af \mathbf{1} : \af \in \mathbb{U}_{2d}^n ),
%$
So we get the formula
\be  \label{uc-def:u}
(\mc{A}\mc{A}^*)^{-1} = \diag( \mathbf{a} ),
\quad \mbox{ where } \quad
\mathbf{a} = ( \,(\mathbf{1}^T A_\af \mathbf{1})^{-1} :
\af \in \mathbb{U}_{2d}^n ).
\ee
Here $\mathbf{1}$ is the vector of all ones.

Algorithm~\ref{alg:Newton-CG} can be applied to solve the SDP \reff{sdp:stdsos}-\reff{sdp:momsos}.
We use $(\mc{A}\mc{A}^*)^{-1}$ via the formula \reff{uc-def:u} as a preconditioner.
Suppose $(X^*,y^*,Z^*)$ is optimal for \reff{sdp:stdsos}-\reff{sdp:momsos}.
Then $-f_{sdp}^{uc}+f_0$ is a lower bound for the minimum $f_{min}^{uc}$.
Furthermore, if $Z^*$ has rank one, then $f_{min}^{uc}=f_{sos}^{uc}$
and a global minimizer for \reff{uc:minf(x)}
can be obtained easily. This can be illustrated as follows.
When $\rank(Z^*)=1$, the constraint in \reff{sdp:momsos}
implies $Z^*=[x^*]_d[x^*]_d^T$ for some $x^*\in \re^n$,
and hence $y^*= -[x^*]_{2d}$. Then, for any $x\in \re^n$,
\[
-f(x^*) = -f_0 + b^T y^* \geq -f_0 + \sum_{\af \in \mathbb{U}_{2d}^n} -b_\af x^{\af}
= -f(x).
\]
The inequality above follows the optimality of $y^*$
and the fact that $Z=[x]_d[x]_d^T$ is always feasible for \reff{sdp:momsos}.
So $x^*$ is a global minimizer.


When $\rank(Z^*)>1$,
several global minimizers for \reff{uc:minf(x)}
could be obtained if a so called {\it flat extension condition (FEC)} holds.
We refer to Curto and Fialkow \cite{CF00} for FEC,
and Henrion and Lasserre \cite{HL05} for
how to extract global minimizers from $Z^*$.
Typically, FEC fails if there are a lot of global minimizers
or the SOS relaxation is not exact.



\subsection{Numerical experiments}
\label{subsec:uc-num}

Now we show some numerical examples of
applying Algorithm~\ref{alg:Newton-CG} to solve the SDP \reff{sdp:stdsos}-\reff{sdp:momsos}
which is equivalent to SOS relaxation~\reff{opt:stdsos}.
The computation is implemented on a Linux Desktop of 3.8GB memory and 2.8GHz CPU frequency.
Algorithm~\ref{alg:Newton-CG} is implemented in Matlab, and its parameters are:
$\eps^{in}=\eps^{out}=10^{-6}$, $\rho=5$, $\sig_0=1$, $\sig_{max} = 10^6$, $\dt = \half$, $K=500$.
The inverse $(\mc{A}\mc{A}^*)^{-1}$ via the formula \reff{uc-def:u} is used as a preconditioner.
Its inner loop is terminated if it does not converge within $25$ steps,
and the outer loop is terminated if it does not converge within $20$ steps.
If the optimal $Z^*$ satisfies FEC,
we extract one or several global minimizers $x^*$ by using the method in \cite{HL05};
otherwise, we just set $x^* = Z^*(2:n+1,1)$.
In either case, the relative error of $x^*$ is measured as
\[
err = \frac{|f(x^*)-f_{sos}^{uc}|}{\max\{1, |f(x^*)|\} }.
\]
Since $f_{sos}^{uc} \leq f_{min}^{uc} \leq f(x^*)$,
$err = 0$ if and only if $x^*$ is a global minimizer.
So the smaller $err$ is, the better the obtained solution is.



\begin{exm} Minimize the quartic polynomial $(x^2)^T A x^2 - x^TBx$.
Here $x^2$ denotes the component-wise square of $x$,
and $A,B$ are symmetric matrices given as
\[
A =
\bbm
n & n-1 & n-2 & \cdots  & 2& 1 \\
n-1 & n & 0 & \ddots  & 0 & 2 \\
n-2 & 0 & n & \ddots  & 0 & 3 \\
\vdots & \vdots & \ddots & \ddots & \vdots & \vdots \\
2 & 0 & \cdots & 0 & n &  n-1 \\
1 & 2 & \cdots & n-2 & n-1 & n
\ebm,
B =
\bbm
1 & 1 & 1 & \cdots & 1 & 1 \\
1 & 1 & 0 & \ddots & 0 & 1 \\
1 & 0 & 1 & \ddots & \vdots &\vdots \\
\vdots & \vdots & \ddots & \ddots & \vdots &\vdots \\
1 & 0 & \cdots & 0 & 1 & 1 \\
1 & 1 & \cdots & 1 &  1 & 1
\ebm .
\]
When $n=50$, the size of the resulting SDP \reff{sdp:stdsos}-\reff{sdp:momsos} is $(N,m)=(1326,316250)$.
Solving \reff{sdp:stdsos}-\reff{sdp:momsos} takes about $10$ hours.
The computed lower bound $f_{sos}^{uc}= -1.2919$ (only four decimal digits are shown).
The optimal $Z^*$ has rank two and satisfies FEC.
So we get two solutions from $Z^*$ (only four decimal digits are shown):
\[
\baray{l}
\pm(0.0000,  0.1504,  0.1480,  0.1456,  0.1432,
\ldots,  0.0584,  0.0572,  0.0560, 0.0548,  0.2188).
\earay
\]
\iffalse
%% the full solution
\[
\baray{l}
\pm(0.0000,  0.1504,  0.1480,  0.1456,  0.1432,  0.1407,  0.1383, 0.1358,  0.1334,  0.1309,\\
\quad\, 0.1285,  0.1260,  0.1235, 0.1211,  0.1186,  0.1162,  0.1138,  0.1114,  0.1090,  0.1066,\\
\quad\, 0.1043,  0.1019, 0.0997,  0.0974, 0.0951,  0.0930, 0.0908,  0.0887,  0.0866,  0.0845,\\
\quad\, 0.0826,  0.0806,  0.0787,  0.0769, 0.0751,  0.0733,  0.0716,  0.0699,  0.0683, 0.0668,\\
\quad\, 0.0653,  0.0638,  0.0624,  0.0610,  0.0597,  0.0584,  0.0572,  0.0560, 0.0548,  0.2188).
\earay
\]
\fi
The relative error is about $ 10^{-8}$.
So they are global minimizers up to a tiny numerical error.
% Minimize (x^2)^T*A x^2 - x^T*B*x
% A is sparse boundsary + diagonal
% B is sparse boundsary + diagonal
%
% n= 50  deg= 4
%
% it    time         gap        lg(rrd)  lg(rrp)     sigma
%  1  07:04:15   -3.24128e+05    2.454   -3.671         1
%  2  08:05:44   -3.38517e+00   -1.110   -6.726         5
%  3  08:49:31   -1.17853e+00   -1.565   -7.060        25
%  4  09:47:27   -4.80567e-01   -2.045   -6.314       125
%  5  10:26:48   -9.15869e-02   -2.811   -7.238       625
%  6  10:33:31   -4.67365e-03   -4.098   -9.105      3125
%  7  10:38:35   -5.17572e-05   -6.036   -9.246     15625
%total time:  38315.551
%
%
%ans =
%
%   -1.2919   -1.2919
%
%
%error =
%
%  -1.4527e-08
%
% there are totally   2 extracted solutions
%ans =
%
%  Columns 1 through 12
%
%    0.0000   -0.1504   -0.1480   -0.1456   -0.1432   -0.1407   -0.1383   -0.1358   -0.1334   -0.1309   -0.1285   -0.1260
%
%  Columns 13 through 24
%
%   -0.1235   -0.1211   -0.1186   -0.1162   -0.1138   -0.1114   -0.1090   -0.1066   -0.1043   -0.1019   -0.0997   -0.0974
%
%  Columns 25 through 36
%
%   -0.0951   -0.0930   -0.0908   -0.0887   -0.0866   -0.0845   -0.0826   -0.0806   -0.0787   -0.0769   -0.0751   -0.0733
%
%  Columns 37 through 48
%
%   -0.0716   -0.0699   -0.0683   -0.0668   -0.0653   -0.0638   -0.0624   -0.0610   -0.0597   -0.0584   -0.0572   -0.0560
%
%  Columns 49 through 50
%
%   -0.0548   -0.2188
%
%
%ans =
%
%  Columns 1 through 12
%
%   -0.0000    0.1504    0.1480    0.1456    0.1432    0.1407    0.1383    0.1358    0.1334    0.1309    0.1285    0.1260
%
%  Columns 13 through 24
%
%    0.1235    0.1211    0.1186    0.1162    0.1138    0.1114    0.1090    0.1066    0.1043    0.1019    0.0997    0.0974
%
%  Columns 25 through 36
%
%    0.0951    0.0930    0.0908    0.0887    0.0866    0.0845    0.0826    0.0806    0.0787    0.0769    0.0751    0.0733
%
%  Columns 37 through 48
%
%    0.0716    0.0699    0.0683    0.0668    0.0653    0.0638    0.0624    0.0610    0.0597    0.0584    0.0572    0.0560
%
%  Columns 49 through 50
%
%    0.0548    0.2188
\end{exm}




\begin{exm}
Solve the quartic polynomial optimization
\[
\min_{x\in\re^n} \quad (x_1^2+\cdots+x_n^2)^2 +
\sum_{ \substack{ 1\leq i < j, j+1 < k \leq n} } \frac{1}{n}  x_ix_jx_k.
\]
When $n=50$, the resulting SDP \reff{sdp:stdsos}-\reff{sdp:momsos}
has size $(N,m)=(1326,316250)$. Solving it takes about $1.5$ hours.
We get $f_{sos}^{uc}=-0.1246$ (only four decimal digits are shown) and
$\rank(Z^*)=1$. The solution extracted (only four decimal digits are shown) is
\[
-(0.1127,  0.1127,  0.1126,  0.1125,  0.1124,
\ldots,  0.1086,  0.1085,  0.1085,  0.1082,  0.1128).
\]
\iffalse
%% the following is the full solution
\[
\baray{l}
-(0.1127,  0.1127,  0.1126,  0.1125,  0.1124,  0.1123,  0.1122,  0.1121,  0.1120,  0.1119, \\
\quad\, 0.1118,  0.1117,  0.1116,  0.1115,  0.1114,  0.1113,  0.1112,  0.1111,  0.1111,  0.1110, \\
\quad\, 0.1109,  0.1108,  0.1107,  0.1106,  0.1105,  0.1104,  0.1103,  0.1102,  0.1101,  0.1100, \\
\quad\, 0.1099,  0.1099,  0.1098,  0.1097,  0.1096,  0.1095,  0.1094,  0.1093,  0.1092,  0.1091, \\
\quad\, 0.1091,  0.1090,  0.1089,  0.1088,  0.1087,  0.1086,  0.1085,  0.1085,  0.1082,  0.1128).
\earay
\]
\fi
Its relative error is about $2 \cdot 10^{-7}$.
So it is a global minimizer up to a small numerical error.
%
% Minimize sum_{i<j, j+1<k} x_i*x_j*x_k/n + (x^Tx)^2
%
% n= 50  deg= 4
%
% it    time         gap        lg(rrd)  lg(rrp)     sigma
%  1  00:50:21   -1.68852e+03    1.313   -7.860         1
%  2  01:13:23   -5.05915e-01   -0.920   -6.767         5
%  3  01:20:15   -2.12447e-01   -1.547   -9.454        25
%  4  01:27:19   -3.17376e-02   -2.465   -7.907       125
%  5  01:31:23   -1.21901e-03   -3.896  -10.396       625
%  6  01:33:51   -9.68787e-06   -5.991   -8.182      3125
%  7  01:43:01    2.03762e-07   -8.782   -7.854     15625
%total time:   6181.001
%
%
%ans =
%
%   -0.1246   -0.1246
%
%
%error =
%
%  -2.1963e-07
%
% there are totally   1 extracted solutions
%ans =
%
%  Columns 1 through 12
%
%   -0.1127   -0.1127   -0.1126   -0.1125   -0.1124   -0.1123   -0.1122   -0.1121   -0.1120   -0.1119   -0.1118   -0.1117
%
%  Columns 13 through 24
%
%   -0.1116   -0.1115   -0.1114   -0.1113   -0.1112   -0.1111   -0.1111   -0.1110   -0.1109   -0.1108   -0.1107   -0.1106
%
%  Columns 25 through 36
%
%   -0.1105   -0.1104   -0.1103   -0.1102   -0.1101   -0.1100   -0.1099   -0.1099   -0.1098   -0.1097   -0.1096   -0.1095
%
%  Columns 37 through 48
%
%   -0.1094   -0.1093   -0.1092   -0.1091   -0.1091   -0.1090   -0.1089   -0.1088   -0.1087   -0.1086   -0.1085   -0.1085
%
%  Columns 49 through 50
%
%   -0.1082   -0.1128
\end{exm}




\begin{exm}
Minimize the following least squares polynomial
\[
\left(\sum_{i=1}^n x_i-1\right)^2 +
\left(\sum_{i=1}^n x_i^2 - 1 \right)^2 +
\left(\sum_{i=1}^n x_i^3 - 1\right)^2 +
\sum_{i=1}^n \left(x_{i-1}^2+x_i^2+x_{i+1}^2 - x_i^3-1\right)^2
\]
where $x_0=x_{n+1}=0$.
For $n=16$, the resulting SDP \reff{sdp:stdsos}-\reff{sdp:momsos} has size $(N,m)=(969, 74612)$.
Solving \reff{sdp:stdsos}-\reff{sdp:momsos} takes about $12$ hours.
The computed lower bound $f_{sos}^{uc} = 7.5586$ (only four decimal digits are shown).
The optimal $Z^*$ has rank two and FEC holds.
Two extracted solutions (only four decimal digits are shown) are
\[
\baray{l}
(0.0039,   0.6285,   0.5370,   0.0259,   -0.4324,   -0.4266,   0.1540,   -0.5108,   0.2172,   -0.4029 \\
\quad 0.4400,  -0.4307,   0.0230,   0.5378,   0.6285,   0.0039), \\
(0.0039,   0.6285,   0.5378,   0.0230,   -0.4307,   0.4400,   -0.4029,   0.2172   -0.5108,   0.1540  \\
\quad -0.4266,  -0.4324,   0.0259,   0.5370,   0.6285,    0.0039).
\earay
\]
The relative error is around $10^{-8}$.
So they are global minimizers up to a tiny numerical error.
% Minimize  (1*x-1)^2 + (1*x^2-1)^2 + (1*x^3-1)^2 + sum_i (x_{i-1}^2+x_i^2+x_{i+1}^2-x_i^3-1)^2
%
% n= 16  deg= 6
%
% it    time         gap        lg(rrd)  lg(rrp)     sigma
%  1  01:44:00   -3.47018e+01    0.469   -4.496         1
%  2  02:37:25   -1.50975e-01   -1.187   -8.609         5
%  3  04:08:50   -1.97183e-02   -2.070   -7.947        25
%  4  06:58:11   -5.50590e-05   -4.265   -4.656       125
%  5  09:46:38   -2.00290e-07   -6.632   -5.465       625
%  6  12:34:58    1.29022e-08   -8.800   -6.052      3125
%total time:  45298.937
%
%
%z =
%
%    0.0039
%    0.6285
%    0.5370
%    0.0259
%   -0.4324
%   -0.4266
%    0.1540
%   -0.5108
%    0.2172
%   -0.4029
%    0.4400
%   -0.4307
%    0.0230
%    0.5378
%    0.6285
%    0.0039
%
%
%ans =
%
%    7.5586    7.5586
%
%
%error =
%
%  -1.3922e-08
%
% there are totally   2 extracted solutions
%ans =
%
%  Columns 1 through 12
%
%    0.0039    0.6285    0.5370    0.0259   -0.4324   -0.4266    0.1540   -0.5108    0.2172   -0.4029    0.4400   -0.4307
%
%  Columns 13 through 16
%
%    0.0230    0.5378    0.6285    0.0039
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0039    0.6285    0.5378    0.0230   -0.4307    0.4400   -0.4029    0.2172   -0.5108    0.1540   -0.4266   -0.4324
%
%  Columns 13 through 16
%
%    0.0259    0.5370    0.6285    0.0039
%
\end{exm}




\begin{exm}
Solve the least square problem
\[
\min_{x\in\re^n} \quad \sum_{i=1}^n
\left(1-\sum_{j=1}^i (x_j+x_j^2)\right)^2 +
\left(1-\sum_{j=i}^n (x_j+x_j^3)\right)^2.
\]
For $n=15$, the resulting SDP \reff{sdp:stdsos}-\reff{sdp:momsos} has size $(N,m)= (816,54263)$.
Solving \reff{sdp:stdsos}-\reff{sdp:momsos} takes about $4$ hours.
The computed lower bound $f_{sos}^{uc} = 2.1011$ (only four decimal digits are shown).
Since $\rank(Z^*)=1$, we get a solution (only four decimal digits are shown)
\[
\baray{l}
(\, 0.5233,    0.2429,    0.0802,    0.0011,   -0.0003,    0.0001,   -0.0000,    0.0000,  -0.0000, \\
\quad   0.0001,   -0.0005,    0.0020,   -0.0088,    0.0361,    0.6011).
\earay
\]
Its relative error is around $10^{-9}$.
So it is a global minimizer up to a tiny numerical error.
%
% Minimize sum_{i=1}^n
%  (x_1+...+x_i + x_1^2+...+x_i^2 - 1)^2 +  (x_i+...+x_n + x_i^3+...+x_n^3 - 1)^2
% n= 15  deg= 6
%
% it    time         gap        lg(rrd)  lg(rrp)     sigma
%  1  01:16:43   -4.41515e-01   -0.479   -3.850         1
%  2  02:33:53   -3.70608e-04   -2.920   -7.441         5
%  3  03:11:25   -3.66005e-06   -4.747   -7.163        25
%  4  04:12:22   -1.35241e-08   -7.126   -7.002       125
%total time:  15142.133
%
%
%ans =
%
%    2.1011    2.1011
%
%
%error =
%
%  -1.4866e-09
%
% there are totally   1 extracted solutions
%ans =
%
%  Columns 1 through 12
%
%    0.5233    0.2429    0.0802    0.0011   -0.0003    0.0001   -0.0000    0.0000   -0.0000    0.0001   -0.0005    0.0020
%
%  Columns 13 through 15
%
%   -0.0088    0.0361    0.6011
%
\end{exm}





%
%  sensor network is deleted
%

\iffalse

\begin{exm}[Sensor Network Localization]
This problem is also known as {\it graph realization} or {\it distance geometry}
in practical applications.
Given a graph $G=(V,E)$ along with a positive number associated with each edge,
the problem is to the coordinates of every vertex so that
the Euclidean distance between any two adjacent vertices
is equal to the positive number associated with that edge connecting them.
This problem has wide applications in Ad Hoc wireless communication,
protein folding, information theory, etc.


The sensor network localization problem can be mathematically formulated as follows:
find a sequence of unknown vectors ({\it sensors})
$u_1,u_2,\ldots,u_s \in \re^k(k=1,2,\ldots)$
such that distances between these sensors and some other fixed vectors
({\it anchors}) $a_1,\ldots,a_\ell$
are equal to some given constants.
Given edge subsets
\[
\mc{E}_{S}\subset \{(i,j):\, 1 \leq i<j \leq s \}, \quad
\mc{E}_{A}=\{(i,j): 1\leq i \leq s, 1\leq j \leq \ell \},
\]
for every $(i,j) \in \mc{E}_{S}$,
let $d_{ij}$ be the given distance between $u_i$ and $u_j$,
and for every $(i,j) \in \mc{E}_{A}$,
let $e_{ij}$ be the given distance between $u_i$ and $a_j$.
Denote $u_i = (x_{ki-k+1},\ldots, x_{ki})$ for $i=1,\ldots,s$.
Then the sensor network localization problem is equivalent
to finding coordinates $x_1,\ldots,x_{ks}$ such that
\begin{align*}
\|u_i - u_j \|_2^2 =
\sum_{t=0}^{k-1} (x_{ki-t}-x_{kj-t})^2 &= d_{ij}^2, \forall\, (i,j) \in\mc{E}_{S},  \\
\|u_i - a_j \|_2^2 =
\sum_{t=0}^{k-1} (x_{ki-t}-a_{kj-t})^2 &= e_{ij}^2, \forall\, (i,j) \in\mc{E}_{A}.
\end{align*}
Usually solving this polynomial system directly is very expensive.
However, the problem is equivalent to finding the
global minimizers of the quartic polynomial
\begin{align} \label{L2err}
f(x):=
\sum_{(i,j)\in \mc{E}_{S}}
\left( \| u_i - u_j \|_2^2 - d_{ij}^2\right)^2 +\,
\sum_{(i,j)\in \mc{E}_{A}}
\left( \| u_i - a_j \|_2^2 - e_{ij}^2 \right)^2.
\end{align}
Thus the SOS methods can be applied to minimize $f(x)$.
The exist other methods for solving the sensor network localization problem,
like semidefinite programming methods, second-order cone methods,
sparse SOS methods.
These methods can applied to solve large scale problems, but might fail to
find the correct sensor locations, especially the sensors are not uniquely
localization in a higher dimensional lifting spaces.
In such situations, minimizing $f(x)$ directly would give better solutions.




We randomly generate test problems
which are similar to those given in \cite{BLTWY05}.
First, we randomly generate $n=500$ sensor locations
$x_1^*,\ldots,x_n^*$ from the unit square
$[-0.5,\,0.5]\times [-0.5,\,0.5]$.
The anchors $\{a_1,a_2,a_3,a_4\}\, (m=4)$ are chosen to be
four fixed points $(\pm 0.45,\, \pm 0.45)$.
Choose edge set $\mathcal{A}$ such that
for every sensor $x_i^*$ there are at most $10$ sensors
$x_j^*\,(j>i)$ with $(i,j)\in \mathcal{A}$ and
$\|x_i^*-x_j^*\|_2 \leq 0.3$.
For every $(i,j)\in \mathcal{A}$,
compute the distance $\|x_i^*-x_j^*\|_2 = d_{ij}$.
Choose edge $\mathcal{B}$ such that every anchor
is connected to all the sensors within distance $0.3$.
For every $(i,k)\in \mathcal{B}$,
compute the distance $\|x_i^*-a_k\|_2 = e_{ik}$.
Then we apply sparse SOS relaxation \reff{sumsos1}--\reff{summom3}
to minimize polynomial function~\reff{L2err}.
The accuracy of computed sensor locations
$\hat x_1,\ldots, \hat x_n$
will\vspace*{-3pt} be measured by the Root Mean Square Distance (RMSD),
which is defined as
$
\mbox{RMSD} = \left(\frac{1}{n}
\sum_{i=1}^n\|\hat x_i -x_i^* \|_2^2 \right)^{\frac{1}{2}}.
$
We use {\it SeDuMi} to solve the sparse SOS
relaxation~\reff{sumsos1}--\reff{summom3}
on a Linux machine with 1.46 GHz CPU and 0.98GB memory.
The problem can be solved within about $18$ CPU minutes
with accuracy $\mathcal{O}(10^{-6})$.
%The maximum, median and minimum of the error RMSD
%are $0.0272,4.7\cdot 10^{-6},2.3\cdot 10^{-7}$
%respectively. For those problems with $RMSD$ greater than
%$10^{-4}$, only a few moment matrices satisfy
%the flat extension condition, and all the others have rank one.
%This is probably because that these few sensors can not be localized
%and the others can be localized.


The sparse SOS method based on chordal extension like in \cite{WKKM}
is usually not practical for solving sensor network localization problems,
because the corresponding csp matrices usually do not have sparse
Cholesky factorizations and the resulting chordal extensions
are often too dense to be useful.
We refer to \cite{njw_senloc} for more details
about sparse SOS methods for sensor network localization.
\qed
\end{exm}

\fi




\begin{exm}[Random polynomials]
We test the performance of Algorithm \ref{alg:Newton-CG}
in solving SOS relaxations for minimizing randomly generated polynomials.
To insure the polynomials have finite global minimum,
we generate $f(x)$ randomly as follows
\[
f(x) = f^T[x]_{2d-1} + [x^d]^T F^TF [x^d],
\]
where $f \in \re^{\binom{n+2d-1}{2d-1}}$ is a random vector of Gaussian distribution
and $F \in \re^{\binom{n+d-1}{d} \times \binom{n+d-1}{d}}$
is a random matrix of Gaussian distribution.
Here $[x^d]$ denotes the vector of monomials of degree equaling $d$.
The generated $f(x)$ almost always has a finite global minimum.
%This class of random polynomials is more general than
%the class of testing polynomials used by Parrilo and Sturmfels \cite{PS01}.
The computational results are  listed in Tables~\ref{UCtab:rand:deg=4}-\ref{UCtab:rand:deg=10}.
The first column is the number of variables;
the second column is the size of the resulting SDP \reff{sdp:stdsos}-\reff{sdp:momsos};
the third column is the number of instances generated;
the fourth column is the maximum, median, and minimum time consumed by the computer;
and the last column is the maximum, median, and minimum
relative errors of obtained solutions.


\begin{table}[htb]
\centering
\begin{tabular}{|r||c||c||r|r|r||c|} \hline
$n$ & (N,m) & \# inst & \multicolumn{3}{c||}{time (min, med, max)} & err (min, med, max) \\ \hline
20 &  (231, 10625) & 20 &   0:00:45  &   0:01:07   &   0:01:59  & (5e-9, 7e-8,  2e-8)\\ \hline
30 &  (496, 46375) & 20 & 0:10:45 & 0:11:44  & 0:12:43 & (6e-9, 5e-8, 1e-8) \\ \hline
40 &  (861, 135750) & 10 & 0:28:42 & 1:02:36 & 1:22:24 & (3e-9, 3e-8, 9e-8)\\ \hline
50 & (1326, 316250) & 5 & 2:04:43 & 3:33:31 & 4:15:38 & (6e-8, 4e-7, 9e-7) \\ \hline
60 & (1891, 635375) & 5 & 5:17:42&  8:03:04& 15:23:23& (7e-8, 1e-7, 5e-7) \\ \hline
70 & (2556, 1150625) & 5 & 10:25:07 & 13:06:38 & 20:44:31 & (2e-9, 7e-8, 3e-8) \\ \hline
80 & (3321, 1929500) & 3 & 30:09:34 & 32:17:43 & 33:30:12  & (6e-8, 8e-8, 1e-7) \\ \hline
90 & (4186, 3049500) & 1 & 63:55:41 & 63:55:41 & 63:55:41 & (4e-7, 4e-7, 4e-7)\\ \hline
100& (5151, 4598125) & 1 & 199:49:33 & 199:49:33 & 199:49:33 & (6e-7, 6e-7, 6e-7)\\ \hline
\end{tabular}
\caption{Computational results for random \reff{uc:minf(x)} of degree $4$}
\label{UCtab:rand:deg=4}

\end{table}



\begin{table}[htb]
\centering
\begin{tabular}{|r||c||c||r|r|r||c|} \hline
$n$ & (N,m) & \# inst & \multicolumn{3}{c||}{time (min, med, max)} & err (min, med, max) \\ \hline
10 &  (286, 8007) & 20 & 0:01:45  &  0:00:54  & 0:03:16     &
(9e-10, 4e-8, 6e-7) \\ \hline
15 &  (816, 54263) & 10  &  0:39:05 &  1:07:42    &   1:33:51   &
 (3e-8,  6e-7, 2e-6) \\ \hline
20 &  (1771, 230229)  & 3  &  11:47:23    &  21:35:17    &   25:36:08   &
  (9e-8, 7e-8,  1e-6)  \\ \hline
25 &  (3276, 736280)  & 1  &  104:41:36 &  104:41:36 & 104:41:36    &  (6,e-7, 6e-7, 6e-7)  \\ \hline
\end{tabular}
\caption{Computational results for random \reff{uc:minf(x)} of degree $6$}
\label{UCtab:rand:deg=6}

\end{table}



\begin{table}[htb]
\centering
\begin{tabular}{|r||c||c||r|r|r||r|} \hline
$n$ & (N,m) & \# inst & \multicolumn{3}{c||}{time (min, med, max)} & err (min, med, max) \\ \hline
8 &   (495, 12869) & 20  &  0:03:44 & 0:07:37&  0:19:42   &
  (1e-8, 3e-7, 8e-6)\\ \hline
10 &  (1001, 43757)  & 10  &  1:14:48 & 1:40:58  &  2:20:22   &
  (9e-9, 3e-8, 4e-6) \\ \hline
12 &  (1820, 125969)  & 3  & 13:04:06 &  13:11:35 & 13:34:40 &
 (8e-8,  1e-7,  6e-7) \\ \hline
15 &  (3876, 490313)  & 1  &  217:46:12   &   217:46:12   &   217:46:12   &
 (9e-8, 9e-8, 9e-8) \\ \hline
\end{tabular}
\caption{Computational results for random \reff{uc:minf(x)} of degree $8$}
\label{UCtab:rand:deg=8}

\end{table}




\begin{table}[htb]
\centering
\begin{tabular}{|r||c||c||r|r|r||r|} \hline
$n$ & (N,m) & \# inst & \multicolumn{3}{c||}{time (min, med, max)} & err (min, med, max) \\ \hline
6 &  (462, 8007) & 20  &  0:03:18 & 0:05:19  &  0:06:22    &
  (2e-8, 1e-7, 2e-7) \\ \hline
8 &  (1287, 43757)  & 10  &  1:51:38 & 3:47:52    & 6:53:22 &
  (5e-8,  3e-7, 4e-6) \\ \hline
9 &  (2002, 92377)  & 3 & 9:18:17 &  12:09:07  & 20:52:36 &
 (2e-8, 2e-7, 3e-6) \\ \hline
10 & (3003, 184755)  & 1  &  72:40:36 &  72:40:36  & 72:40:36  &
 (9e-6, 9e-6, 9e-6) \\ \hline
\end{tabular}
\caption{Computational results for random \reff{uc:minf(x)}of degree $10$}
\label{UCtab:rand:deg=10}

\end{table}



When $f(x)$ has degree $4$ ($d=2$),
the computational results are in Table~\ref{UCtab:rand:deg=4}.
As we can see, SOS relaxation \reff{opt:stdsos} is solved quite well.
For $n=20\sim 30$, solving \reff{sdp:stdsos}-\reff{sdp:momsos} takes a couple of minutes;
for $n=40\sim 50$, it takes a couple of hours;
for $n=60\sim 80$, it takes about five to thirty hours;
for $n=90$, it takes about sixty hours;
for $n=100$, it takes about one week.
%All the computed optimal $Z^*s$ have rank one.
%Good solutions with tiny errors are found.



When $f(x)$ has degree $6$ ($d=3$),
the computational results are in Table~\ref{UCtab:rand:deg=6}.
For $n=15$, solving \reff{sdp:stdsos}-\reff{sdp:momsos} takes about forty to ninety minutes;
for $n=20$, it takes about ten to twenty five hours;
for $n=25$, it takes about one hundred hours.
%All the computed optimal $Z^*s$ have rank one
%and the obtained solutions have small relative errors.





When $f(x)$ has degree $8$ ($d=4$),
the computational results are in Table~\ref{UCtab:rand:deg=8}.
For $n=10$, solving \reff{sdp:stdsos}-\reff{sdp:momsos} takes about one or two hours;
for $n=12$, it takes about thirteen hours;
for $n=25$, it takes about two hundred hours.
%All the computed optimal $Z^*s$ have rank one.
%We easily get accurate solutions from $Z^*$.




When $f(x)$ has degree $10$ ($d=5$),
the computational results are in Table~\ref{UCtab:rand:deg=10}.
For $n=6$, solving \reff{sdp:stdsos}-\reff{sdp:momsos} takes a couple of minutes;
for $n=8$, it takes about a couple of hours;
for $n=9$, it takes about ten to twenty hours;
for $n=10$, it takes about seventy hours.
%All the computed optimal $Z^*s$ have rank one,
%and the obtained solutions have tiny errors.



For polynomials generated in this example,
SOS relaxation \reff{opt:stdsos} is very successful in solving \reff{uc:minf(x)} globally.
For all the generated instances, the optimal $Z^*s$ have rank one,
and highly accurate global optimal solutions with tiny errors are found.
This is consistent with the computational results shown in \cite{PS01}.
The computational results demonstrate that
Algorithm \ref{alg:Newton-CG} is very efficient in
solving SOS relaxations for large scale polynomial optimization.
Significantly bigger problems would be solved very well.
A quartic polynomial optimization with $100$ variables would be solved on a regular computer.
This is almost impossible by applying prior existing SOS solvers
using interior point type methods.

\end{exm}


\iffalse

%% the following is about using sparisty via Newton's polytope.
%% it is droped due to using sparsity in Section 5
%%

\subsection{Exploiting sparsity via Newton polytope}


In unconstrained polynomial optimization, it is typical that the polynomials to be minimized
have sparsity structures. To improve the computational efficiency,
it is important to exploit these features.


Let $f(x)$ be a polynomial and $New(f)$ be its Newton polytope.
Assume $ 0\in New(f)$.
If $f(x)-\gamma = \sum_i f_i(x)^2$,
then it can be shown
\[
New(f_i) \subseteq  \mbox{conv} \left(\half New(f_{even}) \right).
\]
Here $f_{even}$ is the polynomial consisting of monomials of $f(x)$
whose exponents are even integer vectors.
For convenience, denote
\[
\mc{G} = \mbox{lattice points of conv} \left(\half New(f_{even})\right),
\quad \mc{F} = \mc{G} + \mc{G}-\{0\}.
\]
Then $f(x)-\gamma$ is SOS if and only if
there exists $X\succeq 0$ such that
\[
f(x) - \gamma = [x]_{\mc{G}}^T X [x]_{\mc{G}}.
\]
Define constant matrices $C$ and $A_\af$ as
\[
[x]_{\mc{G}}  [x]_{\mc{G}}^T =
C + \sum_{\af \in \mc{F} } A_\af x^\af,
\]
then $f(x)-\gamma$ being SOS is equivalent to the existence of $X$ satisfying
\[
\baray{rcl}
C \bullet X + \gamma  & = & f_0 \\
A_\af \bullet X & = & f_\af, \quad \forall \, \af\in\mc{F}\\
X & \succeq & 0.
\earay
\]
Therefore, the SOS relaxation \reff{opt:stdsos} is equivalent to the SDP
\[
\baray{rl}
\min &  C \bullet X \\
s.t. & A_\af \bullet X = f_\af, \quad \forall \, \af\in\mc{F}\\
& X \succeq  0.
\earay
\]

\begin{exm}[Bilinear least squares (BLS) problem]
Suppose $(x^{I})^TF_kx^{II}$ ($k=1,\ldots,\ell$)
are bi-linear forms in
$(x^{I},x^{II})\in \re^{n_1}\times \re^{n_2}$.
Here $F_k \in \re^{n_1 \times n_2}$ is a real matrix.
Given constants $h_1,\ldots,h_\ell$,
we are interested in solving the bilinear equations
\be
(x^{I})^TF_kx^{II} = h_k, \quad k=1,\ldots,\ell.
\ee
If there are no solutions, then we are interested in
bilinear least squares (BLS) solutions.
Denote $x=(x^{I},x^{II})$, and let $n=n_1+n_2$.
The BLS problem is equivalent to finding the global
minimizers of the polynomial
\[
f(x) = \sum_{k=1}^\ell \Big((x^{I})^TF_kx^{II} - h_k \Big)^2.
\]
For this polynomial, we have the support sets
\[
\baray{rcl}
\mc{G} &=& \{0\} \cup \{ e_i+e_j: 1\leq i \leq n_1, n_1+1 \leq j \leq n \},\\
\mc{F} &=& \big(\mc{G}\backslash\{0\}\big) \cup \{ e_{i_1}+e_{i_2}+e_{j_1}+e_{j_2}:
1\leq i_1,i_2 \leq n_1, n_1+1 \leq j_1,j_2 \leq n \}.
\earay
\]
Here $e_i$ denotes the $i$-th unit vector in the space $\re^n$.

\end{exm}


\fi







\section{Homogeneous polynomial optimization}
\setcounter{equation}{0}
\label{sec:hmgopt}


A frequently encountered polynomial optimization is
\be \label{hmgpop:sphere}
\baray{rl}
f_{min}^{hmg}:=
\underset{x\in\re^n}{\min} & f(x) \\
s.t. & \|x\|_2 = 1
\earay
\ee
where $f(x)$ is a form (homogeneous polynomial).
Assume its degree $\deg(f)=2d$ is even, i.e., $f(x)$ is an even form.
The case that $\deg(f)$ is odd can be transformed to
minimizing an even form over the unit sphere
$\mathbb{S}^{n-1}=\{x\in\re^n: \|x\|_2 =1 \}$, as will be shown later.
In \reff{hmgpop:sphere}, $f_{min}^{hmg}$ denotes the minimum objective value.
Problem~\reff{hmgpop:sphere} is NP-hard when $d>1$.
%because checking co-positivity of a symmetric matrix
%is equivalent to minimizing a quartic form for over the unit sphere.



The standard SOS relaxation for \reff{hmgpop:sphere} is
\be  \label{hmgpop:sos}
\baray{rl}
f_{sos}^{hmg}:=
\max & \gamma \\
s.t. & f(x) - \gamma \cdot (x^Tx)^d \quad \mbox{ is SOS}.
\earay
\ee
Let $[x^d]$ be the vector of monomials having degree $d$
ordered lexicographically, i.e.,
\[
[x^d]^T = \bbm x_1^d & x_1^{d-1}x_2 & x_1^{d-1}x_3 & \cdots & x_{n-1}x_n^{d-1} & x_n^d \ebm.
\]
Denote $\N_{d} = \{\af\in \N^n: |\af| = d\}$. For each $\af \in \N_d$,
define $D_\af = \frac{(|\af|)!}{\af_1! \cdots \af_n!}$.
Let $D=\diag(D_\af)$ be a diagonal matrix.
It holds the relation
\be  \label{def:Daf}
(x^Tx)^d = \sum_{\af \in \N_{d} } D_\af  x_1^{2\af_1}\cdots x_n^{2\af_n}
= [x^d]^TD[x^d] = ([x^d][x^d]^T)\bullet D.
\ee
Define $0/1$ matrices $H_\af$ such that
\be  \label{def:H_af}
[x^d][x^d]^T = \sum_{\af \in \N_{2d} }  H_\af x_1^{\af_1}\cdots x_n^{\af_n}.
\ee
Then $f(x)-\gamma \cdot (x^Tx)^d$ being SOS
is equivalent to the existence of $X$ satisfying
\[
\baray{rl}
f_\af - \gamma D_{\af/2} &= H_\af \bullet X \quad \forall\, \af \in 2\N_{d}, \\
f_\af  & = H_\af \bullet X \quad \forall \, \af \in \N_{2d}\backslash 2\N_{d},\\
X& \succeq 0.
\earay
\]
Letting $\af = 2de_1$ in the above, we get
\[
\gamma = f_{2de_1} -  H_{2de_1} \bullet X.
\]
Denote $\mathbb{H}_{2d}^n = \N_{2d}\backslash \{2de_1\}$,
and set $r=(r_\af)_{\af\in \mathbb{H}_{2d}^n}$ as
\be \label{hmg-def:r}
r_\af =
\bca
D_{\af/2}  & \text{ if } \af \in 2\N_{d}\backslash\{2de_1\}, \\
0 & \text{ if } \af \in \N_{2d}\backslash\{2\N_{d}\}.
\eca
\ee
Define matrices $C, A_\af$ and scalars $b_\af$ as
\be \label{matdef:Homg}
C= H_{2de_1}, \quad
A_\af = H_\af - r_\af H_{2de_1}, \quad
b_\af = f_\af - r_\af f_{2de_1}  \quad \forall \, \af\in \mathbb{H}_{2d}^n.
\ee
Let $b=(b_\af)_{\af\in \mathbb{H}_{2d}^n}$.
Define linear operators
$\mc{H}, \mc{A}: \mc{S}^{|\N_{d}|} \rightarrow \re^{|\mathbb{H}_{2d}^n|}$
and a cone $\mc{K}$ as
\[
\mc{H}(X) = (H_\af \bullet X)_{\af\in \mathbb{H}_{2d}^n}, \quad
\mc{A}(X) = (A_\af \bullet X)_{\af\in \mathbb{H}_{2d}^n}, \quad
\mc{K} = \mc{S}_+^{|\N_{d}|}.
\]
Then SOS relaxation \reff{hmgpop:sos} is equivalent to
\be  \label{hmgsdp:prime}
\baray{rl}
f_{sdp}^{hmg}:=
\min & C \bullet X \\
s.t. & \mc{A}(X) = b, \, X \in \mc{K}.
\earay
\ee
The dual of \reff{hmgsdp:prime} is
\be  \label{hmgsdp:dual}
\baray{rl}
\max & b^Ty \\
s.t. &  \mc{A}^*(y) +Z = C, \, Z \in \mc{K}.
\earay
\ee
Here $\mc{A}^*(y) = \sum_{\af\in \mathbb{H}_{2d}^n} y_\af A_\af$ is the adjoint of $\mc{A}$.



\begin{theorem} \label{thm-bd:hmg}
For notations defined above, it holds
$f_{sos}^{hmg} = -f_{sdp}^{hmg}+f_{2de_1} \leq f_{min}^{hmg}$.
If SOS relaxation \reff{opt:stdsos} is exact,
then $f_{min}^{hmg} = -f_{sdp}^{hmg}+f_{2de_1}$.
\end{theorem}
\noindent
Its proof is the same as for Theorem~\ref{thm-bd:uc}, and hence is omitted.
A general error bound for \reff{hmgpop:sos} approximating \reff{hmgpop:sphere}
is given in \cite{Nie09b}.
\iffalse
% the proof is dropped.
\begin{proof}
Optimization \reff{hmgpop:sos} and \reff{hmgsdp:prime}
are equivalent. Their objectives satisfy the relation
\[
C \bullet X = -\gamma + f_{2de_1}.
\]
Maximizing $\gamma$ in \reff{hmgpop:sos} is equivalent to
minimizing $C\bullet X$ in \reff{hmgsdp:prime}.
The feasible sets of \reff{hmgpop:sos} and \reff{hmgsdp:prime}
are the same. So $f_{sos}^{hmg} = -f_{sdp}^{hmg}+f_{2de_1}$.
The relation $f_{sos}^{hmg} \leq f_{min}^{hmg}$ is obvious
since \reff{hmgpop:sos} is a relaxation for \reff{hmgpop:sphere}.


When SOS relaxation \reff{hmgpop:sos} is exact,
i.e., $f_{sos}^{hmg} = f_{min}^{hmg}$,
we immediately have
$f_{min}^{hmg} = -f_{sdp}^{hmg}+f_{2de_1}$.
\end{proof}
\fi


Now we turn to finding $(\mc{A}\mc{A}^*)^{-1}$.
Note the matrices $H_\af$ defined in \reff{def:H_af} are orthogonal to each other.
Thus $\mc{H}\mc{H}^*$ is diagonal,
%$
%\mc{H}\mc{H}^* = \diag( \mathbf{1}^T H_\af \mathbf{1} : \af \in \mathbb{H}_{2d}^n ).
%$
and
\[
(\mc{H}\mc{H}^*)^{-1} = \diag\big( (\mathbf{1}^T H_\af \mathbf{1})^{-1} :
\af \in \mathbb{H}_{2d}^n \big).
\]
Here $\mathbf{1}$ is the vector of all ones.
Thus, we have for all $\af,\bt \in  \mathbb{H}_{2d}^n$
\begin{align*}
A_\af \bullet A_\bt
& = (H_\af - r_\af H_{2de_1} ) \bullet (H_\bt - r_\bt H_{2de_1}) \\
& = H_\af \bullet H_\bt - H_{2de_1} \bullet ( r_\af H_\bt + r_\bt H_\af ) +
r_\af r_\bt H_{2de_1} \bullet H_{2de_1} \\
& = H_\af \bullet H_\bt + r_\af r_\bt,
\end{align*}
that is, $\mc{A}\mc{A}^* = \mc{H}\mc{H}^* + rr^T$.
By Sherman-Morrison formula, we have
\[
(\mc{A}\mc{A}^*)^{-1} = (\mc{H}\mc{H}^*)^{-1} -
\frac{(\mc{H}\mc{H}^*)^{-1}rr^T(\mc{H}\mc{H}^*)^{-1}}{1+r^T(\mc{H}\mc{H}^*)^{-1}r}.
\]
Setting vectors
\[
\mathbf{h} = (\,(\mathbf{1}^T H_\af \mathbf{1})^{-1} : \af \in \mathbb{H}_{2d}^n ), \quad
\mathbf{p} = \left(1+r^T\diag(\mathbf{h})r\right)^{-1/2}\diag(\mathbf{h})r,
\]
we get a neat formula
\be \label{hmg-def:h-p}
(\mc{A}\mc{A}^*)^{-1} = \diag(\mathbf{h}) - \mathbf{p}\mathbf{p}^T.
\ee




Algorithm~\ref{alg:Newton-CG} can be applied to solve \reff{hmgsdp:prime}-\reff{hmgsdp:dual}.
We use $(\mc{A}\mc{A}^*)^{-1}$ via the formula \reff{hmg-def:h-p} as a preconditioner.
Let $(X^*,y^*,Z^*)$ be optimal for \reff{hmgsdp:prime}-\reff{hmgsdp:dual}.
Then $-f_{sdp}^{hmg}+f_{2de_1}$ is a lower bound for the minimum $f_{min}^{hmg}$.
\iffalse
%% the following is deleted
and \reff{hmgsdp:dual} implies
\[
Z^* = H_{2de_1} - \sum_{\af \in \mathbb{H}_{2d}^n } y_\af^* (H_\af - r_\af H_{2de_1})
= \sum_{\af \in \N_{2d} } y_\af^* H_\af
\]
where $y_{2de_1}^* =1+\sum_{\af \in \mathbb{H}_{2d}^n } y_\af^* r_\af$.
So $Z^*$ is a matrix of homogeneous moments.

If $Z^*$ has rank one, then $f_{min}^{hmg}=f_{sos}^{hmg}$
and a global minimizer for \reff{hmgpop:sphere} can be obtained easily.
To see this, suppose $Z^*=uu^T$ has rank one.
Let $x^* = u/ \|u\|_2$.
The definition of $H_\af$ in \reff{def:H_af} implies $Z^*=[(x^*)^d][(x^*)^d]^T$.
Then, for any $x\in \mathbb{S}^{n-1}:=\{x\in \re^n: \|x\|_2 =1\}$,
\[
-f(x^*) = -f_{2de_1} + b^T y^* \geq -f_{2de_1} + \sum_{\af \in \mathbb{H}_{2d}^n} -b_\af x^{\af}
= -f(x),
\]
where the inequality is implied by the optimality of $y^*$
and the fact that $Z=[(x)^d][(x)^d]^T$ is always feasible for \reff{hmgsdp:dual}.
So $x^*$ is a global minimizer, and hence $f_{min}^{hmg}=f_{sos}^{hmg}$.
\fi
We would also get global minimizers from $Z^*$ when FEC holds.
Note that \reff{def:Daf} and \reff{matdef:Homg} imply
\be  \label{A.D=Z.D=0}
A_\af \bullet D = 0 \quad \forall \af \in \mathbb{H}_{2d}^n  \quad \text{ and } \quad
Z^* \bullet D = H_{2de_1} \bullet D = 1.
\ee
So $Z^*(de_i,de_i)$ ($Z$ is indexed by integer vectors in $\N^n$)
can not vanish for every $i$,
because otherwise we would get $Z^*=0$ contradicting \reff{A.D=Z.D=0}.
Up to a permutation of $x$, we can assume $Z^*(de_1,de_1)\ne 0$,
and normalize $Z^*$ as $\widehat{Z^*}=Z^*/Z^*(de_1,de_1)$.
Then $\widehat{Z^*}$ would be thought of as
a moment matrix of order $d$ in $n-1$ variables (see \cite{CF00}).
If $\widehat{Z^*}$ satisfies FEC, using the method in \cite{HL05},
we can get $v^{(1)},\ldots,v^{(r)}\in\re^{n-1}$ such that
\[
\widehat{Z^*} = \lmd_1 [v^{(1)}]_d[v^{(1)}]_d^T+\cdots +\lmd_r[v^{(r)}]_d[v^{(r)}]_d^T
\]
for some scalars $\lmd_i > 0$. Setting
$x^{(i)} = \big(1+\|v^{(i)}\|_2^2\big)^{-1/2}
\bbm 1 \\ v^{(i)} \ebm \in \mathbb{S}^{n-1}$, we get
\[
Z^*  = \nu_1 [(x^{(1)})^d] [(x^{(1)})^d]^T + \cdots + \nu_r [(x^{(r)})^d] [(x^{(r)})^d]^T
\]
for some scalars $\nu_i > 0$.
Then the relations \reff{def:Daf} and \reff{A.D=Z.D=0} imply $\nu_1+\cdots +\nu_r=1$.
Since every $Z^{(i)}=[(x^{(i)})^d] [(x^{(i)})^d]^T$ is feasible for \reff{hmgsdp:dual},
the optimality of $Z^*$ implies
every $x^{(i)}$ is a global minimizer for \reff{hmgpop:sphere}.


In particular, if $\rank(Z^*)=1$, then $\widehat{Z^*}$ satisfies FEC automatically.
When FEC fails, it is not clear how to extract a minimizer from $Z^*$.
This happens usually because there are a lot of global minimizers
or SOS relaxation \reff{hmgpop:sos} is not exact
($f_{sos}^{hmg} < f_{min}^{hmg}$).




\subsection{Numerical experiments}


Now we present some numerical examples of applying Algorithm~\ref{alg:Newton-CG}
to solve the SDP \reff{hmgsdp:prime}-\reff{hmgsdp:dual}
which is equivalent to SOS relaxation \reff{hmgpop:sos}.
As in Section~\ref{subsec:uc-num},
the computation is implemented in Matlab on the same machine.
The parameters and stopping criteria in Algorithm~\ref{alg:Newton-CG}
are also the same as there.
The inverse $(\mc{A}\mc{A}^*)^{-1}$ via the formula \reff{hmg-def:h-p} is used as a preconditioner.
%$\eps^{in}=\eps^{out}=10^{-6}$, $\rho=5$, $\sig_0=1$, $\sig_{max} = 10^6$, $\dt = \half$, $K=500$.
%We terminate the inner loop if it takes over $25$ steps,
%and terminate the outer loop if it takes more than $20$ steps.
If the optimal $Z^*$ satisfies FEC,
we apply the preceding procedure to extract one or several global minimizers $x^*$;
otherwise, we just choose $x^*$ to be the normalization of $Z^*(1:n,1)$.
The relative error of $x^*$ is measured as
\[
err = \frac{|f(x^*)-f_{sos}^{hmg}|}{\max\{1, |f(x^*)|\} }.
\]
Since $f_{sos}^{hmg} \leq f_{min}^{hmg} \leq f(x^*)$,
$err =0$ if and only if $x^*$ is a global minimizer of \reff{hmgpop:sphere}.




\begin{exm} Minimize the following square free quartic form over $\mathbb{S}^{n-1}$
\[
\sum_{ 1 \leq i < j < k < \ell \leq n}  (-i-j+k+\ell)x_ix_jx_kx_\ell .
\]
For $n=50$, the resulting SDP \reff{hmgsdp:prime}-\reff{hmgsdp:dual}
has size $(N,m)=(1275,292824)$. Solving \reff{hmgsdp:prime}-\reff{hmgsdp:dual} takes about $3.5$ hours.
The computed lower bound $f_{sos}^{hmg}=-140.4051$ (only four decimal digits are shown).
The optimal $Z^*$ has rank two and FEC holds.
There are 2 solutions (only four decimal digits are shown) extracted from $Z^*$:
\[
\baray{c}
(0.6393,  0.1217,  0.0827,  0.0575,  0.0401, \ldots,
-0.1653,  -0.1619,  -0.1581, -0.1541,  -0.1499), \\
(0.1499,  0.1541,   0.1581,   0.1619,   0.1653, \ldots,
-0.0401,  -0.0575,  -0.0827, -0.1217,  -0.6393).
\earay
\]
\iffalse
%% the following is the full solution
\[
\baray{l}
\,\,(0.6393,   0.1217,   0.0827,   0.0575,   0.0401,   0.0274,   0.0180,   0.0105,   0.0045,  -0.0006, -0.0051, \\
-0.0093, -0.0133,  -0.0172, -0.0211,  -0.0252,  -0.0295,  -0.0340,  -0.0387,  -0.0438, -0.0492, \\
-0.0549, -0.0609,  -0.0673, -0.0741,  -0.0811,  -0.0884,  -0.0958,  -0.1035,  -0.1112, -0.1188, \\
-0.1264,  -0.1337, -0.1407,  -0.1472,  -0.1531, -0.1584,  -0.1629,  -0.1664,  -0.1691, -0.1708, \\
-0.1715,  -0.1712,  -0.1700, -0.1680,  -0.1653,  -0.1619,  -0.1581, -0.1541,  -0.1499),
\earay
\]
\[
\baray{l}
(0.1499,  0.1541,   0.1581,   0.1619,   0.1653,   0.1680,   0.1700,   0.1712,   0.1715,   0.1708, 0.1691, \\
\,0.1664,  0.1629,   0.1584,   0.1531,   0.1472,   0.1407,   0.1337,   0.1264,   0.1188, 0.1112,   0.1035, \\
\,0.0958,  0.0884, 0.0811,   0.0741,   0.0673,   0.0609,   0.0549,   0.0492, 0.0438,   0.0387,   0.0340, \\
\,0.0295,  0.0252,   0.0211, 0.0172,   0.0133,   0.0093,   0.0051, 0.0006,  -0.0045,  -0.0105, \\
\, -0.0180, -0.0274,  -0.0401,  -0.0575,  -0.0827, -0.1217,  -0.6393).
\earay
\]
\fi
The error is around $3\cdot 10^{-9}$.
So they are global minimizers up to a tiny numerical error.
%
% Minimize the form:
%  sum_{i<j<k<l} (-i-j+k+el)*x_ix_jx_kx_l
%over the unit sphere
% n= 50  deg= 4
%
% it    time         gap        lg(rrd)  lg(rrp)   sigma
%  1  02:36:40   -1.37656e+04    1.768   -9.990        1
%  2  02:59:23   -2.70815e-02   -2.432   -8.833        5
%  3  03:07:46   -4.10529e-02   -2.590   -8.884       25
%  4  03:16:04   -3.01697e-02   -2.888   -8.894      125
%  5  03:24:20   -9.50329e-03   -3.498  -10.749      625
%  6  03:32:04   -7.47640e-04   -4.676   -9.807     3125
%  7  03:37:47   -1.11284e-05   -6.521   -8.718    15625
%total time:  13067.212
%
%ans =
%
% -140.4051 -140.4051
%
%
%error =
%
%  -3.8854e-09
%
%there are totally   2 extracted solutions
%ans =
%
%  Columns 1 through 12
%
%    0.6393    0.1217    0.0827    0.0575    0.0401    0.0274    0.0180    0.0105    0.0045   -0.0006   -0.0051   -0.0093
%
%  Columns 13 through 24
%
%   -0.0133   -0.0172   -0.0211   -0.0252   -0.0295   -0.0340   -0.0387   -0.0438   -0.0492   -0.0549   -0.0609   -0.0673
%
%  Columns 25 through 36
%
%   -0.0741   -0.0811   -0.0884   -0.0958   -0.1035   -0.1112   -0.1188   -0.1264   -0.1337   -0.1407   -0.1472   -0.1531
%
%  Columns 37 through 48
%
%   -0.1584   -0.1629   -0.1664   -0.1691   -0.1708   -0.1715   -0.1712   -0.1700   -0.1680   -0.1653   -0.1619   -0.1581
%
%  Columns 49 through 50
%
%   -0.1541   -0.1499
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.1499    0.1541    0.1581    0.1619    0.1653    0.1680    0.1700    0.1712    0.1715    0.1708    0.1691    0.1664
%
%  Columns 13 through 24
%
%    0.1629    0.1584    0.1531    0.1472    0.1407    0.1337    0.1264    0.1188    0.1112    0.1035    0.0958    0.0884
%
%  Columns 25 through 36
%
%    0.0811    0.0741    0.0673    0.0609    0.0549    0.0492    0.0438    0.0387    0.0340    0.0295    0.0252    0.0211
%
%  Columns 37 through 48
%
%    0.0172    0.0133    0.0093    0.0051    0.0006   -0.0045   -0.0105   -0.0180   -0.0274   -0.0401   -0.0575   -0.0827
%
%  Columns 49 through 50
%
%   -0.1217   -0.6393
%
\end{exm}


\begin{exm} Minimize the following tridiagonal sexitic form over $\mathbb{S}^{n-1}$
\[
\sum_{1\leq i \leq n} x_i^6 + \sum_{1\leq i \leq n-1}2x_i^3x_{i+1}^3 =
\bbm x_1^3 \\ x_2^3 \\ x_3^3 \\ \vdots \\ x_{n-1}^3 \\ x_n^3 \ebm^T
\bbm
1 & 1 & 0 & \cdots & 0 \\
1 & 1 & 1 & \ddots & 0 \\
0 & 1 & 1 & \ddots & \vdots \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0 & 0 & \cdots & 1 & 1
\ebm
\bbm x_1^3\\ x_2^3 \\ x_3^3 \\ \vdots \\ x_{n-1}^3 \\ x_n^3 \ebm.
\]
For $n=20$, the size of the resulting SDP  \reff{hmgsdp:prime}-\reff{hmgsdp:dual} is $(N,m)=(1540,177099)$.
Solving \reff{hmgsdp:prime}-\reff{hmgsdp:dual} takes about $10$ hours,
and the computed lower bound $f_{sos}^{hmg}=3.446 \cdot 10^{-5}$ (only four leading digits are shown).
So this form is positive definite, which is contrast to the fact that
the associated tridiagonal matrix is not positive semidefinite.
The optimal $Z^*$ has rank one, and the extracted solution (only four decimal digits are shown) is
\[
\baray{l}
(0.1474,   -0.1811,    0.2024,   -0.2176,    0.2290,   -0.2375,    0.2439,   -0.2484,    0.2514,   -0.2528, \\
\,\, 0.2528,   -0.2514,    0.2484,   -0.2439,    0.2375,   -0.2290,    0.2176,   -0.2024,    0.1811,   -0.1474).
\earay
\]
Its relative error is around $8 \cdot 10^{-11}$.
So it is a global minimizer up to a tiny numerical error.
% Minimize the form:
% sum_{1<=i<=n} x_i^6 +  sum_{2<=i<=n-1}2x_i^3x_{i+1}^3
% over the unit sphere
% n= 20  deg= 6
%
% it    time         gap        lg(rrd)  lg(rrp)   sigma
%  1  05:11:02   -7.49318e-01   -0.364   -7.570        1
%  2  07:11:27   -1.15799e-05   -3.694   -7.895        5
%  3  08:55:03   -1.51730e-07   -5.424   -7.796       25
%  4  09:47:26   -5.36498e-10   -7.702   -7.985      125
%total time:  35246.232
%
%ans =
%
%   1.0e-04 *
%
%    0.3446    0.3446
%
%
%error =
%
%  -8.7597e-11
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.1474   -0.1811    0.2024   -0.2176    0.2290   -0.2375    0.2439   -0.2484    0.2514   -0.2528    0.2528   -0.2514
%
%  Columns 13 through 20
%
%    0.2484   -0.2439    0.2375   -0.2290    0.2176   -0.2024    0.1811   -0.1474
%
%%% an ealier computation using flat extension conditions but not correct due to numerical troubles
%
% Minimize the form:
% sum_{1<=i<=n} x_i^6 +  sum_{2<=i<=n-1}2x_i^3x_{i+1}^3
% over the unit sphere
% n= 20  deg= 6
%
% it    time         gap        lg(rrd)  lg(rrp)   sigma
%  1  05:10:50   -7.49318e-01   -0.364   -7.570        1
%  2  07:11:08   -1.15799e-05   -3.694   -7.895        5
%  3  08:54:24   -1.51730e-07   -5.424   -7.796       25
%  4  09:46:27   -5.36498e-10   -7.702   -7.985      125
%total time:  35187.816
%
%ans =
%
%    0.0000    0.0267
%
%
%error =
%
%   -0.0267
%
%%%  the following extracted solutions are not correct due to numerical errors
%%%  the optimal Z^* has rank one ( numerically )
%%%
%
%there are totally  25 extracted solutions
%ans =
%
%  Columns 1 through 12
%
%    0.0056   -0.0540    0.1110    0.0007    0.0232   -0.0030   -0.0107    0.0198   -0.1972    0.4419    0.2386   -0.2678
%
%  Columns 13 through 20
%
%    0.0543   -0.0935    0.1617   -0.2368    0.2697   -0.3772    0.2205   -0.5132
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0151    0.2213   -0.1418    0.2685   -0.3358    0.4074   -0.3541    0.2554    0.0252   -0.1038   -0.4446    0.2963
%
%  Columns 13 through 20
%
%   -0.0139    0.0784   -0.0944    0.0855   -0.0657   -0.1916   -0.0762   -0.1513
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0057    0.1739   -0.1826    0.1246   -0.1823    0.1631   -0.1865    0.2216   -0.0140   -0.2747   -0.5245    0.1662
%
%  Columns 13 through 20
%
%   -0.0143   -0.0656    0.1770   -0.2467    0.2968   -0.3094    0.3128   -0.1426
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0326   -0.1717    0.1014   -0.2582    0.2560   -0.2752    0.2798   -0.3454    0.2720    0.2553    0.4212   -0.4033
%
%  Columns 13 through 20
%
%    0.1045   -0.1799    0.0509    0.0057   -0.0425    0.0787   -0.0321    0.1395
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0354    0.0061   -0.0373   -0.1178    0.0800   -0.1016    0.1570   -0.2866    0.2232    0.3693    0.2550   -0.4296
%
%  Columns 13 through 20
%
%    0.4002   -0.3382    0.2511   -0.2096    0.1383   -0.0445    0.0974    0.1100
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0357    0.1977   -0.1887    0.1249   -0.1902    0.1709   -0.0993   -0.0218   -0.0994   -0.5100   -0.6173    0.1396
%
%  Columns 13 through 20
%
%    0.2354   -0.1837    0.0927   -0.0643    0.0900   -0.1378    0.1068    0.1471
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0166    0.0945   -0.1282    0.0052   -0.0495    0.0311    0.0626   -0.2917    0.1562    0.3974    0.1489   -0.3437
%
%  Columns 13 through 20
%
%    0.3598   -0.3927    0.3317   -0.2594    0.1812    0.0070    0.1291    0.2209
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0566    0.1320   -0.1386    0.1567   -0.1785    0.2287   -0.2465    0.2705   -0.1326   -0.3525   -0.5247    0.3401
%
%  Columns 13 through 20
%
%   -0.2455    0.1834   -0.1383    0.0996   -0.0761    0.0021    0.0571   -0.2295
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0454    0.1365   -0.1364    0.1459   -0.1918    0.2232   -0.2587    0.3488   -0.2321   -0.1625   -0.3625    0.4345
%
%  Columns 13 through 20
%
%   -0.2742    0.2903   -0.2098    0.1311   -0.0786   -0.0282   -0.0371   -0.1864
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.1175    0.1840   -0.1523    0.2253   -0.2168    0.2878   -0.3188    0.2143   -0.1256    0.0723   -0.5467    0.1753
%
%  Columns 13 through 20
%
%   -0.3080    0.1825   -0.1531    0.0794   -0.1058    0.2040   -0.1669    0.1087
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0112    0.0672   -0.1542    0.0613   -0.0747    0.0715   -0.1083    0.1037   -0.1011   -0.3454   -0.1646    0.2397
%
%  Columns 13 through 20
%
%   -0.2928    0.3267   -0.3601    0.3559   -0.3600    0.2484   -0.2773    0.0595
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0061    0.1434   -0.1246    0.1268   -0.1842    0.1995   -0.2065    0.2274   -0.2701    0.2005   -0.0142    0.2397
%
%  Columns 13 through 20
%
%   -0.3161    0.3260   -0.2950    0.2855   -0.3019    0.3033   -0.2256    0.0114
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0636    0.0377    0.0273    0.0864   -0.1328    0.1511   -0.2669    0.2573   -0.2515   -0.0802    0.1268    0.5437
%
%  Columns 13 through 20
%
%    0.0009   -0.0836    0.2006   -0.2391    0.2905   -0.3641    0.2421   -0.2170
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0100   -0.2323    0.2019   -0.1644    0.2038   -0.1610    0.0947   -0.0228   -0.0036   -0.2174   -0.0611   -0.3335
%
%  Columns 13 through 20
%
%   -0.2176    0.2587   -0.3093    0.2947   -0.3213    0.3249   -0.3049    0.2009
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0446    0.0037    0.0239    0.0192    0.0807   -0.0315    0.1323   -0.2088    0.0807    0.0030    0.4081    0.1061
%
%  Columns 13 through 20
%
%    0.2997   -0.3309    0.3762   -0.3574    0.3734   -0.2037    0.2453   -0.1857
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0409    0.2007   -0.0300   -0.0557   -0.1048    0.2058   -0.0791    0.2149    0.0229   -0.1564   -0.8260   -0.2053
%
%  Columns 13 through 20
%
%   -0.0828    0.0001   -0.0852   -0.0135   -0.0702    0.0565   -0.1009    0.2575
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.1903   -0.1704    0.1934   -0.2476    0.2065   -0.2274    0.1908   -0.2091    0.2418   -0.2076    0.2601   -0.1229
%
%  Columns 13 through 20
%
%    0.2736   -0.2888    0.2615   -0.2608    0.2453   -0.2585    0.1774    0.1555
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0156    0.3149   -0.2653    0.2116   -0.3009    0.3452   -0.2459    0.2737   -0.1399    0.1987   -0.2795    0.2901
%
%  Columns 13 through 20
%
%    0.0262   -0.0740    0.1156   -0.1562    0.1929   -0.2955    0.2140   -0.0689
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.1024    0.0023   -0.0016    0.3071    0.1516   -0.2166    0.1950    0.0189    0.2503    0.1826   -0.5752   -0.4721
%
%  Columns 13 through 20
%
%    0.0081   -0.0572   -0.0132    0.0652    0.0424   -0.2578    0.2420   -0.0515
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0104   -0.0096    0.0273   -0.0888    0.0777   -0.0853    0.0785   -0.0562    0.1549   -0.2766   -0.0804    0.0776
%
%  Columns 13 through 20
%
%    0.2508   -0.2623    0.3639   -0.4287    0.3910   -0.2126    0.3345   -0.3126
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0031   -0.2445    0.2263   -0.3335    0.2722   -0.3110    0.3082   -0.2226    0.2149   -0.1229    0.4460   -0.1755
%
%  Columns 13 through 20
%
%    0.1681   -0.1777    0.1693   -0.1793    0.1651   -0.1047    0.1007   -0.0496
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0137   -0.3157    0.3107   -0.2954    0.1378   -0.0580   -0.1092   -0.0193   -0.0984    0.2693    0.3479   -0.0790
%
%  Columns 13 through 20
%
%   -0.2777    0.2340   -0.2396    0.4392   -0.2708    0.0954    0.0405   -0.0767
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0418   -0.0419    0.0483   -0.1017    0.0524   -0.2561    0.2687   -0.2845    0.1910    0.1870    0.6666    0.1294
%
%  Columns 13 through 20
%
%    0.2015   -0.2262    0.1873   -0.0368    0.1543   -0.2620    0.0917    0.0279
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.1460   -0.0092    0.0528   -0.3873    0.1277   -0.2223    0.2457   -0.1639    0.1954   -0.1588    0.2366   -0.2841
%
%  Columns 13 through 20
%
%    0.2216   -0.2473    0.1701    0.3650    0.0825   -0.3752    0.2213   -0.0616
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0102    0.1766   -0.4420    0.3443   -0.0504    0.0719   -0.1395   -0.2009    0.2553   -0.5064   -0.1473   -0.0137
%
%  Columns 13 through 20
%
%    0.2668   -0.2431    0.1587   -0.1700    0.0805   -0.1358   -0.0195    0.1846
%
\end{exm}




\begin{exm} Minimize the following sexitic form over $\mathbb{S}^{n-1}$
\[
\sum_{ 1 \leq i < j < k \leq n} x_i^2x_j^2x_k^2 + x_i^3x_j^2x_k +
x_i^2x_j^3x_k + x_ix_j^3x_k^2.
\]
For $n=20$, the resulting SDP \reff{hmgsdp:prime}-\reff{hmgsdp:dual} has size $(N,m)=(1540,177099)$.
Solving \reff{hmgsdp:prime}-\reff{hmgsdp:dual} takes about $11$ hours.
The computed lower bound $f_{sos}^{hmg}=-0.3827$ (only four decimal digits are shown).
The optimal $Z^*$ has rank one, and the extracted solution (only four decimal digits are shown) is
\[
\baray{l}
-(-0.6551,   -0.6286,   0.1109,   0.1089,   0.1071,   0.1054,  0.1038,  0.1023,  0.1009,  0.0996, \\
\qquad 0.0983, 0.0971,  0.0960,  0.0949,  0.0939,  0.0929,  0.0920,  0.0911,  0.0902,  0.0894).
\earay
\]
Its relative error is around $6\cdot 10^{-11}$.
So it is a global minimizer up to a tiny numerical error.
%
%
% Minimize the form:
% sum_{i<j<k}  x_i^2*x_j^2*x_k^2 + x_i^3*x_j^2*x_k + x_i^2*x_j^3*x_k^1 + x_i^1*x_j^3*x_k^2
% over the unit sphere
% n= 20  deg= 6
%
% it    time         gap        lg(rrd)  lg(rrp)   sigma
%  1  04:44:30   -3.31322e+03    1.459   -8.518        1
%  2  06:42:39   -2.00629e-03   -3.063   -8.787        5
%  3  07:10:44   -2.56144e-03   -3.246   -8.386       25
%  4  07:33:26   -1.32558e-03   -3.786   -8.246      125
%  5  08:30:51   -1.39768e-04   -4.894   -8.719      625
%  6  10:57:30   -2.56005e-06   -6.662   -7.650     3125
%total time:  39450.680
%
%ans =
%
%   -0.3827   -0.3827
%
%
%error =
%
%   6.1316e-11
%
%there are totally   1 extracted solutions
%ans =
%
%  Columns 1 through 12
%
%    0.6551    0.6286   -0.1109   -0.1089   -0.1071   -0.1054   -0.1038   -0.1023   -0.1009   -0.0996   -0.0983   -0.0971
%
%  Columns 13 through 20
%
%   -0.0960   -0.0949   -0.0939   -0.0929   -0.0920   -0.0911   -0.0902   -0.0894
%
%
\end{exm}


An important application of homogenous polynomial optimization \reff{hmgpop:sphere}
is estimating stability numbers of graphs, as shown below.


\begin{exm}[Stability number of graphs]
Let $G=(V,E)$ be a graph with $|V|=n$. The stability number $\af(G)$ is the cardinality of
the biggest stable subset(s) (their vertices are not connected by any edges) of $V$.
It was shown in Motzkin and Straus \cite{MS65}
(also see De Klerk and Pasechnik \cite{dKP02}) that
\[
\af(G)^{-1} = \min_{x\in \Dt_n} \quad x^T(A+I_n)x,
\]
where $\Dt_n$ is the standard simplex in $\re^n$
and $A$ is the adjacency matrix associated with $G$.
If replacing every $x_i\geq 0$ by $x_i^2$, we get
\be \label{af(G):M-S}
\af(G)^{-1} = \min_{\|x\|_2 = 1} \quad \sum_{i=1}^n x_i^4
+ 2 \sum_{(i,j)\in E} x_i^2x_j^2 .
\ee
This is a quartic homogeneous polynomial optimization.
When a lower bound $f_{sos}^{hmg}$ for \reff{af(G):M-S} via its SOS relaxation \reff{hmgpop:sos}
is computed, we round $\left(f_{sos}^{hmg}\right)^{-1}$ to the nearest integer
which will be used as an estimate for $\af(G)$.



We generate random graphs $G$, and solve the SDP  \reff{hmgsdp:prime}-\reff{hmgsdp:dual}
which is equivalent to its SOS relaxation by applying Algorithm~\ref{alg:Newton-CG}.
The generation of random graphs is in the similar way
as in Bomze and De Klerk \cite[Section~6]{BdK02}.
For $n=20, 30, 40, 50, 60$, we generate random graphs $G=(V,E)$ with $|V|=n$.
Select a random subset $M\subset V$ with $|M|=n/2$.
The edges $e_{ij}(\{i,j\} \not\subset M)$ are generated with probability $\half$.
The computational results are in Table \ref{stab:test}.
The first column is the number of vertices,
the second column is the size of the resulting SDP \reff{hmgsdp:prime}-\reff{hmgsdp:dual},
the third column is the number of generated instances,
and the last column is the minimum, median and maximum time consumed by the computer.



As we can see from Table \ref{stab:test}, for $n=20$,
solving \reff{hmgsdp:prime}-\reff{hmgsdp:dual} takes a couple of minutes;
for $n=30$, it takes about four to ninety minutes;
for $n=40$, it takes about twenty minutes to five hours;
for $n=50$, it takes about twenty five hours;
for $n=60$, it takes about seventy five hours.
In all the instances, we get correct stability numbers
via solving the SOS relaxation of \reff{af(G):M-S}.



\begin{table}
\centering
\begin{tabular}{|r||c||c||r|r|r|} \hline
$n$ & (N,m) &\# inst   & \multicolumn{3}{c|}{time (min, med, max)}\\ \hline
20 &  (210, 8854)  & 20  & 0:00:22 & 0:03:04  & 0:06:25\\ \hline
30 &  (465, 40919) & 10  & 0:04:12 & 0:39:58  & 1:33:47 \\ \hline
40 &  (820, 123409)& 10  & 0:20:52 & 1:21:01  & 5:24:12 \\ \hline
50 & (1275, 292824)& 3  & 23:50:02 & 26:27:12 & 27:15:46  \\ \hline
60 & (1830, 595664)& 1  & 75:46:58 & 75:46:58 & 75:46:58 \\ \hline
%70 & (2485, 1088429)&  &  &  &  \\ \hline
%80 & (3240, 1837619)&  &  &  &  \\ \hline
\end{tabular}
\caption{Computational results for stability number of random graphs }
\label{stab:test}

\end{table}


\end{exm}




\begin{exm}[Random forms]
We test the performance of Algorithm~\ref{alg:Newton-CG}
in solving SOS relaxation \reff{hmgpop:sos} for random forms.
Generate a form $f(x)$ randomly as
\[
f(x) = \sum_{\af\in \N^n: |\af|=2d }  f_\af x^\af,
\]
where the coefficients $f_\af$ are Gaussian random variables.
We apply Algorithm~\ref{alg:Newton-CG} to solve the SDP \reff{hmgsdp:prime}-\reff{hmgsdp:dual}
which is equivalent to the corresponding SOS relaxation \reff{hmgpop:sos}.
The random forms of degrees $4,6,8,10$ are tested.
The computational results are in Tables~\ref{hmgtab:rd:dg=4}-\ref{hmgtab:rd:dg=10}.
The meaning of their columns are the same as in
Tables~\ref{UCtab:rand:deg=4}-\ref{UCtab:rand:deg=10}.
%The first column is the number of variables,
%the second column is the size of the resulting SDP \reff{hmgsdp:prime}-\reff{hmgsdp:dual},
%the third column is the number of generated instances,
%and the fourth column is the minimum, median and maximum time consumed by the computer,
%and the last column is the minimum, median and maximum errors of obtained solutions.



When $f(x)$ has degree $4$ ($d=2$),
the computational results are in Table~\ref{hmgtab:rd:dg=4}.
For $n=30$, solving \reff{hmgsdp:prime}-\reff{hmgsdp:dual} takes a couple of minutes;
for $n=40$, it takes about thirty to forty minutes;
for $n=50$, it takes a few hours;
for $n=60$, it takes about six to nineteen hours;
for $n=70$, it takes about fifteen to thirty four hours;
for $n=80$, it takes about one hundred hours.
In all the instances, the optimal $Z^*$ has rank one and FEC holds.
The obtained solutions have tiny relative errors.



When $f(x)$ has degree $6$ ($d=3$),
the computational results are in Table~\ref{hmgtab:rd:dg=6}.
For $n=15$, solving \reff{hmgsdp:prime}-\reff{hmgsdp:dual} takes around
seventeen to twenty six minutes;
for $n=20$, it takes a couple of hours;
for $n=25$, it takes about fifty hours.
The optimal $Z^*$ now does not always satisfy FEC.
When $n=10$ or $15$, there are some instances that
their errors are relatively bigger than others.
%We are not sure whether this is because of the numerical errors occurred during the computation,
%or the inexactness of SOS relaxation \reff{hmgpop:sos}.
When $n=20$, there are a few instances that their relative errors are significantly much bigger.
This is probably due to the inexactness of their SOS relaxations.
For $n=25$, we have generated only one instance since its computation takes more than two days,
and luckily a solution of high accuracy was found.
But we should not expect this generally.


\begin{table}[htb]
\centering
\begin{tabular}{|r||c||c||r|r|r||c|} \hline
$n$ & (N,m) & \# inst & \multicolumn{3}{c||}{time (min, med, max)} & err (min, med, max) \\ \hline
20 &  (210, 8854)   & 20 & 0:00:26 & 0:00:31 & 0:00:39 & (2e-9, 9e-9, 3e-8) \\ \hline
30 &  (465, 40919)  & 20 & 0:06:04 & 0:06:38 & 0:08:41 & (4e-9, 8e-9, 2e-8) \\ \hline
40 &  (820, 123409) & 10 & 0:35:54 & 0:36:03 & 0:41:01 & (4e-10, 7e-9, 3e-8) \\ \hline
50 & (1275, 292824) & 10 & 2:29:23 & 2:35:09 & 3:58:45 & (5e-8, 7e-7, 3e-7) \\ \hline
60 & (1830, 595664) & 3 & 6:31:22 &  9:02:36 & 19:07:32 & (9e-9, 5e-8, 1e-8) \\ \hline
70 & (2485, 1088429)& 3 & 15:48:39 & 24:57:13  & 34:03:29 & (1e-10, 6e-9, 1e-9) \\ \hline
80 & (3240, 1837619)& 1 & 99:26:20 & 99:26:20 & 99:26:20 & (7e-8, 2e-8, 9e-7) \\ \hline
%90 & (4095, 2919734)& &  &  &  & \\ \hline
%100& (5050, 4421274)& &  &  &  & \\ \hline
\end{tabular}
\caption{Computational results for random \reff{hmgpop:sphere} of degree $4$}
\label{hmgtab:rd:dg=4}

\end{table}



\begin{table}[htb]
\centering
\begin{tabular}{|r||c||c||r|r|r||c|} \hline
$n$ & (N,m) & \# inst & \multicolumn{3}{c||}{time (min, med, max)} & err (min, med, max) \\ \hline
10 & (220, 5004)    & 20 & 0:00:21 & 0:00:25 & 0:00:30 & (1e-9,4e-7,5e-5) \\ \hline
15 & (680, 38759)   & 20 &  0:17:02 & 0:18:24 & 0:26:27 &(6e-9,9e-7,3e-4)  \\ \hline
20 & (1540, 177099) & 5 & 5:49:08 & 5:51:44  & 6:12:59 &(5e-11,5e-9,3.63) \\ \hline
25 & (2925, 593774) & 1 & 51:05:51  & 51:05:51 &  51:05:51  &  (4e-9, 4e-9, 4e-9)  \\ \hline
\end{tabular}
\caption{Computational results for random \reff{hmgpop:sphere} of degree $6$}
\label{hmgtab:rd:dg=6}

\end{table}


\begin{table}[htb]
\centering
\begin{tabular}{|r||c||c||r|r|r||c|} \hline
$n$ & (N,m) & \# inst & \multicolumn{3}{c||}{time (min, med, max)} & err (min, med, max) \\ \hline
8 &  (330, 6434)    & 20 & 0:01:56  & 0:02:28  & 0:03:11
& (4e-9,0.85, 2.08)\\ \hline
10 & (715, 24309)   & 20 & 0:29:18  & 0:33:55  & 0:40:24  &
(5e-11,0.02,2.84)\\ \hline
13 & (1820, 125969) & 3 & 10:36:40  & 16:35:53  & 18:14:41  &
(1.64, 1.91, 3.26)\\ \hline
15 & (3060, 319769) & 1 & 104:41:59   & 104:41:59  & 104:41:59   &  (1.24, 1.24, 1.24) \\ \hline
\end{tabular}
\caption{Computational results for random \reff{hmgpop:sphere} of degree $8$}
\label{hmgtab:rd:dg=8}
\end{table}



\begin{table}[htb]
\centering
\begin{tabular}{|r||c||c||r|r|r||c|} \hline
$n$ & (N,m) & \# inst & \multicolumn{3}{c||}{time (min, med, max)} & err (min, med, max) \\ \hline
6 &  (252, 3002)     & 20 & 0:00:43 & 0:00:57 & 0:01:16  &
(3e-9,0.98, 2.24) \\ \hline
8 &  (792, 19447)    & 10 &  0:36:26 &  0:59:51& 1:28:46  &
   (0.14, 1.67, 2.15) \\ \hline
10 & (2002, 92377)   & 3 &  25:58:15 & 29:30:35 & 30:37:08  &
 (3e-3, 0.87,  2.51) \\ \hline
11 & (3003, 184755)  & 1 & 110:48:52  & 110:48:52 & 110:48:52  &  (1.27, 1.27, 1.27) \\ \hline
\end{tabular}
\caption{Computational results for random \reff{hmgpop:sphere} of degree $10$}
\label{hmgtab:rd:dg=10}

\end{table}



When $f(x)$ has degree $8$ ($d=4$),
the computational results are in Table~\ref{hmgtab:rd:dg=8}.
For $n=8$, solving \reff{hmgsdp:prime}-\reff{hmgsdp:dual} takes a few minutes;
for $n=10$, it takes about thirty to forty minutes;
for $n=13$, it takes about ten to eighteen hours;
for $n=15$, it takes about one hundred hours.
The optimal $Z^*$ usually does not satisfy FEC,
and the solutions obtained have very low accuracy.
For a few instances, very accurate solutions are found.
However, for most instances, the obtained solutions have big errors.
This phenomena is probably due to the inexactness of their SOS relaxations.







When $f(x)$ has degree $10$ ($d=5$),
the computational results are in Table~\ref{hmgtab:rd:dg=10}.
For $n=8$, solving \reff{hmgsdp:prime}-\reff{hmgsdp:dual} takes about thirty to ninety minutes;
for $n=10$, it takes about twenty five to thirty hours;
for $n=11$, it takes about one hundred hours. Similar to the case of degree $8$,
the optimal $Z^*$ usually does not satisfy FEC.
For most instances, the solutions obtained have big errors.
This is probably because of the inexactness of their SOS relaxations.




From the above, for randomly generated forms,
we know SOS relaxation \reff{hmgpop:sos} is very successful when $\deg(f)=4$
and also works well when $\deg(f)=6$.
However, when $\deg(f)=8$ or $10$,
SOS relaxation \reff{hmgpop:sos} usually fails to find the true global minimum.
This observation is consistent with the discovery of Blerkherman \cite{Blek}.
On a regular computer, Algorithm~\ref{alg:Newton-CG} would solve
big homogeneous polynomial optimization problems
which are not solvable by prior existing SOS solvers.
This is a big advantage of regularization type methods.


\end{exm}




\subsection{Minimizing odd forms over unit spheres}


In some applications, we might need to optimize an odd form $f(x)$ ($\deg(f)$ is odd)
over $\mathbb{S}^{n-1}$. This problem is usually also very difficult.
Nesterov \cite{Nes03} proved that it is already NP-hard when $\deg(f)=3$.
For odd forms, SOS relaxation \reff{hmgpop:sos} can not be applied directly.
However, after a reformulation,
this problem is equivalent to minimizing certain even forms.

Let $f(x)$ be an odd form of degree $2d-1$,
and $f_{min}^{hmg}$ be its minimum on $\mathbb{S}^{n-1}$.
Let $\hat{f}(x,t) = f(x) t $ be a new even form in $(x,t)$.
Consider homogeneous optimization
\be \label{opt:odd->even}
\widehat{f_{min}^{hmg}} := \min_{\|x\|_2^2 + t^2 = 1}  f(x)t.
\ee
%\[
%\hat{f}_{min} = \min_{\|x\|_2^2 + t^2 = 1}  f(x)t,  \qquad
%\hat{f}_{max} = \max_{\|x\|_2^2 + t^2 = 1}  f(x)t.
%\]
%Note the following relations
%\begin{align*}
%&\hspace{2cm} \min_{ 0\leq t \leq 1}  \min_{ \|x\|_2 = \sqrt{1-t^2}} f(x) t =
%\min_{ 0\leq t \leq 1} t  \min_{ \|x\|_2 = \sqrt{1-t^2}} f(x) = \\
%& \min_{ 0\leq t \leq 1} (t (1-t^2)^{d-1/2})  f_{min} =
%f_{min} \max_{ 0\leq t \leq 1} (t (1-t^2)^{d-1/2}) =
%\frac{1}{\sqrt{2d-1}}\left(1-\frac{1}{2d}\right)^d f_{min},
%\end{align*}
%\begin{align*}
%&\hspace{2cm} \min_{ -1\leq t \leq 0}  \min_{ \|x\|_2 = \sqrt{1-t^2}} f(x) t =
%\min_{ 0\leq t \leq 1} t  \max_{ \|x\|_2 = \sqrt{1-t^2}} f(x) = \\
%& \min_{ -1\leq t \leq 0} (t (1-t^2)^{d-1/2})  f_{max} =
%f_{max} \min_{ -1\leq t \leq 0} (t (1-t^2)^{d-1/2}) =
%\frac{1}{\sqrt{2d-1}}\left(1-\frac{1}{2d}\right)^d f_{min}.
%\end{align*}
%Thus we have
It was shown in \cite{Nie09b} that
\[
f_{min}^{hmg}  = \sqrt{2d-1}\left(1-\frac{1}{2d}\right)^{-d} \widehat{f_{min}^{hmg}}.
\]
Since the form $\hat{f}(x,t) = f(x) t $ is even, SOS relaxation \reff{hmgpop:sos}
can be applied to get a lower bound $\widehat{f_{sos}^{hmg}}$ for $\widehat{f_{min}^{hmg}} $.
Then
\[
f_{sos}^{hmg}  := \sqrt{2d-1}\left(1-\frac{1}{2d}\right)^{-d} \widehat{f_{sos}^{hmg}}
\]
is a lower bound for $f_{min}^{hmg}$.
When a minimizer $(\hat{x}, \hat{t})$ is obtained for \reff{opt:odd->even},
the normalization of $\hat{x}/\hat{t}$
is then a minimizer of $f(x)$ over $\mathbb{S}^{n-1}$.


\begin{exm} Minimize the following cubic form over $\mathbb{S}^{n-1}$
\[
f(x) = \sum_{1\leq i < j < k \leq n} (x_ix_jx_k + x_i^2x_j - x_i^2x_k + x_jx_k^2).
\]
We first transform it to \reff{opt:odd->even}, and then solve
its SOS relaxation by Algorithm~\ref{alg:Newton-CG}.
For $n=49$, the computation takes about four hours,
and the lower bound found is $-124.9645$ (only four decimal digits are shown).
The optimal $Z^*$ has rank one, and the extracted solution (only four decimal digits are shown) is
\[
-(0.0320,   0.0344,   0.0368,   0.0394,   0.0420,
\ldots, 0.2304,   0.2434,   0.2600,   0.2836,   0.3238).
\]
\iffalse
%% the following is the full solution
\[
\baray{l}
-(0.0320,   0.0344,   0.0368,   0.0394,   0.0420,   0.0447,   0.0474,   0.0502,   0.0531,   0.0561, \\
\quad\, 0.0592,   0.0623,   0.0655,   0.0688,   0.0721,   0.0755,   0.0790,   0.0826,   0.0863,   0.0900, \\
\quad\, 0.0938,   0.0977,   0.1017,   0.1057,   0.1099,   0.1141,   0.1184,   0.1228,   0.1274,   0.1320, \\
\quad\, 0.1367,   0.1416,   0.1466,   0.1518,   0.1571,   0.1626,   0.1683,   0.1743,   0.1805,   0.1872, \\
\quad\, 0.1942,   0.2019,   0.2102,   0.2196,   0.2304,   0.2434,   0.2600,   0.2836,   0.3238).
\earay
\]
\fi
Its relative error is about $10^{-9}$.
So it is a global minimizer up to a tiny numerical error.
% Minimize the odd form:
%  sum_{i<j<k}  x_i*x_j*x_k +  x_i^2*x_j - x_i^2*x_k + x_j*x_k^2
%  over the unit sphere
% n= 49  deg= 3
%
% it    time         gap        lg(rrd)  lg(rrp)   sigma
%  1  03:02:23   -1.52086e+03    1.290   -8.685        1
%  2  03:28:52   -7.36084e-03   -2.533   -9.592        5
%  3  03:35:14   -1.22344e-02   -2.666   -8.567       25
%  4  03:47:37   -9.93761e-03   -3.097   -8.831      125
%  5  03:57:44   -1.67619e-03   -4.050   -7.637      625
%  6  04:05:44   -4.78049e-05   -5.646   -8.775     3125
%  7  04:11:11   -2.50723e-07   -7.930   -8.470    15625
%total time:  15071.740
%
%ans =
%
% -124.9645 -124.9645
%
%
%error =
%
%  -1.1277e-09
%
%
%ans =
%
%  Columns 1 through 12
%
%    0.0320    0.0344    0.0368    0.0394    0.0420    0.0447    0.0474    0.0502    0.0531    0.0561    0.0592    0.0623
%
%  Columns 13 through 24
%
%    0.0655    0.0688    0.0721    0.0755    0.0790    0.0826    0.0863    0.0900    0.0938    0.0977    0.1017    0.1057
%
%  Columns 25 through 36
%
%    0.1099    0.1141    0.1184    0.1228    0.1274    0.1320    0.1367    0.1416    0.1466    0.1518    0.1571    0.1626
%
%  Columns 37 through 48
%
%    0.1683    0.1743    0.1805    0.1872    0.1942    0.2019    0.2102    0.2196    0.2304    0.2434    0.2600    0.2836
%
%  Column 49
%
%    0.3238
%
\end{exm}






%
% the following section was skipped
%

\iffalse

\section{Bi-quadratic optimization over unit spheres}
\setcounter{equation}{0}

In this section,
we consider the bi-quadratic optimization
\be \label{bq:sphere}
\baray{rl}
f_{min}^{bqp} :=
\min & f(x) \\
s.t. & \|x^{I}\|_2 = \|x^{II}\|_2 = 1.
\earay
\ee
Here $f(x)$ is a homogenous bi-quadratic form given as
\be \label{bqform}
f(x) = \sum_{(\af,\bt) \in \N^{n_1}\times \N^{n_2}: |\af| = |\bt| =2}
f_{\af,\bt} \cdot (x^{I})^\af \cdot (x^{II})^\bt
\ee
where $x=(x^{I}, x^{II})$, $x^{I} \in \re^{n_1}$, $x^{II} \in \re^{n_2}$
and $n=n_1+n_2$, that is,
\[
x=(x^{I}, x^{II}), \quad
x^{I}=(x_1,\ldots,x_{n_1}), \quad
x^{II}=(x_{n_1+1}, \ldots, x_n).
\]
Note that $f(x)$ is quadratic in $x^{I}$ and $x^{II}$ respectively,
but the total degree is four.
Every bi-quadratic form given by \reff{bqform} defines
a linear map $\mc{F}$ from the symmetric matrix space $\mc{S}_{n_1}$
to the symmetric matrix space $\mc{S}_{n_2}$ as
\be
\mc{F}: \quad
X=(X_{ij})_{1\leq i,j\leq n_1} \mapsto
Y=\left(\sum_{1\leq i,j\leq n_1} f_{e_i+e_j, e_k+e_l} X_{ij}\right)_{1\leq k,l\leq n_2}.
\ee
The linear map $\mc{F}$ is called positive if
$\mc{F}(X)$ is positive semidefinite for any positive semidefinite $X$.
Thus the linear map $\mc{F}$ is positive if and only if
the associated bi-quadratic form $f(x)$ is nonnegative for all $x$.


Bi-quadratic optimization \reff{bq:sphere} has important applications
in quantum science, entanglement problem,
higher-order statistics, and independent component analysis.
Ling et al. \cite{LNQY} proved that problem \reff{bq:sphere} is NP-hard,
and proposed various approximation algorithms for solving \reff{bq:sphere},
as well as analyzing their approximation performance.
We refer to \cite{LNQY} for an introduction in this area.
One efficient method in these approximation approaches is
the standard SOS relaxation
\be  \label{bq:sos}
\baray{rl}
f_{sos}^{bqp} :=
\max & \gamma \\
s.t. & f(x) - \gamma \cdot \|x^{I}\|_2^2\|x^{II}\|_2^2  \quad \mbox{ is SOS}.
\earay
\ee
The optimal value $f_{sos}^{bqp}$ above is a lower bound
for the global minimum $f_{min}^{bqp}$ of \reff{bq:sphere}.




As indicated in \cite{LNQY},
the SOS relaxation usually provides very good approximations to \reff{bq:sphere}.
Among the currently existing approximation methods,
the SOS relaxation is probably the best in the view of approximation quality.
But it is limited to solving small scale problems due to the expensive cost
of interior-point algorithms.
By SDP software like SeDuMi, it is very difficult to solve \reff{bq:sphere}
when $n_1,n_2 > 15$ on a regular computer,
due to the hight cost of memory and computations.
Here we propose Algorithm \ref{alg:Newton-CG} to solve the SOS relaxation.
As illustrated later in this section,
we are able to solve the SOS relaxation quite efficiently
for the bi-quadratic optimization problem \reff{bq:sphere} having much bigger sizes.
From our numerical experience,
we can solve \reff{bq:sos} for $n_1=n_2 = 50$
within  hours.
This significantly broaden the practical applications
of SOS relaxation \reff{bq:sos}.









Let $[x]_{bl}$ denote the column vector of bilinear monomials
$x_ix_j\,(1\leq i\leq n_1, n_1+1\leq j\leq n$)
ordered lexicographically as
\[
[x]_{bl}^T = \bbm x_1 x^{II} & x_2 x^{II} & \cdots & x_{n_1}x^{II} \ebm.
\]
Note that
\[
\|x^{I}\|_2^2\|x^{II}\|_2^2  =
\sum_{1\leq i\leq n_1, n_1+1\leq j \leq n } (x_ix_j)^2 = [x]_{bl}^T[x]_{bl}.
\]
Define the index set
\[
\Omega_d^n = \{(\af,\bt) \in \N^{n_1} \times \N^{n_2}:
|\af| = |\bt| = d \}.
\]
Define $0/1$ matrices $B_\af$ such that
\be  \label{bq:B_afbt}
[x]_{bl}[x]_{bl}^T = \sum_{(\af,\bt) \in \Omega_2^n }
B_{\af,\bt} \cdot (x^{I})^{\af}\cdot (x^{II})^{\bt}.
\ee
Then the condition that $f(x)-\gamma \cdot [x]_{bl}^T[x]_{bl}$ being SOS
is equivalent to the system that
\[
\baray{rcl}
f_{\af,\bt} - \gamma &=& B_{\af,\bt} \bullet X, \quad \forall\,
(\af,\bt) \in 2\Omega_1^n, \\
f_{\af,\bt} & =& B_{\af,\bt} \bullet X, \quad \forall \,
(\af,\bt) \in \Omega_{2}^n\backslash 2\Omega_1^n,\\
X & \succeq & 0.
\earay
\]
Choosing $(\af,\bt) = (2e_1^{I}, 2e_1^{II})$ in the above, we get
\[
\gamma = f_{(2e_1^{I}, 2e_1^{II})} -  B_{(2e_1^{I}, 2e_1^{II})} \bullet X.
\]
Here $e_i^{I}$ (resp. $e_i^{II})$) denotes the $i$-th standard
unit vector in $\re^{n_1}$ (resp. $\re^{n_2}$).
For convenience, denote $\mathbb{B}^n = \Omega_{2}^n\backslash \{(2e_1^{I}, 2e_1^{II})\}$.
Define vector $s=(s_{\af,\bt})_{(\af,\bt)\in \mathbb{B}^n}$ as follows:
\be \label{bq-def:s}
s_{\af,\bt} =
\bca
1  & \text{ if } \af = \bt, \\
0, & \text{ if } \af \ne \bt .
\eca
\ee
Then define matrices $C, A_\af$ and scalars $b_\af$ as follows
\be \label{matdef:bq}
\baray{r}
C= B_{(2e_1^{I}, 2e_1^{II})}, \quad
A_{\af,\bt} = B_{\af,\bt} - s_{\af,\bt} B_{(2e_1^{I}, 2e_1^{II})}, \\
b_{\af,\bt} = f_{\af,\bt} - s_{\af,\bt} f_{(2e_1^{I}, 2e_1^{II})},
\earay \quad \forall \, (\af,\bt)\in \mathbb{B}^n.
\ee
Let $b=(b_\af)_{\af\in \mathbb{B}^n}$.
Define linear operators
$\mc{B}, \mc{A}: \mc{S}^{|\Omega_{d}^n|} \rightarrow \re^{|\mathbb{B}^n|}$
and the cone $\mc{A}$ as
\[
\mc{B}(X) = (B_\af \bullet X)_{\af\in \mathbb{B}^n}, \quad
\mc{A}(X) = (A_\af \bullet X)_{\af\in \mathbb{B}^n}, \quad
\mc{K} = \mc{S}_+^{|\Omega_{d}^n|}.
\]
For $C,\mc{A},b,\mc{K}$ defined above,
we known the SOS relaxation \reff{bq:sos} is equivalent to
the standard primal SDP problem
\be  \label{bqsdp:prime}
\baray{rl}
f_{sdp}^{bqp} :=
\min & C \bullet X \\
s.t. & \mc{A}(X) = b, \, X \succeq 0.
\earay
\ee
The dual optimization problem of \reff{hmgsdp:prime} becomes
\be  \label{bqsdp:dual}
\baray{rl}
\max & b^Ty \\
s.t. &  C-\mc{A}^*(y) =Z,\, Z \succeq 0.
\earay
\ee



\begin{theorem}
For notations defined above, it holds that
$f_{sos}^{bqp} = -f_{sdp}^{bqp}+f_{(2e_1^{I},2e_1^{II})}$.
If the SOS relaxation \reff{bq:sos} is exact,
then $f_{min}^{bqp} = -f_{sdp}^{bqp}+f_{(2e_1^{I},2e_1^{II})}$.
\end{theorem}
\begin{proof}
The proof is the same as for Theorem~\ref{thm-bd:uc}.
%Optimization problems \reff{bq:sos} and \reff{bqsdp:prime}
%are essentially equivalent. Their objectives satisfy the relation
%\[
%C \bullet X = -\gamma + f_{(2e_1^{I},2e_1^{II})}.
%\]
%Maximizing $\gamma$ in \reff{bq:sos} is equivalent to
%minimizing $C\bullet X$ in \reff{bqsdp:prime}.
%The feasible sets of \reff{bq:sos} and \reff{bqsdp:prime}
%are the same.
%So it must hold $f_{sos}^{bqp} = -f_{sdp}^{bqp}+f_{(2e_1^{I},2e_1^{II})}$.
%
%When the SOS relaxation \reff{bq:sos} is exact,
%that is, $f_{sos}^{bqp} = f_{min}^{bqp}$,
%we immediately have
%$f_{min}^{bqp} = -f_{sdp}^{bqp}+f_{(2e_1^{I},2e_1^{II})}$.
\end{proof}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% the following is for boundary point method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\iffalse

Now we turn to find $(\mc{A}\mc{A}^*)^{-1}$ for $\mc{A}$ defined in \reff{matdef:bq}.
Note that the matrices $B_{\af,\bt}$ defined in \reff{bq:B_afbt} are orthogonal to each other.
Thus $\mc{B}\mc{B}^*$ is a diagonal matrix, that is,
\[
\mc{B}\mc{B}^* = \diag( \mathbf{1}^T B_\af \mathbf{1} : \af \in \mathbb{B}^n ),
\]
where $\mathbf{1}$ is the vector of all ones.
The product $\mathbf{1}^T B_\af \mathbf{1}$ is simply
the number of nonzero entries of $B_\af$.
So the inverse $(\mc{B}\mc{B}^*)^{-1}$ is simply the diagonal matrix
\[
(\mc{B}\mc{B}^*)^{-1} = \diag( (\mathbf{1}^T B_\af \mathbf{1})^{-1} :
\af \in \mathbb{B}^n ).
\]
Thus we have for all $\vartheta^{(1)}=(\af^{(1)},\bt^{(1)}) \in  \mathbb{B}^n,
\vartheta^{(2)}=(\af^{(2)},\bt^{(2)})\in  \mathbb{B}^n$,
\begin{align*}
A_{\vartheta^{(1)}} \bullet A_{\vartheta^{(2)}}
& = (B_{\vartheta^{(1)}} - r_{\vartheta^{(1)}} B_{(2e_1^{I}, 2e_1^{II})} ) \bullet
 (B_{\vartheta^{(2)}} - r_{\vartheta^{(2)}} B_{(2e_1^{I}, 2e_1^{II})}) \\
& = B_{\vartheta^{(1)}} \bullet B_{\vartheta^{(2)}}
- B_{(2e_1^{I}, 2e_1^{II})} \bullet
( r_{\vartheta^{(1)}} B_{\vartheta^{(2)}} + r_{\vartheta^{(2)}} B_{\vartheta^{(1)}} ) +
r_{\vartheta^{(1)}} r_{\vartheta^{(2)}} \\
& = B_{\vartheta^{(1)}} \bullet B_{\vartheta^{(2)}} +
r_{\vartheta^{(1)}} r_{\vartheta^{(2)}}.
\end{align*}
In other words, the above is the same as
\[
\mc{A}\mc{A}^* = \mc{B}\mc{B}^* + ss^T.
\]
By Sherman-Morrison formula, we have
\[
(\mc{A}\mc{A}^*)^{-1} = (\mc{B}\mc{B}^*)^{-1} -
\frac{(\mc{B}\mc{B}^*)^{-1}ss^T(\mc{B}\mc{B}^*)^{-1}}{1+s^T(\mc{B}\mc{B}^*)^{-1}s}.
\]
If define vectors $\mathbf{w}$ and $\mathbf{q}$ such that
\be \label{bq-def:h-p}
\mathbf{w} = ( \, (\mathbf{1}^T B_\af \mathbf{1})^{-1} : \af \in \mathbb{B}^n ), \quad
\mathbf{q} = \frac{1}{\left(1+s^T\diag(\mathbf{w})s\right)^{1/2}}\diag(\mathbf{w})s,
\ee
then we get a simple formula
\[
(\mc{A}\mc{A}^*)^{-1} = \diag(\mathbf{w}) - \mathbf{q}\mathbf{q}^T.
\]

\fi
%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%





Algorithm~\ref{alg:Newton-CG} can be applied to solve
\reff{bqsdp:prime}-\reff{bqsdp:dual} in large scale case.
Let $(X^*,y^*,Z^*)$ be an optimal solution pair for \reff{bqsdp:prime}-\reff{bqsdp:dual}.
Then
\[
-C \bullet X^* + f_{(2e_1^{I},2e_1^{II})}= -X^*(1,1)+f_{(2e_1^{I},2e_1^{II})}
\]
is a lower bound for the global minimum $f_{min}^{bqp}$.
By \reff{bqsdp:dual} and the definition of $C,\mc{A}$ in this section, we have
\[
Z^* = B_{2e_1^{I},2e_1^{II}} -
\sum_{ (\af,\bt) \in \mathbb{B}^n }
y_{\af,\bt}^* (B_{(\af,\bt)} - s_{(\af,\bt)} B_{(2e_1^{I},2e_1^{II})} )
= \sum_{(\af,\bt) \in \Omega_2^n } y_{\af,\bt}^* B_{\af,\bt}
\]
where $y_{2e_1^{I},2e_1^{II}}^*$ was set to be
$1+\sum_{(\af,\bt) \in \mathbb{B}^n } y_{\af,\bt}^* s_{\af,\bt}$.



If $Z^*$ has rank one, then $f_{min}^{bqp}=f_{sos}^{bqp}$
and a global minimizer for \reff{bq:sphere} can be extracted.
Suppose $Z^*=uu^T$. Define vectors
\[
(x^{I})^* =  \frac{u(1:n_2:n_1n_2)}{\|u(1:n_2:n_1n_2)\|_2}, \quad
(x^{II})^* =  \frac{u(1:n_2)}{\|u(1:n_2)\|_2}, \quad
x^* = ((x^{I})^*, (x^{II})^*).
\]
Obviously $(x^{I})^*, (x^{II})^*$ both have unit length.
The definition of $B_{\af,\bt}$ in \reff{bq:B_afbt}
implies that $Z^*=[x^*]_{bl}[x^*]_{bq}^T$.
Then for any $x=(x^{I},x^{II})$ with $\|x^{I}\|_2 =\|x^{II}\|_2=1$, it holds
\[
-f(x^*) = -f_{2e_1^{I},2e_1^{II}} + b^T y^* \geq -f_{2e_1^{I},2e_1^{II}} +
\sum_{(\af,\bt) \in \mathbb{B}^n} -b_\af (x^{I})^{\af}(x^{II})^{\af}
= -f(x),
\]
because $y^*$ is optimal and $Z=[x]_{bl}[x]_{bl}^T$ is always feasible for \reff{bqsdp:dual}.
So we get
\[
f(x^*) \leq f(x) \quad \forall \, x =(x^{I},x^{II}): \, \|x^{I}\|_2 =\|x^{II}\|_2=1.
\]
Therefore, $x^*$ is a global minimizer of \reff{bq:sphere}
and thus $f_{min}^{bqp}=f_{sos}^{bqp}$.


When $Z^*$ having rank greater than one,
the method using SVDs in \cite[Section~4.1]{LNQY}
can be applied to find approximate solutions from $Z^*$ for \reff{bq:sphere},
but they are not guaranteed to be global minimizers.
In these cases, the SOS relaxation \reff{bq:sos} might not be exact, that is,
$f_{sos}^{bqp} < f_{min}^{bqp}$.
This is not surprising since problem \reff{bq:sphere} is NP-hard.




Now we give some numerical experiments on solving \reff{bq:sphere}
by applying Algorithm~\ref{alg:Newton-CG}.
They are implemented in a Linux machine
of 3.8 GB memory and CPU frequency 2.83 GHz.



\begin{exm}
Consider the following bi-quadratic optimization
\[
\baray{rl}
\underset{x\in\re^n}{\min} &
\underset{ 1\leq i\leq n_1-1 \, } {\sum} \,
\underset{ \,n_1+1 \leq j \leq n } {\sum}
 x_ix_{i+1} x_jx_{j+1} \\
s.t. & \|x^{I}\| = \| x^{II}\| = 1.
\earay
\]
\end{exm}



\begin{exm}[Randomly generated problems]
To test the numerical performance of the regularized SDP method
for solving the bi-quadratic optimization \reff{bq:sphere},
we apply it to solve some randomly generated problem.
To get a random bi-quadratic form $f(x)$ in \reff{bqform},
its coefficients $f_{\af,\bt}$ are generated randomly
with Gaussian distribution.
To avoid numerical difficulty, we normalize $f(x)$
such that every $f_{\af,\bt}$ has absolute value no greater than one.


\begin{table}
\centering
\begin{tabular}{|c|c|r|r|r|r|} \hline
($n_1$,$n_2$) & iter & \multicolumn{3}{c|}{time (min,ave,max)} & gap \\ \hline
(10, 10)  & (6,6,7) & 5.76 &  7.19   &  10.08  & 5.4e-7 \\ \hline
(20, 20)  & (6,12,15)  & 8:35  &  15:04   &  19:15  & 6.6e-6 \\ \hline
(30, 30)  & (13,16,18)  & 1:30:05 & 1:50:53 & 2:04:45  & 1.5e-5\\ \hline
(40, 40)  &   &  &  &  & \\ \hline
(50, 50)  & (,30,) &  & 38:09:24  &  & 5.7e-2 \\ \hline
\end{tabular}
\caption{Computational results for solving randomly generated
bi-quadratic optimization \reff{bq:sphere} by Algorithm~\ref{alg:Newton-CG}.
In the ``iter" column, the triple $(p,q,r)$ means $p,q,r$ respectively denotes
the minimum, median, and maximum outer iteration numbers
required by Algorithm~\ref{alg:Newton-CG}.
The last column is the gap $C\bullet X^* -b^Ty^*$.
}
\label{BQtab:rand}
\end{table}


When $n_1 =n_2 = 10$,  the optimal $Z^*$ for these random problems
almost always has rank one, and we can easily
extract global minimizers by applying the method described earlier.
However, for $n_1 \geq 20, n_2 \geq 20$,
the optimal $Z^*$ for these random problems usually has rank bigger than one.
\end{exm}

\fi
%%%%%%%%







\section{Constrained polynomial optimization}
\setcounter{equation}{0}
\label{sec:conopt}

Consider general constrained polynomial optimization
\be \label{conpop:gi>=0}
\baray{rl}
f_{min}^{con}:=
\underset{x\in \re^n}{\min} & f(x) \\
s.t. & g_1(x)\geq 0, \ldots, g_\ell(x) \geq 0,
\earay
\ee
where $f(x)$ and $g_1(x),\ldots,g_\ell(x)$ are all multivariate polynomials
and have degrees no greater than $2d$.
%A quite broad class of optimization problems can be formulated in the form \reff{conpop:gi>=0}.
%For instance, when $f(x)$ and every $g_i(x)$ are quadratic,
%\reff{conpop:gi>=0} reduces to the quadratically constrained quadratic optimization.
Problem \reff{conpop:gi>=0} is NP-hard, even when
$f(x)$ is quadratic and the feasible set is a simplex.
Lasserre's relaxation is a typical approach for solving \reff{conpop:gi>=0} approximately.
We refer to \cite{DNP,Las01,Las08,Lau09,Sweg05} for the work in this area.



The $d$-th Lasserre's relaxation ($d$ is also called the relaxation order)
for \reff{conpop:gi>=0} is
\be  \label{consos:deg2d}
\baray{rl}
f_{sos}^{con}:=
\max & \gamma \\
s.t. & f(x) - \gamma = \sig_0(x) + g_1(x)\sig_1(x) + \cdots + g_\ell(x)\sig_\ell(x), \\
& \sig_0(x),\sig_1(x),\ldots,  \sig_\ell(x) \mbox{ are SOS}, \\
& \deg(\sig_0),\deg(\sig_1g_1),\ldots,  \deg(\sig_\ell g_\ell) \leq 2d.
\earay
\ee
Let $N(k) = \binom{n+k}{k}$, $d_i=\lceil\deg(g_i)/2\rceil$ and $g_0(x)=1$.
Then $\gamma$ is feasible for \reff{consos:deg2d} if and only if there exists
$X^{(i)} \in\mc{S}^{N(d-d_i)}\,(i=0,1,\ldots,\ell)$ such that
\[
\baray{c}
f(x)-\gamma = \overset{\ell}{\underset{i=0}{\sum}} g_i(x) [x]_{d-d_i}^TX^{(i)}[x]_{d-d_i}
= \overset{\ell}{\underset{i=0}{\sum}} X^{(i)} \bullet (g_i(x)[x]_{d-d_i}[x]_{d-d_i}^T), \\
\, X^{(0)} \succeq 0, X^{(1)} \succeq 0, \ldots, X^{(\ell)} \succeq 0.
\earay
\]
Define constant symmetric matrices
$A_\af^{(0)}, A_\af^{(1)}, \ldots, A_\af^{(\ell)}$ in the way that
\be \label{consos:diagA}
g_i(x)[x]_{d-d_i}[x]_{d-d_i}^T = \sum_{\af\in\N^n: |\af| \leq 2d} A_\af^{(i)} x^\af,
\quad i=0,1,\ldots,\ell.
\ee
Denote $A_\af = (A_\af^{(0)}, A_\af^{(1)}, \ldots, A_\af^{(\ell)})$,
$X = (X^{(0)}, X^{(1)}, \ldots, X^{(\ell)})$, and define
\[
\mc{K} = \mc{S}_+^{N(d-d_0)} \times \mc{S}_+^{N(d-d_1)}
\times \cdots \times \mc{S}_+^{N(d-d_\ell)}.
\]
Recall that $\mathbb{U}_{2d}^n=\{\af\in \N^n: 0 < |\af| \leq 2d\}$.
If $ f(x) =f_0+ \sum_{\af\in\mathbb{U}_{2d}^n} f_\af x^\af$,
then $\gamma$ is feasible for \reff{consos:deg2d}
if and only if there exists $X$ satisfying
\[
\baray{rcl}
A_0 \bullet X + \gamma  & = & f_0, \\
A_\af \bullet X & = & f_\af \quad \forall \, \af\in\mathbb{U}_{2d}^n, \\
X & \in & \mc{K}.
\earay
\]
Now define $\mc{A},b,C$ as
\be  \label{consos:data-def}
\mc{A}(X) = (A_\af \bullet X)_{\af\in\mathbb{U}_{2d}^n}, \quad
b=(f_\af)_{\af\in\mathbb{U}_{2d}^n}, \quad
C=A_0.
\ee
The vector $b$ has dimension $m =N(2d)-1$.
Then SOS relaxation \reff{consos:deg2d} is equivalent to
\be  \label{prmsdp:consos}
\baray{rl}
f_{sdp}^{con}:=
\min &  C \bullet X \\
s.t. & \mc{A}(X) = b, \, X \in \mc{K}.
\earay
\ee
The dual of \reff{prmsdp:consos} is
\be  \label{dualsdp:consos}
\baray{rl}
\max &  b^Ty  \\
s.t. & \mc{A}^*(y) + Z = C,\, Z \in \mc{K}.
\earay
\ee
Here the adjoint $\mc{A}^*(y)$ has the form
\[
\mc{A}^*(y) = \left(
\sum_{\af\in\mathbb{U}_{2d}^n} y_\af A_\af^{(0)},
\sum_{\af\in\mathbb{U}_{2d}^n} y_\af A_\af^{(1)}, \ldots,
\sum_{\af\in\mathbb{U}_{2d}^n} y_\af A_\af^{(\ell)} \right).
\]

\begin{theorem} \label{thm-bd:con}
For notations defined above, it holds
$f_{sos}^{con} = -f_{sdp}^{con}+f_0 \leq f_{min}^{con}$.
If SOS relaxation \reff{consos:deg2d} is exact,
then $f_{min}^{con} = -f_{sdp}^{con}+f_0$.
\end{theorem}
\noindent
The proof is the same as for Theorem~\ref{thm-bd:uc}, and hence is omitted.
A general approximation bound analysis for Lasserre's relaxation is given in \cite{Nie09c}.
\iffalse
% the proof is dropped
\begin{proof}
Optimization \reff{consos:deg2d} and \reff{prmsdp:consos}
are equivalent. Their objectives satisfy
\[
C \bullet X = -\gamma + f_0.
\]
Maximizing $\gamma$ in \reff{consos:deg2d} is equivalent to
minimizing $C\bullet X$ in \reff{prmsdp:consos}.
The feasible sets of \reff{consos:deg2d} and \reff{prmsdp:consos}
are the same. So $f_{sos}^{con} = -f_{\mc{K}}^{con}+f_0$.
The relation $f_{sos}^{con} \leq f_{min}^{con}$ is obvious
because \reff{consos:deg2d} is a relaxation for \reff{conpop:gi>=0}.

When SOS relaxation \reff{consos:deg2d} is exact,
that is, $f_{sos}^{con} = f_{min}^{con}$,
we immediately have
$f_{min}^{con} = -f_{\mc{K}}^{con}+f_0$.
\end{proof}
\fi


Algorithm~\ref{alg:Newton-CG} can be applied to solve
the conic SDP \reff{prmsdp:consos}-\reff{dualsdp:consos}.
Unlike in Sections~\ref{sec:unopt} and \ref{sec:hmgopt},
there is usually no explicit formula for $(\mc{A}\mc{A}^*)^{-1}$ in Lasserre's relaxation,
and it is typically quite difficult to compute $(\mc{A}\mc{A}^*)^{-1}$.
So we just simply use $\big(\diag(\mc{A}\mc{A}^*)\big)^{-1}$ as a preconditioner.
Suppose $(X^*,y^*,Z^*)$ is optimal for \reff{prmsdp:consos}-\reff{dualsdp:consos}.
Then $-f_{sdp}^{con} + f_0$ is a lower bound for
the minimum $f_{min}^{con}$.
The information for minimizers could be obtained from $Z^*$.
Note $Z^*=(Z_0^*,Z_1^*,\ldots,Z_\ell^*)$.
Since $Z^*\in \mc{K}$, every $Z_i^* \succeq 0$.
When $Z_0^*$ has rank one, $f_{min}^{con}=f_{sos}^{con}$
and a global minimizer for \reff{conpop:gi>=0} can be extracted.
Actually, if $\rank(Z_0^*)=1$, \reff{consos:diagA} implies
$Z_0^*=[x^*]_d[x^*]_d^T$ for some $x^*\in \re^n$,
and hence $y_\af^*= -(x^*)^{\af}$ for every $\af \in \mathbb{U}_{2d}^n$.
Furthermore, the structures of $Z_i^*$ implies that
$Z_i^* = g_i(x^*) [x^*]_{d-d_i} [x^*]_{d-d_i}^T$.
Hence, every $g_i(x^*) \geq 0$ and $x^*$ is feasible for \reff{conpop:gi>=0}.
Then, for every $x\in \re^n$ feasible in \reff{conpop:gi>=0}, we have
\[
-f(x^*) = -f_0 + b^T y^* \geq -f_0 + \sum_{\af \in \mathbb{U}_{2d}^n} -b_\af x^{\af}
= -f(x),
\]
because $y^*$ is optimal and for every $x$ feasible in \reff{conpop:gi>=0}
\[
Z=([x]_d[x]_d^T,g_1(x)[x]_{d-d_1}[x]_{d-d_1}^T,\ldots,[x]_{d-d_\ell}[x]_{d-d_\ell}^T )
\]
is feasible for \reff{dualsdp:consos}.
So $x^*$ is a global minimizer of \reff{conpop:gi>=0}.



When $\rank(Z_0^*)>1$, if FEC holds,
several global minimizers could be extracted, e.g., by the method in \cite{HL05}.
Typically, FEC fails if there are a lot of global minimizers
or Lasserre's relaxation \reff{consos:deg2d} is not exact.
If this happens, it is not clear how to get minimizers.





\subsection{Numerical experiments}


Now we present some numerical examples of applying Algorithm~\ref{alg:Newton-CG}
to solve the conic SDP \reff{prmsdp:consos}-\reff{dualsdp:consos}
which is equivalent to Lasserre's relaxation \reff{consos:deg2d}.
As in Section~\ref{subsec:uc-num},
the computations are implemented in Matlab on the same machine.
%a Linux Desktop of 3.8GB memory and 2.8 GHz CPU frequency.
The parameters and stopping criteria of Algorithm~\ref{alg:Newton-CG}
are also the same as there.
The $\big(\diag(\mc{A}\mc{A}^*)\big)^{-1}$ is used as a preconditioner.
%$\eps^{in}=\eps^{out}=10^{-6}$, $\rho=5$, $\sig_0=1$, $\sig_{max} = 10^6$, $\dt = \half$, $K=500$.
%We terminate the inner loop if it takes over $25$ steps,
%and terminate the outer loop if it takes more than $20$ steps.
If FEC holds, we use the method in \cite{HL05} to get one or several global minimizers $x^*$;
otherwise, just set $x^* = Z^*(2:n+1,1)$.
In either case, the relative error of $x^*$ is measured as
\[
err = \frac{|f(x^*)-f_{sos}^{con}|}{\max\{1, |f(x^*)|\} }.
\]
Note we always have $f_{sos}^{con} \leq f_{min}^{con}$.
If the obtained $x^*$ is feasible for \reff{conpop:gi>=0},
then $err = 0$ if and only if $x^*$ is a global minimizer.




\begin{exm} Minimize the polynomial
\[
\sum_{i=1}^n x_i + \sum_{1 \leq i < j \leq n}^n x_ix_j +
 \sum_{1 \leq i < j < k \leq n}^n x_ix_jx_k +
 \sum_{1 \leq i < j < k < l \leq n}^n x_ix_jx_kx_l
\]
over the unit ball $B(0,1)=\{x \in \re^n: \|x\|_2 \leq 1\}$.
It has degree $4$ and is a summation of square free monomials.
We apply the $2^{nd}$ Lasserre's relaxation \reff{consos:deg2d} to solve this problem.
The resulting cone $\mc{K}$ has $2$ blocks.
When $n=50$, solving \reff{prmsdp:consos}-\reff{dualsdp:consos} takes about $5$ hours.
The computed lower bound $f_{sos}^{con} = -1.5240$ (only four decimal digits are shown).
The optimal $Z^*$ has rank $50$ and FEC holds. There are $50$ solutions extracted.
One of them (only four decimal digits are shown) is
\[
( 0.8852, \quad -0.0665, \quad  -0.0665, \quad \ldots, \quad -0.0665, \quad -0.0665  ).
\]
It has unit norm, and hence is feasible. Its relative error is $7\cdot 10^{-7}$.
So it is a global minimizer up to a tiny numerical error.
Actually all its permutations are also solutions of the same accuracy,
since the objective polynomial is symmetric in $x$.
% square free polys over unit ball
%
% n= 50  deg= 4
%
% it    time         gap        lg(rrd)  lg(rrp)   sigma
%  1  03:39:34   -4.90816e+04    1.963  -11.373        1
%  2  03:46:20   -1.39749e+00   -1.056  -10.106        5
%  3  03:55:38   -5.22466e-01   -1.569  -11.971       25
%  4  04:03:01   -3.18366e-02   -2.114  -10.944      125
%  5  04:09:29   -1.03068e-01   -2.514  -10.860      625
%  6  04:17:31   -1.39421e-01   -2.845  -10.158     3125
%  7  04:29:59   -4.16668e-02   -3.599   -8.907    15625
%  8  04:45:28   -9.91632e-03   -4.276  -11.954    15625
%  9  04:52:23   -2.22474e-03   -4.934   -9.340    15625
% 10  05:04:07   -4.95717e-04   -5.588  -11.408    15625
% 11  05:08:43   -1.10313e-04   -6.241  -10.276    15625
%
%ans =
%
%  Columns 1 through 12
%
%   -0.0664   -0.0665   -0.0664   -0.0664   -0.0664   -0.0664   -0.0664   -0.0664   -0.0664   -0.0664   -0.0664   -0.0664
%
%  Columns 13 through 24
%
%   -0.0664   -0.0664   -0.0664    0.8852   -0.0664   -0.0664   -0.0664   -0.0664   -0.0664   -0.0664   -0.0664   -0.0664
%
%  Columns 25 through 36
%
%   -0.0666   -0.0664   -0.0664   -0.0664   -0.0664   -0.0665   -0.0664   -0.0664   -0.0664   -0.0664   -0.0664   -0.0664
%
%  Columns 37 through 48
%
%   -0.0665   -0.0664   -0.0666   -0.0665   -0.0665   -0.0664   -0.0664   -0.0664   -0.0664   -0.0664   -0.0665   -0.0664
%
%  Columns 49 through 50
%
%   -0.0664   -0.0664
%
%
%ans =
%
%   -1.5240   -1.5240
%
%
%error =
%
%  -7.6339e-07
\end{exm}









\begin{exm}
Minimize the polynomial
\[
\sum_{1\leq i < j \leq n} x_i^2x_j^2+ix_i^3x_j + j x_ix_j^3
\]
over the hypercube $[-1,1]^n=\{x \in \re^n:  x_i^2 \leq 1\}$. It has degree $4$.
So we apply the $2^{nd}$ Lasserre's relaxation \reff{consos:deg2d} to solve this problem.
The resulting cone $\mc{K}$ has $n+1$ blocks.
When $n=40$, solving \reff{prmsdp:consos}-\reff{dualsdp:consos} takes about $22$ hours,
and the computed lower bound $f_{sos}^{con}=-2751.5$ (only the first 5 digits are shown).
The optimal $Z^*$ has rank two and FEC holds.
The two extracted solutions (only 4 decimal digits are shown) are:
\[
\baray{l}
\pm (
1,    1,  1,    1,    1,    1,    1,  0.9434, 0.8920, 0.8483,  0.8105,  0.7775,
0.7482,  0.7220,  0.6985,  0.6771, \\
\quad\, 0.6576,  0.6397,  0.6232,  0.6080,0.5938,  0.5806,  0.5682,  0.5566, 0.5457,  0.5354, 0.5257,  \\
\quad\, 0.5165,  0.5077,  -1,   -1,   -1,   -1,   -1,   -1,   -1, -1,   -1,   -1,   -1).
\earay
\]
They are obviously feasible solutions.
The relative error is around $6\cdot 10^{-7}$. So they are global minimizers
up to a small numerical error.
% Minimize the polynomial
% sum_{i<j} x_i^2x_j^2+i*x_i^3x_j + j*x_ix_j^3
% over the hypercube [-1,1]^n
% n= 40  deg= 4
%
% it    time         gap        lg(rrd)  lg(rrp)   sigma
%  1  01:43:22   -1.27569e+06    2.183   -3.320        1
%  2  03:27:02   -1.55792e+02   -0.353   -4.865        5
%  3  05:20:44   -4.57539e+02   -0.415   -4.565       25
%  4  07:14:48   -9.62296e+02   -0.564   -3.613      125
%  5  09:08:22   -6.74244e+02   -0.993   -2.253      625
%  6  11:03:08   -2.33525e+02   -1.610   -1.298     3125
%  7  12:58:04   -7.10337e+01   -2.479   -1.042    15625
%  8  14:52:23    3.80197e+00   -3.090   -0.985    15625
%  9  16:48:05   -7.77686e+00   -3.265   -1.403    15625
% 10  18:44:52   -3.77295e+01   -3.241   -0.894    15625
% 11  20:39:02   -5.71849e-01   -3.239   -5.681    15625
% 12  22:34:09   -1.01934e-03   -6.203   -8.064    15625
%
%ans =
%
%   1.0e+03 *
%
%   -2.7515   -2.7515
%
%
%error =
%
%  -6.7261e-07
%
%there are totally   2 extracted solutions
%ans =
%
%  Columns 1 through 12
%
%   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -0.9434   -0.8920   -0.8483   -0.8105   -0.7775
%
%  Columns 13 through 24
%
%   -0.7482   -0.7220   -0.6985   -0.6771   -0.6576   -0.6397   -0.6232   -0.6080   -0.5938   -0.5806   -0.5682   -0.5566
%
%  Columns 25 through 36
%
%   -0.5457   -0.5354   -0.5257   -0.5165   -0.5077    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000
%
%  Columns 37 through 40
%
%    1.0000    1.0000    1.0000    1.0000
%
%
%ans =
%
%  Columns 1 through 12
%
%    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    0.9434    0.8920    0.8483    0.8105    0.7775
%
%  Columns 13 through 24
%
%    0.7482    0.7220    0.6985    0.6771    0.6576    0.6397    0.6232    0.6080    0.5938    0.5806    0.5682    0.5566
%
%  Columns 25 through 36
%
%    0.5457    0.5354    0.5257    0.5165    0.5077   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000
%
%  Columns 37 through 40
%
%   -1.0000   -1.0000   -1.0000   -1.0000
\end{exm}



\begin{exm}
Minimize the polynomial
\[
\sum_{1\leq i < j \leq n} x_ix_j+x_i^3+x_j^3-x_i^2x_j^2
\]
over the hypercube $[-1,1]^n=\{x \in \re^n:  x_i^2 \leq 1\}$. The degree is $4$.
So we apply the $2^{nd}$ Lasserre's relaxation \reff{consos:deg2d}.
The resulting cone $\mc{K}$ has $n+1$ blocks.
For $n=50$, solving \reff{prmsdp:consos}-\reff{dualsdp:consos} takes about $18$ hours.
The computed lower bound $f_{sos}^{con}=-3600$ (rounded to integers).
The optimal $Z^*$ has rank one, and the solution extracted is (rounded to integers)
\[
\baray{c}
( -1,   -1 ,   1 ,  -1,   -1 ,   1 ,   1  , -1 ,  -1 ,   1,    1,   -1
-1,   -1 ,  -1 ,   1 ,  -1 ,   1 ,   1 ,  -1 ,   1 ,   1 ,  -1,   -1, \\
 1,  -1,   -1,   -1 ,   1 ,  -1,   -1 ,  -1 ,   1 ,   1 ,  -1 ,  -1,
 1,    1,    1 ,  -1 ,  -1 ,   1,    1 ,   1,    1 ,  -1 ,  -1,    1,
1,    1).
\earay
\]
It is clearly a feasible solution.
Its relative error is zero, so it is a global minimizer.
%
% bi-quadrac cubic polys over hypercubes
%
% n= 50  deg= 4
%
% it    time         gap        lg(rrd)  lg(rrp)   sigma
%  1  06:16:53   -2.91343e+05    1.821   -1.172        1
%  2  11:58:07   -1.52829e+03   -0.325   -2.341        5
%  3  15:15:29   -1.71293e+01   -2.130   -5.148       25
%  4  16:50:28   -3.10495e-04   -4.732  -10.639      125
%  5  18:39:02   -5.53510e-06   -6.747   -8.559      625
%
%ans =
%
%  Columns 1 through 12
%
%   -1.0000   -1.0000    1.0000   -1.0000   -1.0000    1.0000    1.0000   -1.0000   -1.0000    1.0000    1.0000   -1.0000
%
%  Columns 13 through 24
%
%   -1.0000   -1.0000   -1.0000    1.0000   -1.0000    1.0000    1.0000   -1.0000    1.0000    1.0000   -1.0000   -1.0000
%
%  Columns 25 through 36
%
%    1.0000   -1.0000   -1.0000   -1.0000    1.0000   -1.0000   -1.0000   -1.0000    1.0000    1.0000   -1.0000   -1.0000
%
%  Columns 37 through 48
%
%    1.0000    1.0000    1.0000   -1.0000   -1.0000    1.0000    1.0000    1.0000    1.0000   -1.0000   -1.0000    1.0000
%
%  Columns 49 through 50
%
%    1.0000    1.0000
%
%
%ans =
%
%   1.0e+03 *
%
%   -3.6000   -3.6000
%
%
%ans =
%
%   4.2151e-05
%
\end{exm}




\begin{exm}
Minimize the polynomial
\begin{align*}
\sum_{i=1}^n ix_i^3 + \sum_{1 \leq i < j \leq n} (i+j) x_i^3x_j^3
\end{align*}
over the unit ball $B(0,1)$. It has degree $6$.
So we apply the $3^{rd}$ Lasserre's relaxation \reff{consos:deg2d}.
The resulting cone $\mc{K}$ has $2$ blocks.
When $n=20$, solving \reff{prmsdp:consos}-\reff{dualsdp:consos} takes about $57$ hours.
The computed lower bound $f_{sos}^{con}=-20$ (rounded to integers).
The optimal $Z^*$ has rank one, and the solution extracted (rounded to integers) is
\[
(0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
 0,   0,   0,   0,   0,   0,   0,  -1).
\]
It is feasible and has zero relative error.
So it is a global minimizer.
%
% minimize sum_{i=1}^n i*x_i^3 + sum_{1 <=i < j <= n}^n (i+j) x_i^3*x_j^3 over unit ball
%
% n= 20  deg= 6
%
% it    time         gap        lg(rrd)  lg(rrp)   sigma
%  1  13:55:32   -4.09393e+04    1.923   -4.825        1
%  2  27:18:48   -1.25262e+01   -0.530   -7.656        5
%  3  34:28:27   -8.06239e+00   -1.077   -9.788       25
%  4  42:12:44   -1.67247e+00   -1.742   -7.786      125
%  5  47:52:15   -2.44128e-01   -2.545  -10.361      625
%  6  54:31:17   -1.01000e-02   -3.925   -8.426     3125
%  7  56:29:33   -8.90688e-05   -5.979   -8.942    15625
%  8  57:11:40   -7.87223e-07   -8.033   -9.127    15625
%
%ans =
%
%  Columns 1 through 12
%
%   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000
%
%  Columns 13 through 20
%
%   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -1.0000
%
%
%ans =
%
%  -20.0000  -20.0000
%
%
%error =
%
%   2.8099e-11
\end{exm}



\begin{exm}
Consider the polynomial optimization
\begin{align*}
\min_{x\in\re^n}  & \quad
\sum_{1 \leq i < j < k \leq n}
(i+j)x_ix_jx_k + (j+k)x_i^2x_j^2x_k^2 \\
s.t. & \quad x_1^4+\cdots+x_n^4 \leq 1.
\end{align*}
The constraint is the unit ball defined by the standard $4$-norm. The degree is $6$.
We apply the $3^{rd}$ Lasserre's relaxation \reff{consos:deg2d} to solve this problem.
The resulting cone $\mc{K}$ has $2$ blocks.
When $n=15$, solving \reff{prmsdp:consos}-\reff{dualsdp:consos} takes about $1.5$ hours.
The computed lower bound $f_{sos}^{con} = -575.5928$ (only four decimal digits are shown).
The optimal $Z^*$ has rank one, and the solution extracted (only four decimal digits are shown) is
\[
\baray{l}
-(0.4034,  0.4274,  0.4486,  0.4674,  0.4839,  0.4983,  0.5107, 0.5211, 0.5296,  0.5363, \\
\quad\,  0.5410, 0.5437, 0.5444,  0.5430,  0.5393).
\earay
\]
It has unit $4$-norm, and hence is feasible.
The relative error is around $4\cdot 10^{-7}$.
So it is a global minimizer up to a small numerical error.
% Minimize the polynomial over the 4-ball
% sum_{i < j < k} (i+j)*x_i*x_j*x_k +(j+k)*x_i^2*x_j^2*x_k^2
% n= 15  deg= 6
%
% it    time         gap        lg(rrd)  lg(rrp)   sigma
%  1  00:32:55   -3.89821e+05    2.413  -10.677        1
%  2  00:42:22   -4.05930e+00   -0.904  -10.946        5
%  3  00:55:18   -8.43494e+00   -0.982  -10.878       25
%  4  01:01:37   -1.36385e+01   -1.210  -10.872      125
%  5  01:10:55   -6.81807e+00   -1.785   -8.767      625
%  6  01:17:25   -6.84172e-01   -2.887  -11.053     3125
%  7  01:25:44   -1.27816e-02   -4.636  -10.888    15625
%  8  01:28:06   -2.29251e-04   -6.384   -9.312    15625
%
%ans =
%
% -575.5925 -575.5928
%
%
%error =
%
%   4.7607e-07
%
%there are totally   1 extracted solutions
%ans =
%
%  Columns 1 through 12
%
%   -0.4034   -0.4274   -0.4486   -0.4674   -0.4839   -0.4983   -0.5107   -0.5211   -0.5296   -0.5363   -0.5410   -0.5437
%
%  Columns 13 through 15
%
%   -0.5444   -0.5430   -0.5393
%
\end{exm}


\begin{exm}
Consider the polynomial optimization
\begin{align*}
\min_{x\in\re^n}  & \quad
\sum_{1 \leq i < j <k \leq n} x_ix_jx_k(1+x_i+x_j+x_k)+
ix_i^6+jx_j^6+kx_k^6 \\
s.t. & \quad x_1^4+\cdots+x_{\frac{n}{2}}^4 \leq 1, \quad
x_{\frac{n}{2}+1}^4+\cdots+x_{n}^4 \leq 1
\end{align*}
where $n$ is an even integer. The degree is $6$,
and the constraint is a bi-unit ball under the $4$-norm.
We apply the $3^{rd}$ Lasserre's relaxation \reff{consos:deg2d} to solve this problem.
The resulting cone $\mc{K}$ has $3$ blocks.
For $n=16$, solving \reff{prmsdp:consos}-\reff{dualsdp:consos} takes about $19$ hours,
and the computed lower bound $f_{sos}^{con}=-1.1078$ (only four decimal digits are shown).
The optimal $Z^*$ has rank one, and the extracted solution (only four decimal digits are shown) is
\[
\baray{l}
-(0.2418,  0.2208,  0.2085,  0.2000,  0.1934,  0.1882,  0.1838,  0.1800, 0.1767,  0.1738, \\
\quad\,  0.1712, 0.1688,  0.1667,  0.1647,  0.1629,  0.1612).
\earay
\]
It is a feasible solution, and the relative error is around $4\cdot 10^{-9}$.
So it is a global minimizer up to a tiny numerical error.
% Minimize the polynomial over the bi 4-ball
%
% sum_{i<j<k} x_i*x_j*x_k(1+x_i+x_j+x_k)+i*x_i^6+j*x_j^ +k*x_k^6
% n= 16  deg= 6
%
% it    time         gap        lg(rrd)  lg(rrp)   sigma
%  1  01:58:59   -1.64941e+07    3.172   -6.206        1
%  2  04:15:54   -8.15341e+00   -0.525   -6.758        5
%  3  06:35:50   -2.36137e+00   -1.129   -6.796       25
%  4  08:43:20   -9.33947e-01   -1.672   -6.744      125
%  5  11:02:37   -5.61002e-01   -2.118   -6.444      625
%  6  13:27:21   -1.58096e-01   -2.815   -6.539     3125
%  7  15:48:58   -1.23524e-02   -3.951   -5.962    15625
%  8  17:53:13   -1.01255e-03   -5.026   -7.462    15625
%  9  19:11:55   -8.56760e-05   -6.089   -8.258    15625
%
%ans =
%
%   -1.1078   -1.1078
%
%
%error =
%
%  -4.7453e-09
%
%there are totally   1 extracted solutions
%ans =
%
%  Columns 1 through 12
%
%   -0.2418   -0.2208   -0.2085   -0.2000   -0.1934   -0.1882   -0.1838   -0.1800   -0.1767   -0.1738   -0.1712   -0.1688
%
%  Columns 13 through 16
%
%   -0.1667   -0.1647   -0.1629   -0.1612
%
%
\end{exm}


\begin{exm}
Consider the polynomial optimization
\begin{align*}
\min_{x\in\re^n}  & \quad
\sum_{1 \leq i < j <k \leq n/2}
ix_ix_jx_k+ jx_{\frac{n}{2}+i}x_{\frac{n}{2}+j}x_{\frac{n}{2}+k} +
kx_ix_jx_kx_{\frac{n}{2}+i}x_{\frac{n}{2}+j}x_{\frac{n}{2}+k} \\
s.t. & \quad x_1^4+\cdots+x_{\frac{n}{2}}^4 \leq 1, \quad
x_{\frac{n}{2}+1}^4+\cdots+x_{n}^4 \leq 1
\end{align*}
where $n$ is even. Since the degree is $6$,
we apply the  $3^{rd}$ Lasserre's relaxation \reff{consos:deg2d}.
The cone $\mc{K}$ has $3$ blocks.
For $n=20$, solving \reff{prmsdp:consos}-\reff{dualsdp:consos} takes about $32$ hours.
The computed lower bound $f_{sos}^{con} =-153.6180$ (only four decimal digits are shown).
The optimal $Z^*$ has rank one, and the solution extracted (only four decimal digits are shown) is
\[
\baray{l}
-(-0.3642,  0.3955,  0.5042,  0.5589,  0.5892,  0.6049,  0.6109,  0.6104,  0.6057,  0.5991,  \\
\qquad 0.5828,  0.5173,  0.5193,  0.5306,  0.5459,  0.5619,  0.5763,  0.5869,  0.5919,  0.5896).
\earay
\]
The first and second subvector both have unit $4$-norm,
so it is a feasible solution. Its relative error is around $5 \cdot 10^{-8}$.
So it is a global minimizer up to a tiny numerical error.
% Minimize the polynomial over the bi 4-ball
% sum_{i<j<k} i*x_i*x_j*x_k+j*x_{n/2+i}*x_{n/2+j}*x_{n/2+k}+ k*x_i*x_j*x_k*x_{n/2+i}*x_{n/2+j}*x_{n/2+k}
% n= 20  deg= 6
%
% it    time         gap        lg(rrd)  lg(rrp)   sigma
%  1  13:48:22   -2.33611e+04    1.748   -2.744        1
%  2  18:04:15   -1.93663e+00   -1.306  -10.060        5
%  3  21:49:48   -8.99972e-01   -1.844   -8.111       25
%  4  26:00:15   -9.84502e-02   -2.934   -8.584      125
%  5  29:38:36   -1.87978e-03   -4.684   -8.107      625
%  6  31:59:22   -6.87451e-06   -7.125   -8.080     3125
%
%ans =
%
% -153.6180 -153.6180
%
%
%error =
%
%   5.6398e-08
%
%there are totally   1 extracted solutions
%ans =
%
%  Columns 1 through 12
%
%    0.3642   -0.3955   -0.5042   -0.5589   -0.5892   -0.6049   -0.6109   -0.6104   -0.6057   -0.5991   -0.5828   -0.5173
%
%  Columns 13 through 20
%
%   -0.5193   -0.5306   -0.5459   -0.5619   -0.5763   -0.5869   -0.5919   -0.5896
%
%
\end{exm}




\subsection{Sparse SOS relaxation for constrained optimization}


In many applications, polynomials in optimization are often sparse.
It is important to exploit sparsity for computational efficiency.
We refer to \cite{KimKW05,KojKW05,KKW09,Las06,ND08,Par03,WKKM}.



Suppose each polynomial $g_j$ only involves a subvector $x_{I_j}=\{x_i: i \in I_j\}$
for an index set $I_j \subseteq \{1,2,\ldots,n\}$.
Consider the sparse polynomial optimization
\be  \label{pop:spar-g}
\baray{rl}
\underset{x\in \re^n }{\min} & f(x) \\
s.t. & g_1(x_{I_1})\geq 0, \ldots, g_\ell(x_{I_\ell}) \geq 0.
\earay
\ee
The sparse generalized Lagrangian function (see Waki et al. \cite{WKKM}) for the above is
\[
L(x,\phi) = f(x) - \overset{\ell}{ \underset{j=1}{\sum} } \phi_j(x_{I_j})g_j(x_{I_j}),
\]
where $\phi_j(x_{I_j})$ are SOS polynomials in $x_{I_j}$. Let $\{C_1,\ldots,C_r\}$
be the set of maximal cliques of a chordal extension of
correlative sparsity pattern (csp) matrix of $L(x,\phi)$ (see Waki et al. \cite{WKKM}).
A sparse Lasserre's relaxation of order $d$ proposed in \cite{WKKM} for \reff{pop:spar-g} is
\be \label{spar-sos}
\baray{rl}
\max &  \gamma \\
s.t. &  f(x) - \gamma =
\overset{r}{ \underset{i=1}{\sum} } \sig_i(x_{C_i})+
\overset{\ell}{ \underset{j=1}{\sum} } \phi_j(x_{I_j})g_j(x_{I_j}), \\
& \deg(\sig_i), \deg(\phi_j g_j) \leq 2d, \,\sig_i, \phi_j  \mbox{ are SOS},\\
&  \, i = 1,\ldots, r, j=1,\ldots,\ell.
\earay
\ee


%For convenience, we still denote $d_i=\lceil\deg(g_i)/2\rceil$, $g_0=1$ and $d_0=0$.
Like the general dense case, \reff{spar-sos} can also be formulated
in the form of conic semidefinite optimization \reff{Ksdp:prm}.
For convenience, we still use the same notations as before.
Note $\gamma$ is feasible for \reff{spar-sos} if and only if there exists
$X^{(i)}\,(i=1,\ldots,r+\ell)$ satisfying
\[
\baray{c}
f(x)-\gamma = \overset{r}{\underset{i=1}{\sum}}  [x_{C_i}]_d^TX^{(i)}[x_{C_i}]_d +
\overset{\ell}{\underset{j=1}{\sum}} g_j(x) [x_{I_j}]_{d-d_j}^TX^{(r+j)}[x_{I_j}]_{d-d_j}, \\
\, X^{(1)} \succeq 0, \ldots, X^{(r)} \succeq 0,
X^{(r+1)} \succeq 0, \ldots, X^{(r+\ell)} \succeq 0.
\earay
\]
Let $\mathbb{W}$ be the set of exponents of nonconstant monomials appearing in \reff{spar-sos}.
For $\af \in \mathbb{W}$, define constant symmetric matrices
$ A_\af^{(1)}, \ldots, A_\af^{(r+\ell)}$ in the way that
\be
\baray{c}
[x_{C_i}]_d [x_{C_i}]_d^T = \underset{\af\in \{0\} \cup \mathbb{W}}{\sum}
 A_\af^{(i)} x^\af, i=1,\ldots, r, \\
g_j(x)[x]_{d-d_j}[x]_{d-d_j}^T = \underset{\af\in \{0\} \cup \mathbb{W}}{\sum}
 A_\af^{(r+j)} x^\af, j=1,\ldots,\ell.
\earay
\ee
It is obvious that
\[
\baray{c}
A_\af^{(i)} \in \mc{S}^{\binom{|C_i|+d}{d}}, i=1,\ldots, r, \quad
A_\af^{(r+j)} \in \mc{S}^{\binom{|I_j|+d-d_j}{d-d_j}}, j=1,\ldots, \ell.
\earay
\]
Denote $A_\af = (A_\af^{(1)}, \ldots, A_\af^{(r+\ell)})$,
$X = (X^{(1)}, \ldots,\ldots, X^{(r+\ell)})$, and define
\[
\mc{K} = \mc{S}_+^{\binom{|C_1|+d}{d}} \times \cdots \times \mc{S}_+^{\binom{|C_r|+d}{d}} \times
\mc{S}_+^{\binom{|I_1|+d-d_1}{d-d_1}} \times \cdots \times \mc{S}_+^{\binom{|I_\ell|+d-d_\ell}{d-d_\ell}}.
\]
If $f(x) = f_0 + \sum_{\af \in \mathbb{W}} f_\af x^\af$,
the constraint of \reff{spar-sos} is equivalent to
\[
\baray{rcl}
A_0 \bullet X + \gamma  & = & f_0, \\
A_\af \bullet X & = & f_\af \quad \forall \, \af \in \mathbb{W}, \\
X & \in & \mc{K}.
\earay
\]
Similar to \reff{consos:data-def}, define $\mc{A},b,C$ as
\be
\mc{A}(X) = (A_\af \bullet X)_{\af \in \mathbb{W}}, \quad
b=(f_\af)_{\af \in \mathbb{W}}, \quad
C=A_0.
\ee
The dimension of $b$ is $m = |\mathbb{W}|$, the cardinality of $\mathbb{W}$.
So \reff{spar-sos} is equivalent to
\be \label{sprLas:SDP:prm}
\baray{rl}
\min &  C \bullet X \\
s.t. & \mc{A}(X) = b, \, X \in \mc{K}.
\earay
\ee
The dual of the above is
\be \label{sprLas:SDP:dul}
\baray{rl}
\max &  b^Ty  \\
s.t. & \mc{A}^*(y) + Z = C,\, Z \in \mc{K}.
\earay
\ee
The adjoint $\mc{A}^*(y)$ now takes the form
\[
\mc{A}^*(y) = \left(
\sum_{\af\in\mathbb{W}} y_\af A_\af^{(1)},
\sum_{\af\in\mathbb{W}} y_\af A_\af^{(2)}, \ldots,
\sum_{\af\in\mathbb{W}} y_\af A_\af^{(r+\ell)} \right).
\]
The information about minimizers would be obtained from $Z^*$.
We refer to Laurent and Mourrain \cite{LauMou09}.

\iffalse
%% the following example is droped since CG does not work well
\begin{exm} Consider the Broyden banded polynomial function \cite{WKKM}
\[
\baray{rl}
\underset{k\in \{1,\ldots,n\} }{\sum} \quad
\left( x_k(2+5x_k^2) + 1 - \sum_{i\in J_k} (1+x_i)x_i   \right)^2
\earay
\]
where $ J_k = \{ j: j \ne k,  \, \max(1,k-5) \leq j \leq \min(n,k+1)\}$.
It is a sparse polynomial. Every summand only involves at most $7$ variables.
We want to minimize it in the whole space.
Since the total degree is $6$, we apply the sparse Lasserre's relaxation \reff{spar-sos} of order $3$.
The csp matrix of this polynomial is banded,
its csp graph is chordal. The set of maximal cliques is
$\{ J_6, J_7, \ldots, J_{n-1} \}$.
So the resulting cone $\mc{K}$ for SDP \reff{sprLas:SDP:prm}-\reff{sprLas:SDP:dul}
coming from sparse Lasserre's relaxation \reff{spar-sos} has totally $n-6$ blocks.
\end{exm}
\fi



\begin{exm} Consider the following sparse polynomial optimization
\[
\baray{rl}
\underset{x\in\re^n}{\min} &
\underset{k\in \{1,\ldots,34\} }{\sum} \quad  \left( \underset{i\in J_k}{\sum} ( x_i +x_i^6)
-\underset{i,j \in J_k: i< j }{\sum}  x_i^2x_j^2 \right) \\
s.t. & 1-\|x_{J_k}\|_4^4 \geq 0, k=1,\ldots,34.
\earay
\]
Here the index set
$
J_k = \{ j: \, \max(1,3k-3) \leq j \leq \min(n,3k+2)\},
$
and $n=100$ is the total number of variables. The degree is $6$.
We apply the sparse Lasserre's relaxation \reff{spar-sos} of order $3$.
Since the csp matrix here is banded, its csp graph is chordal.
So the set of maximal cliques is
$\{J_1,J_2,\ldots, J_{33} \}$ (note $J_{34}\subset J_{33}$).
We apply Algorithm~\ref{alg:Newton-CG} to solve
the resulting conic SDP \reff{sprLas:SDP:prm}-\reff{sprLas:SDP:dul}.
It takes about $15$ minutes,
and the computed lower bound is $-194.8056$ (only four decimal digits are shown).
To get a minimizer, we just simply set $x_i^*=y_{e_i}^*$ for $i=1,\ldots,n$,
and get the solution (only four decimal digits are shown)
\[
-\big(0.7066*2, ( \, 0.6394*3, 0.6385*3 \,)*16,  0.7076*2\big).
%
%-0.7066   -0.7066
%-0.6394   -0.6394   -0.6394
%-0.6385   -0.6385   -0.6385
%-0.6394   -0.6394   -0.6394
%-0.6385   -0.6385   -0.6385
%-0.6394   -0.6394   -0.6394
%-0.6385   -0.6385   -0.6385
%-0.6394   -0.6394   -0.6394
%-0.6385   -0.6385   -0.6385
%-0.6394   -0.6394   -0.6394
%-0.6385   -0.6385   -0.6385
%-0.6394   -0.6394   -0.6394
%-0.6385   -0.6385   -0.6385
%-0.6394   -0.6394   -0.6394
%-0.6385   -0.6385   -0.6385
%-0.6394   -0.6394   -0.6394
%-0.6385   -0.6385   -0.6385
%-0.6394   -0.6394   -0.6394
%-0.6385   -0.6385   -0.6385
%-0.6394   -0.6394   -0.6394
%-0.6385   -0.6385   -0.6385
%-0.6394   -0.6394   -0.6394
%-0.6385   -0.6385   -0.6385
%-0.6394   -0.6394   -0.6394
%-0.6385   -0.6385   -0.6385
%-0.6394   -0.6394   -0.6394
%-0.6385   -0.6385   -0.6385
%-0.6394   -0.6394   -0.6394
%-0.6385   -0.6385   -0.6385
%-0.6394   -0.6394   -0.6394
%-0.6385   -0.6385   -0.6385
%-0.6394   -0.6394   -0.6394
%-0.6385   -0.6385   -0.6385
%-0.7076   -0.7076
\]
In the above, the product $a*t$ means $a$ is repeated for $t$ times in a row.
The above solution has relative error around $8\cdot 10^{-9}$,
so it is a global minimizer up to a tiny numerical error.
% Minimize the polynomial
% sum_{i<j}
% over multi 4-balls via Sparse Lasserre relaxation
% n= 100  deg= 6
%
% it    time         gap        lg(rrd)  lg(rrp)   sigma
%  1  00:04:16   -2.40030e+03    0.724   -9.311        1
%  2  00:05:11   -7.87764e+00   -1.152   -9.279        5
%  3  00:06:04   -1.99874e+00   -1.850   -9.456       25
%  4  00:07:23   -1.81195e-01   -2.871   -8.045      125
%  5  00:09:27   -8.21338e-03   -3.722   -9.354      625
%  6  00:11:51   -1.44066e-03   -4.439   -9.644     3125
%  7  00:14:00   -8.46724e-05   -5.741   -9.456    15625
%  8  00:15:23   -4.38689e-06   -7.038   -9.011    15625
%
%ans =
%
% -194.8056 -194.8056
%
%
%error =
%
%   8.8104e-09
%
%
%ans =
%
%  Columns 1 through 12
%
%   -0.7066   -0.7066   -0.6394   -0.6394   -0.6394   -0.6385   -0.6385   -0.6385   -0.6394   -0.6394   -0.6394   -0.6385
%
%  Columns 13 through 24
%
%   -0.6385   -0.6385   -0.6394   -0.6394   -0.6394   -0.6385   -0.6385   -0.6385   -0.6394   -0.6394   -0.6394   -0.6385
%
%  Columns 25 through 36
%
%   -0.6385   -0.6385   -0.6394   -0.6394   -0.6394   -0.6385   -0.6385   -0.6385   -0.6394   -0.6394   -0.6394   -0.6385
%
%  Columns 37 through 48
%
%   -0.6385   -0.6385   -0.6394   -0.6394   -0.6394   -0.6385   -0.6385   -0.6385   -0.6394   -0.6394   -0.6394   -0.6385
%
%  Columns 49 through 60
%
%   -0.6385   -0.6385   -0.6394   -0.6394   -0.6394   -0.6385   -0.6385   -0.6385   -0.6394   -0.6394   -0.6394   -0.6385
%
%  Columns 61 through 72
%
%   -0.6385   -0.6385   -0.6394   -0.6394   -0.6394   -0.6385   -0.6385   -0.6385   -0.6394   -0.6394   -0.6394   -0.6385
%
%  Columns 73 through 84
%
%   -0.6385   -0.6385   -0.6394   -0.6394   -0.6394   -0.6385   -0.6385   -0.6385   -0.6394   -0.6394   -0.6394   -0.6385
%
%  Columns 85 through 96
%
%   -0.6385   -0.6385   -0.6394   -0.6394   -0.6394   -0.6385   -0.6385   -0.6385   -0.6394   -0.6394   -0.6394   -0.6385
%
%  Columns 97 through 100
%
%   -0.6385   -0.6385   -0.7076   -0.7076
%
%
\end{exm}



\section{Conclusions and some discussions}
\label{sec:cncl-dis}

This paper proposes regularization type methods for solving SOS and Lasserre's relaxations
in large scale polynomial optimization.
After a description of these methods for general conic semidefinite optimization,
we show how to apply them in solving
unconstrained, homogeneous and constrained polynomial optimization.
As demonstrated by extensive numerical experiments, significantly larger problems
would be solved efficiently by them.
For a quartic polynomial optimization problem with $100$ variables,
its SOS relaxation would be solved on a regular computer,
which is almost impossible by applying prior existing SOS solvers.


The proposed regularization methods for solving large scale polynomial optimization
are still in their very early stages.
Much work need to be done for improving their practical performance.
From our experience, the following issues are important.

\bit

\item The performance of Algorithm~\ref{alg:Newton-CG} depends very much on
CG iterations for solving linear system \reff{linsys:Newton}.
In the inner loop, if \reff{linsys:Newton} is solved efficiently by CG,
then Algorithm~\ref{alg:Newton-CG} usually converges fast;
otherwise, a lot of time would be spent on solving \reff{linsys:Newton},
which finally makes Algorithm~\ref{alg:Newton-CG} run very slow.
So the conditioning of \reff{linsys:Newton} seriously affects the performance.
Generally, the accuracy parameters in Algorithm~\ref{alg:Newton-CG}
can not be set too small.

\item For unconstrained and homogeneous polynomial optimization,
there are explicit formulae for the inverse $(\mc{A}\mc{A}^*)^{-1}$
which is used as a preconditioner for solving \reff{linsys:Newton} by CG.
For constrained polynomial optimization,
there is no explicit formula for $(\mc{A}\mc{A}^*)^{-1}$,
and we just use $\big(\diag (\mc{A}\mc{A}^*) \big)^{-1}$ as a preconditioner.
Generally, is there a better preconditioner in polynomial optimization?
This is not clear for us.

\item It is critical to update $\sig_k$ in Algorithm~\ref{alg:Newton-CG}.
%In our numerical experiments, we just simply multiply it by a constant $\rho$,
%which might not always be the best strategy.
From our numerical experience, it is better to choose relatively bigger $\rho$
if \reff{linsys:Newton} is solved accurately by CG;
otherwise, we should choose relatively small $\rho$.
Generally, what is the best strategy for updating $\sig_k$?

\iffalse
%% this might not be very correct
\item Generally Algorithm~\ref{BPM:Ksdp} has slower convergence
than Algorithm~\ref{alg:Newton-CG} does,
since the former is basically a first order method while the latter is a second order method.
However, when linear system \reff{linsys:Newton} is ill-conditioned,
the search direction generated by CG process might work badly,
then Algorithm~\ref{BPM:Ksdp} might be a good choice if a proper steplength is chosen,
because it always generates a descending direction at each step.
\fi

\item Algorithms \ref{BPM:Ksdp} and \ref{alg:Newton-CG} are basically special cases of
Algorithm \ref{regalg:MPRW}.
Is there a more efficient specification of Algorithm \ref{regalg:MPRW}
for solving large scale SOS or Lasserre's relaxations?
Or does polynomial optimization have any special features in applying
Algorithm \ref{regalg:MPRW} or its specifications efficiently?
This is not clear for us.


\eit

The regularization methods in Section~\ref{sec:unopt} are
designed for solving general big conic semidefinite optimization \reff{Ksdp:prm}.
A big advantage of them is the
low memory requirement for storing matrices and vectors during the computation.
Generally, the computer memory could be used to a big extent for storing
problem data and the computation does not require much more extra memory.
Therefore, if a big SDP fits the computer memory,
we can always run the algorithm while possibly waiting for a long time.


\bigskip
\bigskip
\noindent
{\bf Acknowledgement} \quad
The author would like very much to thank Bill Helton and Igor Klep
for fruitful discussions on this work.

\begin{thebibliography}{99}


\bibitem{Blek}
G. Blekherman. There are significantly more nonnegative polynomials than sums of squares.
{\it Israel J. Math.} 153 (2006), 355-380.


\bibitem{BdK02}
I.M. Bomze and E. de Klerk.
Solving standard quadratic optimization problems
via linear, semidefinite and copositive programming.
{\it Journal of Global Optimization}, 24(2), 163-185.

%
%\bibitem{CGT}
%A.R. Conn, N.I.M. Gould and P.L. Toint.
%Testing a class of methods for solving minimization
%problems with simple bounds on the variables.
%{\it Math. Comp.}, 50 (1988), 399-430.


\bibitem{CF00}
R.~Curto and L.~Fialkow.
The truncated complex K-moment problem.
{\it Trans. Am. Math. Soc. 352}, 2825�2855 (2000).

\bibitem{demmel}
J.~Demmel.
{\it Applied numerical linear algebra}.
Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1997.


\bibitem{DNP}
J. Demmel, J. Nie and V. Powers
Representations of positive polynomials on non-compact semialgebraic sets via KKT ideals.
{\it Journal of Pure and Applied Algebra}, Vol. 209, No. 1, pp. 189-200, 2007.


\bibitem{dKP02}
E. de Klerk, and D.V. Pasechnik, D.V.
Approximating of the stability number of a graph via copositive programming.
{\it SIAM Journal on Optimization}, 12(4), 875-892.

%\bibitem{Globa}
%Globallib: \url{http://www.gamsworld.org/global/globallib.htm}.

\bibitem{sdpa}
K. Fujisawa, Y. Futakata, M. Kojima, S. Matsuyama, S. Nakamura, K. Nakata and M. Yamashita.
SDPA-M (SemiDefinite Programming Algorithm in MATLAB),
\url{http://homepage.mac.com/klabtitech/sdpa-homepage/download.html}


\bibitem{HL05}
D.~Henrion and J.~Lasserre.
Detecting global optimality and extracting solutions in GloptiPoly.
{\it Positive polynomials in control}, D. Henrion, A. Garulli Eds.,
{Lecture Notes on Control and Information Sciences},
Vol. 312, Springer, Berlin, 2005, pp. 293�310.

\bibitem{JL05}
D.~Jibetean and M.~Laurent.
Semidefinite approximations for global unconstrained polynomial optimization.
{\it SIAM J. Optim.} 16 (2005), no. 2, 490--514


\bibitem{KimKW05}
S. Kim, M. Kojima and H. Waki.
Generalized lagrangian duals and sums of squares relaxations
of sparse polynomial optimization problems.
{\it SIAM Journal on Optimization},  Vol. 15 (3) 697-719 (2005).

\bibitem{KKW09}
Kim, M. Kojima and H. Waki.
Exploiting sparsity in SDP relaxation for sensor network localization.
{\it SIAM Journal on Optimization},  Vol.~20 (1) 192--215, 2009.


\bibitem{KojKW05}
M. Kojima, S. Kim and H. Waki.
Sparsity in Sums of Squares of Polynomials.
{\it Mathematical Programming}, Vol.103 (1) 45-62 (2005).



\bibitem{Las01}
J. Lasserre. Global optimization with polynomials and the problem of moments.
{\it SIAM J. on Optim.}, Vol. 11, No. 3, 796-817, 2001.

\bibitem{Las06}
J. Lasserre. Convergent SDP-relaxations in polynomial optimization with sparsity.
{\it SIAM Journal on Optimization}, Vol. 17, No. 3, 822-843, 2006.

	
\bibitem{Las06b}
J. Lasserre.
A sum of squares approximation of nonnegative polynomials.
{\it SIAM Journal on Optimization}, Vol.~16, No. 3, pp. 751-765.

\bibitem{Las08}
J.~Lasserre.
A semidefinite programming approach to the generalized problem of moments.
{\it Math. Program.} 112 (2008), no. 1, Ser. B, 65--92.


\bibitem{Lau09}
M.~Laurent.
Sums of squares, moment matrices and optimization over polynomials.
{\it Emerging Applications of Algebraic Geometry,
IMA Volumes in Mathematics and its Applications}, Vol. 149,
M. Putinar and S. Sullivant (eds.),
Springer, pp. 157--270, 2009.


\bibitem{LauMou09}
M.~Laurent and B.~Mourrain.
A generalized flat extension theorem for moment matrices.
{\it Archiv der Mathematik}, 93(1), 2009.


%
%\bibitem{LNQY}
%C. Ling, J. Nie, L. Qi, and Y. Ye.
%Bi-quadratic optimization over unit spheres and semidefinite programming relaxations.
%{\it Preprint}, 2008.


\bibitem{MPRW}
J. Malick, J. Povh, F. Rendl and A. Wiegele.
Regularization methods for semidefinite programming.
%Technical Report, Alpen-Adria-Universit\"{a}t Klagenfurt, Austria, October, 2007.
{\it SIAM Journal on Optimization},
Vol.~20, No.~1, pp. 336--356, 2009.


%\bibitem{MGH}
%J. Mor\'{e}, B. Garbow and K. Hillstrom.
%Testing unconstrained optimization software.
%{\it ACM Trans. Math. Soft.}, 7 (1981), 17-41.

\bibitem{MS65}
T.S. Motzkin and E.G. Straus.
Maxima for graphs and a new proof of a theorem of Tur\'{a}n.
{\it Canadian J. Math.}, 17 (1965), pp. 533�540.

\bibitem{Nes00}
Y.~Nesterov.
Squared functional systems and optimization problems.
{\it High Performance Optimization} H.Frenk et al. (eds.), Kluwer Academic Publishers,
2000, pp. 405�440.

\bibitem{Nes03}
Y.~Nesterov.
Random walk in a simplex and quadratic optimization over convex polytopes.
{\it CORE Discussion Paper}, UCL, Louvain-la-Neuve, Belgium, 2003.

\bibitem{NDS}
J.~Nie, J.~Demmel and B.~Sturmfels.
Minimizing polynomials via sum of squares over the gradient ideal.
{\it Math. Program.} 106 (2006), no. 3, Ser. A, 587--606.

\bibitem{Nie09b}
J. Nie.
Sum of squares methods for minimizing polynomial forms over spheres and hypersurfaces.
{\it Preprint}, 2009. \url{http://www.math.ucsd.edu/~njw}

\bibitem{Nie09c}
J. Nie.
An approximation bound analysis for Lasserre's relaxation in multivariate polynomial optimization.
{\it Preprint}, 2009. \url{http://www.math.ucsd.edu/~njw}


\bibitem{ND08}
J.~Nie and J.~Demmel.
Sparse SOS relaxations for minimizing functions that are summations of small polynomials.
{\it SIAM Journal On Optimization}, Vol.~19, No.~4, pp. 1534-1558 (2008).


\bibitem{ParMP}
P. Parrilo. Semidefinite programming relaxations for semialgebraic problems.
{\it Math. Program.}, Ser. B 96 (2003), No. 2, 293-320.


\bibitem{PS01}
P. Parrilo and B. Sturmfels. Minimizing polynomial functions.
{\it Proceedings of the DIMACS
Workshop on Algorithmic and Quantitative Aspects of Real Algebraic Geometry in Mathe-
matics and Computer Science (March 2001)} (eds. S. Basu and L. Gonzalez-Vega), American
Mathematical Society, 2003, pp. 83-100.

\bibitem{Par03}
P.A. Parrilo. Exploiting structure in sum of squares programs.
{\it Proceedings for the 42nd IEEE Conference on Decision and Control}, Maui, Hawaii, 2003.

\bibitem{PRW06}
J. Povh, F. Rendl, and A. Wiegele.
A boundary point method to solve semidefinite programs.
{\it Computing} 78, 277-286 (2006).


\bibitem{Roc76CO}
R.T.~Rockafellar.
Monotone operators and the proximal point algorithm.
{\it SIAM J. Control and Optim.}, 14 (1976), No.5, pp. 877-898.


\bibitem{Roc76MOR}
R.T.~Rockafellar.
Augmented Lagrangians and applications of the proximal point algorithm in convex programming.
{\it Math. Oper. Res. }, 1 (1976), no. 2, 97--116.


\bibitem{sedumi}
J.F. Sturm. SeDuMi 1.02 � a MATLAB toolbox for optimization over symmetric cones.
{\it Optimization Methods and Software}, 11 \& 12 (1999), pp. 625--653.

\bibitem{Sweg05}
M.~Schweighofer.
Optimization of polynomials on compact semialgebraic sets.
{\it SIAM J. Optim.} 15 (2005), no. 3, 805--825.

\bibitem{Tod01}
M. Todd. Semidefinite Optimization.
{\it Acta Numerica} 10 (2001), pp. 515--560.

\bibitem{sdpt3}
K.C. Toh, M.J. Todd, and R.H. Tutuncu.
SDPT3 � a Matlab software package for semidefinite programming.
{\it Optimization Methods and Software}, 11 (1999), pp. 545--581.


%
%\bibitem{Toi78}
%Ph. L. Toint.
%On sparse and symmetric matrix updating subject to a linear equation.
%{\it Math. Comp.} 31 (1977), no 140, 954--961.


\bibitem{VB96}
L. Vandenberghe and S. Boyd. Semidefinite Programming.
{\it SIAM Review} 38 (1996), pp. 49--95.



\bibitem{WKKM}
H. Waki, S. Kim, M. Kojima and M. Muramatsu.
Sums of squares and semidefinite programming relaxations
for polynomial optimization problems with Structured Sparsity.
{\it SIAM Journal on Optimization}, Vol.17 (1) 218-242 (2006).


\bibitem{WSV}
H. Wolkowicz, R. Saigal, and L. Vandenberghe, editors.
{\it Handbook of semidefinite programming}, Kluwer�s Publisher, 2000.



\bibitem{ZST08}
X.Y. Zhao, D.F. Sun, and K.C. Toh.
A Newton-CG Augmented Lagrangian method for semidefinite programming.
{\it Preprint}, National University of Singapore, March 2008.


\end{thebibliography}





\end{document}

0909.3551.extracted_bib
\begin{thebibliography}{99}


\bibitem{Blek}
G. Blekherman. There are significantly more nonnegative polynomials than sums of squares.
{\it Israel J. Math.} 153 (2006), 355-380.


\bibitem{BdK02}
I.M. Bomze and E. de Klerk.
Solving standard quadratic optimization problems
via linear, semidefinite and copositive programming.
{\it Journal of Global Optimization}, 24(2), 163-185.

%
%\bibitem{CGT}
%A.R. Conn, N.I.M. Gould and P.L. Toint.
%Testing a class of methods for solving minimization
%problems with simple bounds on the variables.
%{\it Math. Comp.}, 50 (1988), 399-430.


\bibitem{CF00}
R.~Curto and L.~Fialkow.
The truncated complex K-moment problem.
{\it Trans. Am. Math. Soc. 352}, 2825�2855 (2000).

\bibitem{demmel}
J.~Demmel.
{\it Applied numerical linear algebra}.
Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1997.


\bibitem{DNP}
J. Demmel, J. Nie and V. Powers
Representations of positive polynomials on non-compact semialgebraic sets via KKT ideals.
{\it Journal of Pure and Applied Algebra}, Vol. 209, No. 1, pp. 189-200, 2007.


\bibitem{dKP02}
E. de Klerk, and D.V. Pasechnik, D.V.
Approximating of the stability number of a graph via copositive programming.
{\it SIAM Journal on Optimization}, 12(4), 875-892.

%\bibitem{Globa}
%Globallib: \url{http://www.gamsworld.org/global/globallib.htm}.

\bibitem{sdpa}
K. Fujisawa, Y. Futakata, M. Kojima, S. Matsuyama, S. Nakamura, K. Nakata and M. Yamashita.
SDPA-M (SemiDefinite Programming Algorithm in MATLAB),
\url{http://homepage.mac.com/klabtitech/sdpa-homepage/download.html}


\bibitem{HL05}
D.~Henrion and J.~Lasserre.
Detecting global optimality and extracting solutions in GloptiPoly.
{\it Positive polynomials in control}, D. Henrion, A. Garulli Eds.,
{Lecture Notes on Control and Information Sciences},
Vol. 312, Springer, Berlin, 2005, pp. 293�310.

\bibitem{JL05}
D.~Jibetean and M.~Laurent.
Semidefinite approximations for global unconstrained polynomial optimization.
{\it SIAM J. Optim.} 16 (2005), no. 2, 490--514


\bibitem{KimKW05}
S. Kim, M. Kojima and H. Waki.
Generalized lagrangian duals and sums of squares relaxations
of sparse polynomial optimization problems.
{\it SIAM Journal on Optimization},  Vol. 15 (3) 697-719 (2005).

\bibitem{KKW09}
Kim, M. Kojima and H. Waki.
Exploiting sparsity in SDP relaxation for sensor network localization.
{\it SIAM Journal on Optimization},  Vol.~20 (1) 192--215, 2009.


\bibitem{KojKW05}
M. Kojima, S. Kim and H. Waki.
Sparsity in Sums of Squares of Polynomials.
{\it Mathematical Programming}, Vol.103 (1) 45-62 (2005).



\bibitem{Las01}
J. Lasserre. Global optimization with polynomials and the problem of moments.
{\it SIAM J. on Optim.}, Vol. 11, No. 3, 796-817, 2001.

\bibitem{Las06}
J. Lasserre. Convergent SDP-relaxations in polynomial optimization with sparsity.
{\it SIAM Journal on Optimization}, Vol. 17, No. 3, 822-843, 2006.

	
\bibitem{Las06b}
J. Lasserre.
A sum of squares approximation of nonnegative polynomials.
{\it SIAM Journal on Optimization}, Vol.~16, No. 3, pp. 751-765.

\bibitem{Las08}
J.~Lasserre.
A semidefinite programming approach to the generalized problem of moments.
{\it Math. Program.} 112 (2008), no. 1, Ser. B, 65--92.


\bibitem{Lau09}
M.~Laurent.
Sums of squares, moment matrices and optimization over polynomials.
{\it Emerging Applications of Algebraic Geometry,
IMA Volumes in Mathematics and its Applications}, Vol. 149,
M. Putinar and S. Sullivant (eds.),
Springer, pp. 157--270, 2009.


\bibitem{LauMou09}
M.~Laurent and B.~Mourrain.
A generalized flat extension theorem for moment matrices.
{\it Archiv der Mathematik}, 93(1), 2009.


%
%\bibitem{LNQY}
%C. Ling, J. Nie, L. Qi, and Y. Ye.
%Bi-quadratic optimization over unit spheres and semidefinite programming relaxations.
%{\it Preprint}, 2008.


\bibitem{MPRW}
J. Malick, J. Povh, F. Rendl and A. Wiegele.
Regularization methods for semidefinite programming.
%Technical Report, Alpen-Adria-Universit\"{a}t Klagenfurt, Austria, October, 2007.
{\it SIAM Journal on Optimization},
Vol.~20, No.~1, pp. 336--356, 2009.


%\bibitem{MGH}
%J. Mor\'{e}, B. Garbow and K. Hillstrom.
%Testing unconstrained optimization software.
%{\it ACM Trans. Math. Soft.}, 7 (1981), 17-41.

\bibitem{MS65}
T.S. Motzkin and E.G. Straus.
Maxima for graphs and a new proof of a theorem of Tur\'{a}n.
{\it Canadian J. Math.}, 17 (1965), pp. 533�540.

\bibitem{Nes00}
Y.~Nesterov.
Squared functional systems and optimization problems.
{\it High Performance Optimization} H.Frenk et al. (eds.), Kluwer Academic Publishers,
2000, pp. 405�440.

\bibitem{Nes03}
Y.~Nesterov.
Random walk in a simplex and quadratic optimization over convex polytopes.
{\it CORE Discussion Paper}, UCL, Louvain-la-Neuve, Belgium, 2003.

\bibitem{NDS}
J.~Nie, J.~Demmel and B.~Sturmfels.
Minimizing polynomials via sum of squares over the gradient ideal.
{\it Math. Program.} 106 (2006), no. 3, Ser. A, 587--606.

\bibitem{Nie09b}
J. Nie.
Sum of squares methods for minimizing polynomial forms over spheres and hypersurfaces.
{\it Preprint}, 2009. \url{http://www.math.ucsd.edu/~njw}

\bibitem{Nie09c}
J. Nie.
An approximation bound analysis for Lasserre's relaxation in multivariate polynomial optimization.
{\it Preprint}, 2009. \url{http://www.math.ucsd.edu/~njw}


\bibitem{ND08}
J.~Nie and J.~Demmel.
Sparse SOS relaxations for minimizing functions that are summations of small polynomials.
{\it SIAM Journal On Optimization}, Vol.~19, No.~4, pp. 1534-1558 (2008).


\bibitem{ParMP}
P. Parrilo. Semidefinite programming relaxations for semialgebraic problems.
{\it Math. Program.}, Ser. B 96 (2003), No. 2, 293-320.


\bibitem{PS01}
P. Parrilo and B. Sturmfels. Minimizing polynomial functions.
{\it Proceedings of the DIMACS
Workshop on Algorithmic and Quantitative Aspects of Real Algebraic Geometry in Mathe-
matics and Computer Science (March 2001)} (eds. S. Basu and L. Gonzalez-Vega), American
Mathematical Society, 2003, pp. 83-100.

\bibitem{Par03}
P.A. Parrilo. Exploiting structure in sum of squares programs.
{\it Proceedings for the 42nd IEEE Conference on Decision and Control}, Maui, Hawaii, 2003.

\bibitem{PRW06}
J. Povh, F. Rendl, and A. Wiegele.
A boundary point method to solve semidefinite programs.
{\it Computing} 78, 277-286 (2006).


\bibitem{Roc76CO}
R.T.~Rockafellar.
Monotone operators and the proximal point algorithm.
{\it SIAM J. Control and Optim.}, 14 (1976), No.5, pp. 877-898.


\bibitem{Roc76MOR}
R.T.~Rockafellar.
Augmented Lagrangians and applications of the proximal point algorithm in convex programming.
{\it Math. Oper. Res. }, 1 (1976), no. 2, 97--116.


\bibitem{sedumi}
J.F. Sturm. SeDuMi 1.02 � a MATLAB toolbox for optimization over symmetric cones.
{\it Optimization Methods and Software}, 11 \& 12 (1999), pp. 625--653.

\bibitem{Sweg05}
M.~Schweighofer.
Optimization of polynomials on compact semialgebraic sets.
{\it SIAM J. Optim.} 15 (2005), no. 3, 805--825.

\bibitem{Tod01}
M. Todd. Semidefinite Optimization.
{\it Acta Numerica} 10 (2001), pp. 515--560.

\bibitem{sdpt3}
K.C. Toh, M.J. Todd, and R.H. Tutuncu.
SDPT3 � a Matlab software package for semidefinite programming.
{\it Optimization Methods and Software}, 11 (1999), pp. 545--581.


%
%\bibitem{Toi78}
%Ph. L. Toint.
%On sparse and symmetric matrix updating subject to a linear equation.
%{\it Math. Comp.} 31 (1977), no 140, 954--961.


\bibitem{VB96}
L. Vandenberghe and S. Boyd. Semidefinite Programming.
{\it SIAM Review} 38 (1996), pp. 49--95.



\bibitem{WKKM}
H. Waki, S. Kim, M. Kojima and M. Muramatsu.
Sums of squares and semidefinite programming relaxations
for polynomial optimization problems with Structured Sparsity.
{\it SIAM Journal on Optimization}, Vol.17 (1) 218-242 (2006).


\bibitem{WSV}
H. Wolkowicz, R. Saigal, and L. Vandenberghe, editors.
{\it Handbook of semidefinite programming}, Kluwer�s Publisher, 2000.



\bibitem{ZST08}
X.Y. Zhao, D.F. Sun, and K.C. Toh.
A Newton-CG Augmented Lagrangian method for semidefinite programming.
{\it Preprint}, National University of Singapore, March 2008.


\end{thebibliography}