# Extended Linear Stiefel Manifold Optimization ††thanks: This research was
supported by the National Natural Science Foundation of China (Grant Nos.
12171021 and 11822103), and the Fundamental Research Funds for the Central
Universities.

Mengmeng Song LMIB of the Ministry of Education, School of Mathematical
Sciences, Beihang University, Beijing 100191, People’s Republic of China,
{songmengmeng<EMAIL_ADDRESS>Yong Xia22footnotemark: 2 Corresponding
author.

###### Abstract

The extended linear Stiefel Manifold optimization (ELS) is studied in the
first time. It aims at minimizing a linear objective function over the set of
all $p$-tuples of orthonormal vectors in ${\mathbb{R}}^{n}$ satisfying $k$
additional linear constraints. Despite the classical polynomial-time solvable
case $k=0$, (ELS) is in general NP-hard. According to the Shapiro-Barvinok-
Pataki theorem, (ELS) admits an exact semidefinite programming (SDP)
relaxation when $p(p+1)/2\leq n-k$, which is tight when $p=1$. Surprisingly,
we can greatly strengthen this sufficient exactness condition to be $p\leq
n-k$, which covers the extreme case $n-p=0=k$. Regarding (ELS) as a smooth
nonlinear programming problem, we reveal a nice property that under LICQ, the
standard first- and second-order local necessary optimality conditions are
sufficient for global optimality when $p+1\leq n-k$.

AMS Subject Classification:  90C26, 90C22, 90C46, 90C20

##  1 Introduction

The Stiefel manifold [29] is the set of all $p$-tuples of orthonormal vectors
in ${\mathbb{R}}^{n}$. We denote it by

$St_{n,p}=\\{X\in{\mathbb{R}}^{n\times p}:~{}X^{T}X=I_{p}\\},$

where $I_{p}$ is the identity matrix of order $p\in\\{1,\cdots,n\\}$. The two
extreme cases $St_{n,1}$ and $St_{n,n}$ correspond to the unit sphere and
“orthogonal group”, respectively. $St_{n,p}$ is connected if $p<n$, and has
two connected components if $p=n$ [24].

The fundamental linear Stiefel manifold optimization [1] reads as

$\min_{X\in St_{n,p}}{\rm tr}(A_{0}X),$ (LS)

where $A_{0}\in{\mathbb{R}}^{p\times n}$ and ${\rm tr}(\cdot)$ is the trace of
matrix $(\cdot)$. It can be globally solved by the singular value
decomposition [33]. In fact, the optimal value is equal to minus the sum of
all singular values of $A_{0}$. Moreover, (LS) is equivalent to the
minimization over the convex hull of $St_{n,p}$, which can be exactly
characterized by linear matrix inequalities [35]. Problem (LS) has an
application in the projection onto $St_{n,p}$ [20, 25], more precisely,

$\arg\min_{X\in St_{n,p}}\|A_{0}^{T}-X\|_{F}^{2}=\arg\min_{X\in St_{n,p}}{\rm
tr}(-A_{0}X),$ (1.1)

where $\|\cdot\|_{F}$ is the Frobenius norm.

In this paper, we consider the following extended linear Stiefel manifold
optimization:

$\begin{array}[]{cl}\min&{\rm tr}(A_{0}X)\\\ {\rm s.t.}&X\in
E_{n,p,k}\triangleq\\{X\in St_{n,p}:~{}a_{i}^{l}\leq{\rm tr}(A_{i}X)\leq
a_{i}^{u},~{}i=1,\cdots,k\\},\end{array}$ (ELS)

where $A_{0}\in{\mathbb{R}}^{p\times n}$, $A_{i}\in{\mathbb{R}}^{p\times n}$,
$a_{i}^{l},~{}a_{i}^{u}\in{\mathbb{R}}$, for $i=1,\cdots,k$, and we assume
$E_{n,p,k}\neq\emptyset$. To the best of our knowledge, there are no studies
aimed at the new model (ELS), except for the extreme case $k=0$, which
corresponds to (LS). We summarize plentiful applications of (ELS) in Section
2, where the first one implies the NP-hardness of (ELS).

Based on quadratic equalities representation of $St_{n,p}$, (ELS) is
formulated as a quadratic matrix programming (QMP) problem [4]. By using the
classical Shapiro-Barvinok-Pataki theorem [3, 22, 26], Beck [6] showed that
QMP with a limited number of constraints admits an exact semidefinite
programming (SDP) relaxation. Moreover, in that case, the globally optimal
solution of QMP can be extracted from that of the SDP relaxation. Problem with
this property is called hidden convex optimization [34]. As a direct
application, one can reveal the hidden convexity of (ELS) under the condition
$p(p+1)/2\leq n-k$. This condition is tight for the case of $p=1$, but fails
to hold true when $p\geq\sqrt{2n}$.

To our surprise, by carefully analyzing the structure of $E_{n,p,k}$, we can
greatly strengthen the condition to be $p\leq n-k$. The new result implies the
hidden convexity of (LS) by setting $p=n$ and $k=0$.

For hidden convex optimization without any local nonglobal minimizer, it could
be expected to characterize global optimality based on local information. It
is known that the first- and second-order local necessary optimality
conditions are sufficient for global optimality of the homogeneous (TRS) [24],
the problem of minimizing the sum of heterogeneous quadratic forms on
$St_{n,p}$ [23, 24], QMP problems arising from high-rank Burer-Monteiro
factorizations [10] of SDP problems [9, 30], and so on. For hidden convex
optimization with at least one local nonglobal minimizer, the local optimality
conditions are characterized case by case. The nonconvex trust-region
subproblem has at most one local nonglobal minimizer [21], at which the first-
and second-order sufficient optimality conditions are necessary [31]. As an
extension, this property holds at any strict local nonglobal minimizer of the
homogeneous quadratic optimization with two quadratic constraints [27].
Another line of extension is considered in [28].

Our second contribution in this paper is to show that under the linear
independence constraint qualification (LICQ), the first- and second-order
local necessary optimality conditions are sufficient for global optimality of
(ELS) when $k\leq n-p-1$.

The remainder of this paper is organized as follows. Section 2 gives some
applications of (ELS) and then shows the NP-hardness of (ELS). Section 3
presents an SDP relaxation of (ELS) and the exactness inheriting from the
result about QMP. Section 4 establishes a much tighter bound on the exactness
of the same SDP relaxation. Section 5 characterizes global optimality by the
first- and second-order necessary local optimality conditions. Conclusions are
made in Section 6.
Notation. For an optimization problem $(P)$, denote the feasible region and
the optimal value of $(P)$ by $\mathcal{F}(P)$ and $v(P)$, respectively. The
extreme points of a given set $\Omega$ form the set $\mathcal{E}(\Omega)$.
Denote by ${\mathbb{S}}^{n}$ the set of symmetric matrices of dimension
$n\times n$. For $A\in{\mathbb{S}}^{n}$, $A\succeq(\preceq)\ 0$ means that $A$
is positive (negative) semidefinite. For matrix $X\in{\mathbb{R}}^{n\times
p}$, $X_{[i_{1}:i_{2}][j_{1}:j_{2}]}$ is the submatrix formed by the entries
in both $\\{i_{1},\cdots,i_{2}\\}$-th rows and $\\{j_{1},\cdots,j_{2}\\}$-th
columns. Notation $X_{[i_{1}:i_{2}]}$ is short for
$X_{[i_{1}:i_{2}][i_{1}:i_{2}]}$. Denote by ${\rm vec}(X)$ the vector obtained
by stacking the columns of $X$ one underneath the other. The Kronecker product
of matrices $A$ and $B$ is denoted by $A\otimes B$. Denote by $0_{n}$ the zero
matrix of dimension $n\times n$. Let $E^{ij}$ be the square matrix of a proper
dimension with all entries being $0$ except the $i$-th row $j$-th column entry
being $1$. Let $e_{i}$ be the $i$-th column of the identity matrix $I_{n}$.
For $\lambda\in{\mathbb{R}}$, define $\lambda^{+}:={\rm max}\\{0,\lambda\\}$
and $\lambda^{-}:={\rm max}\\{0,-\lambda\\}$.

##  2 Applications

We present some applications of the new model (ELS) and show that (ELS) is in
general NP-hard.

######  Example 2.1 (Binary linear programming).

Let $A\in{\mathbb{R}}^{m\times n}$, $a\in{\mathbb{R}}^{n}$ and
$b\in{\mathbb{R}}^{m}$. The binary linear programming problem

$\min\\{a^{T}x:~{}Ax\leq b,~{}x\in\\{-1,1\\}^{n}\\}$

is equivalent to the following case of (ELS)

$\min_{X\in St_{n,1}}\left\\{\sqrt{n}a^{T}X:~{}\sqrt{n}AX\leq b,\
-\frac{1}{\sqrt{n}}\leq X_{i}\leq\frac{1}{\sqrt{n}},~{}i=1,\cdots,n\right\\},$

since it holds that

$\\{-1,1\\}^{n}=\left\\{\sqrt{n}X:~{}-\frac{1}{\sqrt{n}}\leq
X_{i}\leq\frac{1}{\sqrt{n}},~{}i=1,\cdots,n,~{}X\in St_{n,1}\right\\}.$ (2.1)

######  Example 2.2 (Hadamard conjecture).

A Hadamard matrix is a square matrix whose entries are either $+1$ or $-1$ and
whose rows are mutually orthogonal. The long-standing Hadamard conjecture is
on the existence of the Hadamard matrix of any order $n$ divisible by $4$, see
[37]. It amounts to the solvability of the following case of (ELS):

$\min_{X\in St_{n,n}}\left\\{{\rm tr}(0_{n}\cdot X):~{}-\frac{1}{\sqrt{n}}\leq
X_{ij}\leq\frac{1}{\sqrt{n}},~{}i,j=1,\cdots,n\right\\},$

by a reformulation similar to (2.1).

######  Example 2.3 (Linear sum assignment problem).

Let $e=(1,\cdots,1)^{T}\in{\mathbb{R}}^{n}$. The linear sum assignment problem
[12]

$\displaystyle\min$ $\displaystyle{\rm tr}(A_{0}X)$ (2.2) $\displaystyle{\rm
s.t.}$ $\displaystyle X\in\Pi_{n}\triangleq\\{X\in\\{0,1\\}^{n\times
n}:~{}Xe=X^{T}e=e\\}$

can be reformulated in (ELS):

$\displaystyle\min\limits_{X\in St_{n,n}}$ $\displaystyle{\rm tr}(A_{0}X)$
(2.3) $\displaystyle{\rm s.t.}$ $\displaystyle X_{ij}\geq
0,~{}i,j=1,\cdots,n.$

It is further equivalent to the linear programming

$\min\\{{\rm tr}(A_{0}X):~{}Xe=X^{T}e=e,~{}X_{ij}\geq 0,~{}i,j=1,\cdots,n\\},$

by the well-known Birkhoff-von Neumann theorem [8].

######  Example 2.4 (Check optimality in the extended trust-region
subproblem).

Let $Q\in{\mathbb{S}}^{n}$ and $c\in{\mathbb{R}}^{n}$. The nonconvex trust-
region subproblem (TRS) [38] reads as

$\min\\{x^{T}Qx+2c^{T}x:~{}x^{T}x\leq 1\\},$ (TRS)

where $Q\not\succeq 0$. (TRS) has various convex reformulations [34]. The set
of all global minimizers of (TRS), denoted by $X^{*}$, is either a singleton
or a sphere with dimension $s$ for some $s\leq n$. Consider the extended (TRS)
with additional $m$ linear constraints. It is polynomial-time solvable when
$m$ is fixed [7, 16]. The first step in the global optimization algorithm [16]
is to solve (TRS) and then check whether a system of $m$ linear constraints
has a solution in $X^{*}$, which corresponds to a special case of (ELS) with
$p=1$ and $A_{0}=0$, see [16, Section 2.2].

######  Example 2.5 (Unbalanced orthogonal Procrustes problem).

Let $A\in{\mathbb{R}}^{m\times n}$ and $B\in{\mathbb{R}}^{m\times p}$. Solving
the linear matrix equation $AX=B$ on the Stiefel manifold $St_{n,p}$ belongs
to (ELS) with $A_{0}=0$. The corresponding least squares problem is known as
the unbalanced orthogonal Procrustes problem [14, 35]. It reduces to the
quadratic constrained least squares problem [15] when $p=1$.

######  Example 2.6 (Minimax Stiefel manifold optimization).

For $A_{i}\in{\mathbb{R}}^{p\times n}$ and $c_{i}\in{\mathbb{R}}$
($i=1,\cdots,m$), consider the minimax problem

$\min_{X\in E_{n,p,k}}\max_{i=1,\cdots,m}\\{{\rm tr}(A_{1}X)+c_{1},\cdots,{\rm
tr}(A_{m}X)+c_{m}\\}.$ (M-M)

When $k=0$, (M-M) is fundamental in the nonsmooth Stiefel manifold
optimization [13, 39]. Problem (M-M) is equivalent to the minimum of the
following $m$ instances of (ELS) for $i=1,\cdots,m$:

$\min\limits_{X\in E_{n,p,k}}\\{{\rm tr}(A_{i}X):~{}{\rm
tr}(A_{i}X)+c_{i}\geq{\rm tr}(A_{j}X)+c_{j},~{}\forall j\neq i\\}.$

We list some applications of (M-M) as follows.

  * •

The sphere-constrained weighted maximin dispersion problem [32] reads as

$\max_{X\in St_{n,1}}\min_{i=1,\cdots,m}w_{i}\|X-d_{i}\|^{2}=-\min_{X\in
St_{n,1}}\max_{i=1,\cdots,m}-w_{i}\|X-d_{i}\|^{2},$

where $\|\cdot\|$ is the Euclidean norm, $d_{i}\in{\mathbb{R}}^{n}$ and
$w_{i}>0$ for $i=1,\cdots,m$. It can be rewritten as a special case of (M-M)
similar to (1.1).

  * •

For any $D\in{\mathbb{R}}^{m\times n}$, the $\ell_{1}$-norm minimization
problem, with an application in orthogonal dictionary learning [36], can be
reformulated as a case of (M-M)

$\min_{X\in
St_{n,p}}\left\\{\|DX\|_{1}=\sum_{i=1}^{m}\sum_{j=1}^{p}|(DX)_{ij}|=\max_{Y\in\\{\pm
1\\}^{m\times p}}\sum_{i=1}^{m}\sum_{j=1}^{p}Y_{ij}(DX)_{ij}\right\\}.$

  * •

The linear bottleneck assignment problem [11]

$\min_{X\in\Pi_{n}}\max_{i,j=1,\cdots,n}C_{ij}X_{ij}$

is another application of (M-M) based on the equivalence between (2.2) and
(2.3).

The NP-hardness of the binary linear programming in Example 2.1 already
implies that of the general (ELS). We present in the following a complete
proof.

######  Theorem 2.1.

Problem (ELS) with $k\geq n+1$ is NP-hard.

###### Proof.

Given any integer vector $a\in{\mathbb{R}}^{n}$, the partitioning problem asks
whether

$a^{T}y=0,\ y\in{\\{-1,1\\}}^{n}$ (PP)

has a solution. According to (2.1), (PP) is feasible if and only if problem

$\min_{X\in St_{n,1}}\left\\{a^{T}X:~{}a^{T}X\geq 0,\ -\frac{1}{\sqrt{n}}\leq
e_{i}^{T}X\leq\frac{1}{\sqrt{n}},~{}i=1,\cdots,n\right\\}$ (2.4)

has an optimal value $0$. Problem (2.4) can be expanded to the following (ELS)
with $k=n+1$ and any $p\geq 1$:

$\min_{X\in St_{n,p}}\left\\{{\rm tr}(A_{0}X):~{}{\rm tr}(A_{0}X)\geq 0,\
-\frac{1}{\sqrt{n}}\leq{\rm
tr}(E_{i}X)\leq\frac{1}{\sqrt{n}},~{}i=1,\cdots,n\right\\},$

where $A_{0}^{T}=[a~{}0~{}\cdots 0]\in{\mathbb{R}}^{n\times p}$ and
$E_{i}^{T}=[e_{i}~{}0~{}\cdots 0]\in{\mathbb{R}}^{n\times p}$. We complete the
proof by the NP-completeness of the partitioning problem. ∎

##  3  SDP relaxation

In this section, we first present a direct SDP relaxation for (ELS). It can be
reobtained based on the quadratic matrix programming (QMP) reformulation of
(ELS). The exactness of the SDP relaxation follows from relevant results of
QMP.

The nonconvexity of (ELS) comes from the equality quadratic constraints due to
$St_{n,p}$. We start from the convex hull of $St_{n,p}$.

######  Lemma 3.1.

([35, Lemma 1]) The convex hull of $St_{n,p}$ is

$\\{X\in{\mathbb{R}}^{n\times p}:~{}X^{T}X\preceq I_{p}\\}.$ (3.1)

Replacing $St_{n,p}$ with (3.1) yields the conic convex optimization
relaxation:

$\displaystyle\min\limits_{X\in{\mathbb{R}}^{n\times p}}$ $\displaystyle{\rm
tr}(A_{0}X)$ (CR) $\displaystyle\ {\rm s.t.}$ $\displaystyle a_{i}^{l}\leq{\rm
tr}(A_{i}X)\leq a_{i}^{u},~{}i=1,\cdots,k,$ $\displaystyle X^{T}X\preceq
I_{p}.$

For $X\in{\mathbb{R}}^{n\times p}$, define

$Y=\begin{bmatrix}I_{n}&X\\\ X^{T}&I_{p}\end{bmatrix}\in{\mathbb{S}}^{n+p}.$
(3.2)

We have the following congruent transformation:

$Y=\begin{bmatrix}I_{n}&0\\\ X^{T}&I_{p}\end{bmatrix}\begin{bmatrix}I_{n}&0\\\
0&I_{p}-X^{T}X\end{bmatrix}\begin{bmatrix}I_{n}&X\\\ 0&I_{p}\end{bmatrix}.$
(3.3)

Then it holds that

$X^{T}X\preceq I_{p}\Longleftrightarrow Y\succeq 0.$ (3.4)

By (3.2) and (3.4), the conic convex optimization (CR) is equivalent to the
following SDP relaxation:

$\displaystyle\min\limits_{Y\in{\mathbb{S}}^{n+p}}$ $\displaystyle{\rm
tr}(B_{0}Y)$ (SDR) $\displaystyle\ {\rm s.t.}$ $\displaystyle
a_{i}^{l}\leq{\rm tr}(B_{i}Y)\leq a_{i}^{u},~{}i=1,\cdots,k,$ $\displaystyle
Y_{[1:n]}=I_{n},\ Y_{[n+1:n+p]}=I_{p},$ $\displaystyle Y\succeq 0,$

where

$B_{i}=\frac{1}{2}\begin{bmatrix}0_{n}&A_{i}^{T}\\\
A_{i}&0_{p}\end{bmatrix}{\rm\ for\ }i=0,1,\cdots,k.$ (3.5)

We have the following relations between (ELS), (CR), and (SDR).

######  Lemma 3.2.

Suppose $\mathcal{F}\eqref{eq:LPS-SDP3}\neq\emptyset$. Then (SDR) has an
optimal solution and $v\eqref{eq:LPS-R}=v\eqref{eq:LPS-SDP3}$. Moreover, $X$
solves (CR) if and only if $Y$ defined in (3.2) solves (SDR).

###### Proof.

By (3.2), (3.4), and (3.5), $\mathcal{F}\eqref{eq:LPS-SDP3}\neq\emptyset$ if
and only if $\mathcal{F}\eqref{eq:LPS-R}\neq\emptyset$. Let $X$ be any
feasible solution of (CR). Define $Y$ as in (3.2). Then $Y$ is a feasible
solution of (SDR) with the same objective function value as $X$ for (CR). The
converse is also true. Consequently, $v\eqref{eq:LPS-R}=v\eqref{eq:LPS-SDP3}$.
On the other hand, according to $X^{T}X\preceq I_{p}$, we have

$X_{i}^{T}X_{i}\leq 1,~{}i=1,\cdots,p,$

where $X_{i}\in{\mathbb{R}}^{n}$ is the $i$-th column vector of $X$. That is,
$\mathcal{F}\eqref{eq:LPS-R}$ is compact. Hence, $\eqref{eq:LPS-R}$ has an
optimal solution and so does (SDR). ∎

######  Proposition 3.1.

If (SDR) has an optimal solution of rank at most $n$, then (SDR) is tight,
i.e., $v\eqref{eq:LPS-SDP3}=v\eqref{eq:LP_s}$.

###### Proof.

Let $Y$ be any feasible solution of (SDR). Rewrite $Y$ as in (3.2). According
to (3.3), we have

${\rm rank}(Y)={\rm rank}(I_{n})+{\rm rank}(I_{p}-X^{T}X)=n+{\rm
rank}(I_{p}-X^{T}X),$ (3.6)

which implies that

$\displaystyle{\rm rank}(Y)\geq n,$ (3.7) $\displaystyle{\rm
rank}(Y)=n\Longleftrightarrow X^{T}X=I_{p}.$ (3.8)

Let $Y^{*}$ be an optimal solution of (SDR) and ${\rm rank}(Y^{*})\leq n$. It
follows from (3.7) that ${\rm rank}(Y^{*})=n$. By Lemma 3.2,
$X^{*}=Y^{*}_{[1:n][n+1:n+p]}$ is an optimal solution of (CR). Since
$v$(CR)$\leq v$(ELS) and $X^{*}$ remains a feasible solution to (ELS) by
(3.8), we conclude that $v\eqref{eq:LPS-R}=v\eqref{eq:LP_s}$ and then complete
the proof by Lemma 3.2. ∎

The second approach to establish the same SDP relaxation (SDR) is to first
reformulate (ELS) as a QMP problem [4] and then employ its SDP relaxation
proposed in [6]. QMP is the quadratically constrained quadratic optimization
in terms of a matrix variable, say $X\in{\mathbb{R}}^{n\times p}$, where each
involved quadratic function is of the form

${\rm tr}(XQX^{T})+2{\rm tr}(AX)+a,$

where $Q\in{\mathbb{R}}^{p\times p}$, $A\in{\mathbb{R}}^{p\times n}$, and
$a\in{\mathbb{R}}$. The QMP reformulation of (ELS) reads as

$\begin{array}[]{cl}\min\limits_{X\in{\mathbb{R}}^{n\times p}}&{\rm
tr}(A_{0}X)\\\ {\rm s.t.}&a_{i}^{l}\leq{\rm tr}(A_{i}X)\leq
a_{i}^{u},~{}i=1,\cdots,k,\\\ &{\rm tr}(XE^{ij}X^{T})=\delta_{ij},~{}1\leq
i\leq j\leq p,\end{array}$ (3.9)

where $\delta_{ij}=e_{i}^{T}e_{j}$ is the Kronecker delta. The primal lifted
SDP relaxation for the QMP problem (3.9) given in [6] is

$\begin{array}[]{cl}\min\limits_{X\in{\mathbb{R}}^{n\times
p},~{}Z\in{\mathbb{S}}^{p}}&{\rm tr}(A_{0}X)\\\ {\rm s.t.}&a_{i}^{l}\leq{\rm
tr}(A_{i}X)\leq a_{i}^{u},~{}i=1,\cdots,k,\\\ &{\rm
tr}(E^{ij}Z)=\delta_{ij},~{}1\leq i\leq j\leq p,\\\ &\begin{bmatrix}I_{n}&X\\\
X^{T}&Z\end{bmatrix}\succeq 0,\end{array}$ (3.10)

which coincides with (SDR) since the linear constraints ${\rm
tr}(E^{ij}Z)=\delta_{ij}~{}(1\leq i\leq j\leq p)$ imply $Z=I_{p}$. The
exactness of the SDP relaxation for QMP is characterized in [6, Theorem 2.2],
where the following fundamental Shapiro-Barvinok-Pataki rank reduction theorem
[3, 22, 26] plays a key role.

######  Lemma 3.3 ([3, 22, 26]).

Let $n_{1}>n_{2}$. Consider $\mathcal{B}:=\\{Y\in{\mathbb{S}}^{n_{1}}:\
Y\succeq 0,\ a_{i}^{l}\leq{\rm tr}(B_{i}Y)\leq\ a_{i}^{u},~{}i=1,\cdots,m\\}$
and assume $\mathcal{B}\neq\emptyset$. If $m\leq(n_{2}+2)(n_{2}+1)/2-1$, then
${\rm rank}(Y^{e})\leq n_{2}$ for any $Y^{e}\in\mathcal{E}(\mathcal{B})$.

######  Remark 3.1.

The original Shapiro-Barvinok-Pataki theorem works for $\mathcal{B}$ with
$a_{i}^{l}=-\infty,~{}i=1,\cdots,m$. However, its proof allows the slight
extension to any $a_{i}^{l}\leq a_{i}^{u}$.

As a corollary of [6, Theorem 2.2], we immediately have the following
exactness result for (SDR). We present a proof based on Shapiro-Barvinok-
Pataki theorem (Lemma 3.3) for completeness.

######  Proposition 3.2.

Suppose $\mathcal{F}\eqref{eq:LPS-SDP3}\neq\emptyset$. Then $v\eqref{eq:LPS-
SDP3}=v\eqref{eq:LP_s}$ if

$p(p+1)/2\leq n-k.$ (3.11)

###### Proof.

The number of linear constraints of $\eqref{eq:LPS-SDP3}$ is

$m=n(n+1)/2+p(p+1)/2+k.$

It follows from (3.11) that

$m\leq(n+1)(n+2)/2-1.$

According to Lemma 3.3, any $Y^{e}\in\mathcal{E}(\mathcal{F}\eqref{eq:LPS-
SDP3})$ satisfies ${\rm rank}(Y^{e})\leq n$. Since the objective function is
linear, (SDR) has an optimal solution in
$\mathcal{E}(\mathcal{F}\eqref{eq:LPS-SDP3})$. That is, one optimal solution
of (SDR) is of rank less than or equal to $n$. By Proposition 3.1, we complete
the proof. ∎

##  4 Exactness of the SDP relaxation

In this section, we greatly strengthen the condition for the exactness of
(SDR) presented in Proposition 3.2. Then, we use two examples to show that our
new condition is tight.

First, we use the following example to show that Proposition 3.2 already
provides a tight condition for the exactness of (SDR) when $p=1$.

######  Example 4.1.

Consider an instance of (ELS) with $p=1$, $n=2$, and $k=2$:

$\min_{x\in{\mathbb{R}}^{2}}\\{-x_{1}-x_{2}:~{}x_{1}\leq 0,\ x_{2}\leq
0,~{}x^{T}x=1\\}.$ (4.1)

The relaxation (CR) is not exact as one can verify that

$v\eqref{eq:LP_s0}=1>0=\min_{x\in{\mathbb{R}}^{2}}\\{-x_{1}-x_{2}:~{}x_{1}\leq
0,\ x_{2}\leq 0,~{}x^{T}x\leq 1\\}.$

By Lemma 3.2, (SDR) is not exact.

On the other hand, according to Lemmas 3.1 and 3.2, (SDR) is exact for (ELS)
with $n-p=0=k$, which does not satisfy condition (3.11) in Proposition 3.2.
Moreover, when $p\geq\sqrt{2n}$, the condition (3.11) fails to hold even when
$k=0$. One can reasonably expect an improvement of Proposition 3.2.

######  Theorem 4.1.

When $p\leq n-k$, $v\eqref{eq:LPS-SDP3}=v\eqref{eq:LP_s}$.

###### Proof.

Without loss of generalization, we assume that $\mathcal{F}\eqref{eq:LPS-
SDP3}\neq\emptyset$. Otherwise, $v\eqref{eq:LPS-
SDP3}=v\eqref{eq:LP_s}=+\infty$. Since the objective function of (SDR) is
linear, there exists an optimal solution $Y$ of (SDR) satisfies
$Y\in\mathcal{E}(\mathcal{F}\eqref{eq:LPS-SDP3})$. According to Proposition
3.1, it is sufficient to show that ${\rm rank}(Y)\leq n$. Suppose, on the
contrary, it holds that ${\rm rank}(Y)=n+s$ for some $s\in\\{1,\cdots,p\\}$.

Since $Y$ is feasible for (SDR), we can rewrite $Y$ as in (3.2) for some
$X\in{\mathbb{R}}^{n\times p}$. According to (3.4) and $Y\succeq 0$, we have

$I_{p}-X^{T}X\succeq 0.$ (4.2)

It follows from (3.6) and the assumption ${\rm rank}(Y)=n+s$ that

$\ {\rm rank}(I_{p}-X^{T}X)=s.$ (4.3)

Based on (4.2)-(4.3), there exists a full-column-rank matrix
$C\in{\mathbb{R}}^{p\times s}$ such that $I_{p}-X^{T}X=CC^{T}$. Define

$U=\begin{bmatrix}I_{n}&0_{n\times s}\\\
X^{T}&C\end{bmatrix}\in{\mathbb{R}}^{(n+p)\times(n+s)}.$ (4.4)

Then, $U$ is of full column rank, and

$Y=UU^{T}.$ (4.5)

Consider the following linear system in terms of matrix $D$:

$\displaystyle{\rm tr}(B_{i}UDU^{T})=0,~{}i=1,\cdots,k,$ (4.6)
$\displaystyle{UDU^{T}}_{[1:n]}=0_{n},$ (4.7)
$\displaystyle{UDU^{T}}_{[n+1:n+p]}=0_{p},$ (4.8) $\displaystyle
D\in{\mathbb{S}}^{n+s}.$ (4.9)

If (4.6)-(4.9) has a nonzero solution $D$, define

$Y(\epsilon):=Y+\epsilon UDU^{T}=U(I+\epsilon
D)U^{T},~{}\epsilon\in{\mathbb{R}}.$ (4.10)

Since $Y$ is feasible for (SDR), for any $\epsilon$, we have

$\displaystyle{\rm tr}(B_{i}Y(\epsilon))={\rm tr}(B_{i}U(I+\epsilon
D)U^{T})={\rm tr}(B_{i}Y)\in[a_{i}^{l},a_{i}^{u}],~{}i=1,\cdots,k,$ (4.11)
$\displaystyle Y(\epsilon)_{[1:n]}={U(I+\epsilon
D)U^{T}}_{[1:n]}={UU^{T}}_{[1:n]}=I_{n},$ (4.12) $\displaystyle
Y(\epsilon)_{[n+1:n+p]}={U(I+\epsilon
D)U^{T}}_{[n+1:n+p]}={UU^{T}}_{[n+1:n+p]}=I_{p},$ (4.13) $\displaystyle
Y(\epsilon)=U(I+\epsilon D)U^{T}=Y(\epsilon)^{T}.$ (4.14)

Moreover, by the definition of $Y(\epsilon)$, for any $\epsilon$ such that
$|\epsilon|$ is sufficiently small, it holds that

$Y(\epsilon)=U(I+\epsilon D)U^{T}\succeq 0.$ (4.15)

Since $U$ has the full-column rank, it follows from $D\neq 0$ that
$UDU^{T}\neq 0$. By (4.10)-(4.15), $Y(\epsilon)\neq Y$ is feasible for (SDR)
when $\epsilon\neq 0$ and $|\epsilon|$ is sufficiently small. Therefore, we
obtain a contradiction that $Y\not\in\mathcal{E}(\mathcal{F}\eqref{eq:LPS-
SDP3})$, since $Y=(Y(\epsilon)+Y(-\epsilon))/2$.

Now, it suffices to show that the linear system (4.6)-(4.9) has a nonzero
solution. First, we rewrite (4.7) as

${\rm tr}(E^{ij}UDU^{T})=0,\ \forall\ 1\leq i,\ j\leq n.$ (4.16)

By the definition of $U$ (4.4), one can verify that

$U^{T}E^{ij}U=E^{ij},~{}1\leq i,\ j\leq n,$

and hence

$D_{ij}={\rm tr}(E^{ij}D)={\rm tr}(U^{T}E^{ij}UD)={\rm tr}(E^{ij}UDU^{T}).$

Thus, equations (4.16) are equivalent to

$D_{ij}=0,\ \forall\ 1\leq i,j\leq n.$ (4.17)

Second, we rewrite (4.8) as

${\rm tr}(U^{T}E^{(n+i)(n+j)}UD)={\rm tr}(E^{(n+i)(n+j)}UDU^{T})=0,\ \forall\
1\leq i,j\leq p.$ (4.18)

Let $x_{i}\in{\mathbb{R}}^{n}$ be the $i$-th column of $X$ and
$c_{i}\in{\mathbb{R}}^{s}$ be the $i$-th column of $C^{T}$ for $i=1,\cdots,p$.
Then one can verify that

$U^{T}E^{(n+i)(n+j)}U=\begin{bmatrix}x_{i}x_{j}^{T}&x_{i}c_{j}^{T}\\\
c_{i}x_{j}^{T}&c_{i}c_{j}^{T}\end{bmatrix},~{}1\leq i\leq j\leq p.$ (4.19)

Substituting (4.19) into (4.18) yields

${\rm tr}\left(\begin{bmatrix}x_{i}x_{j}^{T}&x_{i}c_{j}^{T}\\\
c_{i}x_{j}^{T}&c_{i}c_{j}^{T}\end{bmatrix}D\right)=0,\ \forall 1\leq i,j\leq
p,$

which are equivalent to

$\displaystyle 0=$ $\displaystyle\ 2{\rm
tr}\left(\begin{bmatrix}x_{i}x_{j}^{T}&x_{i}c_{j}^{T}\\\
c_{i}x_{j}^{T}&c_{i}c_{j}^{T}\end{bmatrix}D\right)$ $\displaystyle=$
$\displaystyle\ 2{\rm tr}\left(\begin{bmatrix}0_{n}&x_{i}c_{j}^{T}\\\
c_{i}x_{j}^{T}&c_{i}c_{j}^{T}\end{bmatrix}D\right)~{}~{}~{}~{}({\rm
since}~{}\eqref{eq:D222})$ $\displaystyle=$ $\displaystyle\ {\rm
tr}\left(\begin{bmatrix}0_{n}&x_{i}c_{j}^{T}\\\
c_{i}x_{j}^{T}&c_{i}c_{j}^{T}\end{bmatrix}D\right)+{\rm
tr}\left(D^{T}\begin{bmatrix}0_{n}&x_{i}c_{j}^{T}\\\
c_{i}x_{j}^{T}&c_{i}c_{j}^{T}\end{bmatrix}^{T}\right)$ $\displaystyle=$
$\displaystyle\ {\rm tr}\left(\begin{bmatrix}0_{n}&x_{i}c_{j}^{T}\\\
c_{i}x_{j}^{T}&c_{i}c_{j}^{T}\end{bmatrix}D\right)+{\rm
tr}\left(\begin{bmatrix}0_{n}&x_{j}c_{i}^{T}\\\
c_{j}x_{i}^{T}&c_{j}c_{i}^{T}\end{bmatrix}D\right)~{}~{}~{}~{}({\rm
since}~{}\eqref{eq:D4})$ $\displaystyle=$ $\displaystyle\ {\rm
tr}\left(\begin{bmatrix}0_{n}&x_{i}c_{j}^{T}\\\
c_{j}x_{i}^{T}&c_{i}c_{j}^{T}\end{bmatrix}D\right)+{\rm
tr}\left(\begin{bmatrix}0_{n}&x_{j}c_{i}^{T}\\\
c_{i}x_{j}^{T}&c_{j}c_{i}^{T}\end{bmatrix}D\right),~{}\forall 1\leq i,j\leq
p.$ (4.20)

Rewrite $D\in{\mathbb{S}}^{n+s}$ as

$\displaystyle D=\begin{bmatrix}D_{11}&D_{12}\\\ D_{21}&D_{22}\end{bmatrix},$
(4.21)

where $D_{11}\in{\mathbb{S}}^{n}$, $D_{12}\in{\mathbb{R}}^{n\times s}$,
$D_{21}=D_{12}^{T}$, $D_{22}\in{\mathbb{S}}^{s}$. Let ${\rm vt}(\cdot)=({\rm
vec}(\cdot))^{T}$. Define

$\displaystyle d$ $\displaystyle=$ $\displaystyle\begin{bmatrix}{\rm
vt}(D_{11})&{\rm vt}(D_{12})&{\rm vt}(D_{21})&{\rm
vt}(D_{22})\end{bmatrix}^{T}\in{\mathbb{R}}^{(n+s)^{2}},$ $\displaystyle
M_{1}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}\cdots&\cdots&\cdots&\cdots\\\ {\rm
vt}(0_{n})&{\rm vt}(x_{i}c_{j}^{T})&{\rm vt}(c_{j}x_{i}^{T})&{\rm
vt}(c_{i}c_{j}^{T})\\\ \cdots&\cdots&\cdots&\cdots\end{bmatrix}_{1\leq i,j\leq
p}\in{\mathbb{R}}^{p^{2}\times(n+s)^{2}},$ (4.22) $\displaystyle M_{2}$
$\displaystyle=$ $\displaystyle\begin{bmatrix}\cdots&\cdots&\cdots&\cdots\\\
{\rm vt}(0_{n})&{\rm vt}(x_{j}c_{i}^{T})&{\rm vt}(c_{i}x_{j}^{T})&{\rm
vt}(c_{j}c_{i}^{T})\\\ \cdots&\cdots&\cdots&\cdots\end{bmatrix}_{1\leq i,j\leq
p}\in{\mathbb{R}}^{p^{2}\times(n+s)^{2}}.$ (4.23)

Since ${\rm tr}(AB^{T})={\rm vt}(A)({\rm vt}(B))^{T}$ for matrices $A$ and
$B$, we can reformulate (4.20) as the following equation

$(M_{1}+M_{2})d=0.$ (4.24)

As $M_{1}$ can be obtained by rearranging the order of the rows of $M_{2}$, we
have

${\rm rank}(M_{1})={\rm rank}(M_{2})={\rm rank}(M_{1}+M_{2}).$ (4.25)

For any $1\leq i\leq p$, define the $i$-th submatrix of $M_{1}$ as

$\displaystyle M_{1}^{i}=$
$\displaystyle\begin{bmatrix}\cdots&\cdots&\cdots&\cdots\\\ {\rm
vec}(0_{n})&{\rm vec}(x_{i}c_{j}^{T})&{\rm vec}(c_{j}x_{i}^{T})&{\rm
vec}(c_{i}c_{j}^{T})\\\ \cdots&\cdots&\cdots&\cdots\end{bmatrix}_{1\leq j\leq
p}\in{\mathbb{R}}^{p\times(n+s)^{2}}.$

Since rank$([c_{1}~{}c_{2}~{}\cdots~{}c_{p}]^{T})=$rank$(C)=s$, we have

${\rm rank}(M_{1}^{i})\leq s,~{}1\leq i\leq p.$ (4.26)

Then, it follows from (4.25)-(4.26) that

${\rm rank}(M_{1}+M_{2})={\rm rank}(M_{1})\leq\sum_{i=1}^{p}{\rm
rank}(M_{1}^{i})\leq ps.$ (4.27)

In sum, the linear equations (4.7)-(4.8) are equivalent to (4.17) and (4.24),
respectively. The rank of the coefficient matrix of (4.17) is $n^{2}$, while
the coefficient matrix of (4.24) has a rank no more than $ps$. The symmetry
condition (4.9) is equivalent to a system with $(n+s)(n+s-1)/2$ linear
equalities, where $n(n-1)/2$ of them are already contained in (4.17).

To sum up the above analysis, (4.6)-(4.9) amounts to a linear system in terms
of $D\in{\mathbb{R}}^{(n+s)^{2}}$ with a coefficient matrix whose rank $R$
satisfies that

$\displaystyle R\leq\ $ $\displaystyle
k+n^{2}+ps+\frac{(n+s)(n+s-1)}{2}-\frac{n(n-1)}{2}$ $\displaystyle=$
$\displaystyle\ k+\frac{n(n+1)}{2}+ps+\frac{(n+s)(n+s-1)}{2}.$

Under the assumption $p\leq n-k$, we have

$\displaystyle(n+s)^{2}-R=$
$\displaystyle\frac{(n+s)(n+s+1)}{2}-\left[k+\frac{n(n+1)}{2}+ps\right]$
(4.28) $\displaystyle\geq$ $\displaystyle\frac{s(s+1)}{2}+(n-p)(s-1)$
$\displaystyle\geq$ $\displaystyle 1,~{}~{}~{}~{}~{}({\rm since}\ s\geq
1,~{}p\leq n).$

Consequently, (4.6)-(4.9) must have a nonzero solution. The proof is complete.
∎

######  Remark 4.1.

According to the proofs of Proposition 3.1 and Theorem 4.1, when $p\leq n-k$,
(SDR) has an optimal solution of rank $n$, and (CR) has an optimal solution in
$E_{n,p,k}$.

######  Remark 4.2.

Suppose $p\leq n-k$. The proof of Theorem 4.1 indicates a rank reduction
algorithm to obtain an optimal rank-$n$ solution of (SDR) similarly as in [4].
Let $Y$ be any optimal solution to (SDR). If rank$(Y)>n$, we first solve
(4.7)-(4.9) to get a nonzero $D$ and define $Y(\epsilon)$ as in (4.10). Since
$Y$ is optimal and $Y(\epsilon)$ is feasible for sufficiently small
$|\epsilon|$, we have

${\rm tr}(B_{0}Y)={\rm tr}(B_{0}Y(\epsilon)),~{}\forall\epsilon.$ (4.29)

Then we can find a suitable $\widetilde{\epsilon}$ satisfying (4.15) and
rank$(Y(\widetilde{\epsilon}))\leq$rank$(Y)-1$. Set
$Y:=Y(\widetilde{\epsilon})$, it remains optimal by (4.29). This procedure is
repeated until rank$(Y)=n$.

When $p=1$, Theorem 4.1 coincides with Proposition 3.2, and the condition is
tight as shown in Example 4.1. When $p>1$, Theorem 4.1 strictly improves
Proposition 3.2. We illustrate the tightness of condition $p\leq n-k$
presented in Theorem 4.1 by two more examples.

######  Example 4.2.

Consider an instance of (ELS) with $p=2$, $n=3$, and $k=2$:

$\min_{X\in{\mathbb{R}}^{3\times
2}}\\{X_{3,2}:~{}X_{1,1}=X_{2,1}=0,~{}X^{T}X=I_{2}\\}.$ (4.30)

One can easily verify that $v\eqref{eq:LP_s1}=0$ as

$\mathcal{F}\eqref{eq:LP_s1}=\\{X\in{\mathbb{R}}^{3\times
2}:X_{1,1}=X_{2,1}=X_{3,2}=0,\ X_{3,1}=\pm 1,\ X_{1,2}^{2}+X_{2,2}^{2}=1\\}.$

The relaxation (CR) is given by

$\min_{X\in{\mathbb{R}}^{3\times
2}}\\{X_{3,2}:~{}X_{1,1}=X_{2,1}=0,~{}X^{T}X\preceq I_{2}\\},$ (4.31)

where

$\bar{X}=\begin{bmatrix}0&0\\\ 0&0\\\ 0&-1\end{bmatrix}.$

is feasible. It follows that

$v\eqref{eq:LP_s1d}\leq\bar{X}_{3,2}=-1<0=v\eqref{eq:LP_s1}.$

That is, (4.31) is not an exact relaxation for (4.30), so is the (SDR)
relaxation due to Lemma 3.2.

######  Example 4.3.

Consider an instance of (ELS) with $p=3$, $n=3$, and $k=1$:

$\min_{X\in{\mathbb{R}}^{3\times
3}}\\{X_{2,2}+X_{3,3}:~{}X_{1,1}=0,~{}X^{T}X=I_{3}\\}.$ (4.32)

The (CR) relaxation

$\min_{X\in{\mathbb{R}}^{3\times
3}}\\{X_{2,2}+X_{3,3}:~{}X_{1,1}=0,~{}X^{T}X\preceq I_{3}\\}$ (4.33)

has a feasible solution

$\bar{X}=\begin{bmatrix}0&0&0\\\ 0&-1&0\\\ 0&0&-1\end{bmatrix}.$ (4.34)

Then it holds that $v\eqref{eq:LP_s2d}\leq-2$. On the other hand, for any
$X^{T}X\preceq I_{3}$, we have $X_{2,2}^{2}\leq 1$ and $X_{3,3}^{2}\leq 1$,
which imply that $v\eqref{eq:LP_s2d}\geq-2$. Therefore,
$v\eqref{eq:LP_s2d}=-2$ and $\bar{X}$ is an optimal solution of (4.33).

It is not difficult to verify that

$X=\begin{bmatrix}0&X_{1,2}&X_{1,3}\\\ X_{2,1}&-1&X_{2,3}\\\
X_{3,1}&X_{3,2}&-1\end{bmatrix}\in{\mathbb{R}}^{3\times 3}$ (4.35)

cannot be a feasible solution to (4.32). It follows that
$v\eqref{eq:LP_s2}>-2=v\eqref{eq:LP_s2d}$. That is, (4.33) is not an exact
relaxation for (4.32), so is the (SDR) relaxation due to Lemma 3.2.

We can apply Theorem 4.1 to reveal the hidden convexity of Example 2.6 under
certain conditions.

######  Corollary 4.1.

Suppose $\mathcal{F}\eqref{eq:LPS-SDP3}\neq\emptyset$. When $p\leq n-k-m+1$,
the problem (M-M) admits the following tight SDP relaxation

$\displaystyle\min\limits_{Y\in{\mathbb{S}}^{n+p},\,t\in{\mathbb{R}}}$
$\displaystyle t$ (4.36) $\displaystyle\ {\rm s.t.}$ $\displaystyle{\rm
tr}(A_{i}Y_{[1:n][n+1:n+p]})+c_{i}\leq t,~{}i=1,\cdots,m,$ $\displaystyle
Y\in\mathcal{F}\eqref{eq:LPS-SDP3}.$

###### Proof.

Clearly, (4.36) is an SDP relaxation for (M-M). Suppose
$\mathcal{F}\eqref{eq:LPS-SDP3}\neq\emptyset$, we have
$\mathcal{F}\eqref{eq:mm:h1}\neq\emptyset$ by selecting $t$ large enough.
Based on a proof similar to that of the last part of Lemma 3.2, we can show
that (4.36) has an optimal solution. According to (3.8), it suffices to show
that there exists an optimal solution of rank $n$. Let
$(Y^{*},t^{*})\in{\mathbb{S}}^{n+p}\times{\mathbb{R}}$ be a global minimizer
of (4.36). There is at least one index $q\in\\{1,\cdots,m\\}$ such that

${\rm tr}(A_{q}Y^{*}_{[1:n][n+1:n+p]})+c_{q}=t^{*}\geq{\rm
tr}(A_{i}Y^{*}_{[1:n][n+1:n+p]})+c_{i},~{}\forall\,i\in\\{1,\cdots,m\\}\setminus\\{q\\}.$

Then, $Y^{*}$ is also an optimal solution of the following SDP problem

$\displaystyle\min$ $\displaystyle{\rm tr}(A_{q}Y_{[1:n][n+1:n+p]})$ (4.37)
$\displaystyle\ {\rm s.t.}$ $\displaystyle{\rm
tr}(A_{q}Y_{[1:n][n+1:n+p]})+c_{q}\geq{\rm
tr}(A_{j}Y_{[1:n][n+1:n+p]})+c_{j},~{}\forall j\neq q,$ $\displaystyle
Y\in\mathcal{F}\eqref{eq:LPS-SDP3},$

which contains $k+m-1$ linear constraints. According to Theorem 4.1 and Remark
4.1, under the assumption $p\leq n-(k+m-1)$, (4.37) has an optimal solution of
rank $n$. Since $v\eqref{eq:mm:h1}=v$(4.37), the proof is complete. ∎

As a byproduct of Theorem 4.1, we characterize the convexity result of a joint
numerical range.

######  Theorem 4.2.

If $p\leq n-k$, the set

$G_{1}=\\{[{\rm tr}(A_{1}X),\ \cdots,\ {\rm
tr}(A_{k}X)]:~{}X^{T}X=I_{p},~{}X\in{\mathbb{R}}^{n\times p}\\}$ (4.38)

is convex, and is equal to

$G_{2}=\\{[{\rm tr}(A_{1}X),\ \cdots,\ {\rm tr}(A_{k}X)]:~{}X^{T}X\preceq
I_{p},~{}X\in{\mathbb{R}}^{n\times p}\\}.$ (4.39)

###### Proof.

It is easy to see $G_{1}\subseteq G_{2}$ and $G_{2}$ is convex. It suffices to
show $G_{2}\subseteq G_{1}$. For any $a\in G_{2}$, we have

$\mathcal{X}_{a}:=\\{X\in{\mathbb{R}}^{n\times p}:~{}X^{T}X\preceq I_{p},\
{\rm tr}(A_{i}X)=a_{i},~{}i=1,\cdots,k\\}\neq\emptyset.$ (4.40)

Then, it holds that

$\min\\{{\rm tr}(0\cdot X):~{}X\in\mathcal{X}_{a}\\}=0.$

Consider this special case of (CR). According to Theorem 4.1, Remark 4.1, and
the assumption $p\leq n-k$, there exists a matrix
$\bar{X}\in{\mathbb{R}}^{n\times p}$ such that

$\bar{X}^{T}\bar{X}=I_{p},\ {\rm tr}(A_{i}\bar{X})=a_{i},~{}i=1,\cdots,k.$

That is, $a\in G_{1}$. We complete the proof. ∎

######  Remark 4.3.

Beck [5] established convexity of the following joint numerical range

$\\{[x^{T}x+a_{0}^{T}x,\ a_{1}^{T}x,\ \cdots,\
a_{k}^{T}x]:~{}x\in{\mathbb{R}}^{n}\\}\subseteq{\mathbb{R}}^{m+1}$ (4.41)

where $a_{0},\cdots,a_{k}\in{\mathbb{R}}^{n}$ and $k\leq n-1$. Up to a
translation transformation, we can set $a_{0}\equiv 0$ in (4.41). The
convexity of (4.41) implies that of $G_{1}$ (4.38) with $p=1$.

######  Remark 4.4.

Define $B_{1},\cdots,B_{k}$ as in (3.5). By (3.2), (3.4) and (3.8), $G_{1}$
(4.38) is convex if and only if

$\displaystyle\\{[{\rm tr}(B_{1}Y)$ $\displaystyle,\ \cdots,\ {\rm
tr}(B_{k}Y)]:~{}Y\in{\mathbb{S}}^{n+p},~{}Y\succeq 0,\ {\rm rank}(Y)=n,$
(4.42) $\displaystyle{\rm tr}(E^{i,j}Y)=\delta_{i,j},~{}\forall 1\leq i\leq
j\leq n,~{}n+1\leq i\leq j\leq n+p\\}$

is convex. According to Theorem 4.2, (4.42) is convex if $p\leq n-k$.
Convexity of (4.42) can also be obtained from that of

$\displaystyle\\{[{\rm tr}(B_{1}Y),\cdots,{\rm tr}(B_{k}Y),$
$\displaystyle~{}{\rm tr}(E_{i,j}Y)\ (1\leq i\leq j\leq n,~{}n+1\leq i\leq
j\leq n+p)]:$ (4.43) $\displaystyle Y\in{\mathbb{S}}^{n+p},~{}Y\succeq
0,~{}{\rm rank}(Y)=n\\},$

which is closely related to the following classical result [2].

######  Lemma 4.1.

[2, Theorem 3] If $1\leq n_{2}\leq n_{1}-1$ and $l\leq(n_{2}+2)(n_{2}+1)/2-1$,
then for $B_{1},\cdots,B_{l}\in{\mathbb{S}}^{n_{1}}$, the set

$\\{[{\rm tr}(B_{1}Y),\cdots,{\rm tr}(B_{l}Y)]:~{}Y\in{\mathbb{S}}^{n_{1}},\
Y\succeq 0,\ {\rm rank}(Y)=n_{2}\\}$

is convex, and is equal to

$\\{[{\rm tr}(B_{1}Y),\cdots,{\rm tr}(B_{l}Y)]:~{}Y\in{\mathbb{S}}^{n_{1}},\
Y\succeq 0\\}.$

In the case (4.43), $n_{1}=n+p,\ n_{2}=n,\ l=k+p(p+1)/2+n(n+1)/2$. Under the
additional assumption (3.11), one can verify that

$\displaystyle l-((n_{2}+1)(n_{2}+2)/2-1)$ $\displaystyle=$ $\displaystyle
k+p(p+1)/2+n(n+1)/2-((n+1)(n+2)/2-1)$ $\displaystyle\leq$ $\displaystyle
n+n(n+1)/2-(n+1)(n+2)/2+1$ $\displaystyle=$ $\displaystyle 0.$

Therefore, (4.42) is convex by Lemma 4.1. Comparing with Theorem 4.2, the
condition (3.11) is too restrictive when $p>1$.

##  5 Sufficient global optimality based on local necessary conditions

In this section, we regard (ELS) as a nonlinear programming problem, and focus
on studying whether the local necessary optimality conditions could guarantee
global optimality.

We have revealed hidden convexity of (ELS) when $p\leq n-k$ in the previous
section. However, as shown by the following example, it does not mean that any
local minimizer of (ELS) is globally optimal.

######  Example 5.1.

Consider an instance of (ELS) with $p=1$, $n=2$ and $k=1$:

$\min_{x\in{\mathbb{R}}^{2}}\\{-x_{1}-2x_{2}:~{}x_{1}\leq 0,~{}x^{T}x=1\\}.$
(5.1)

It is not difficult to verify that $(0,1)^{T}$ is the unique global minimizer,
and $(0,-1)^{T}$ is a local non-global minimizer.

As a well-known result, if the KKT condition holds and the corresponding
Lagrangian function is convex with respect to the primal variables, then the
global optimality holds. For completeness, we present the corresponding result
for (ELS) with a detailed proof.

######  Lemma 5.1.

Let $X^{*}\in E_{n,p,k}$. If there exist $\lambda_{i}\in{\mathbb{R}}$ for
$i=1,\cdots,k$ and $\Lambda\succeq 0$ such that

$\displaystyle A_{0}^{T}+\sum_{i=1}^{k}\lambda_{i}A_{i}^{T}+X^{*}\Lambda=0,$
(5.2) $\displaystyle{\rm tr}(A_{i}X^{*})-a_{i}\leq 0,\ i=1,\cdots,k,\
X^{*T}X^{*}=I_{p},$ (5.3) $\displaystyle\lambda_{i}^{+}({\rm
tr}(A_{i}X^{*})-a_{i}^{u})=\lambda_{i}^{-}({\rm tr}(A_{i}X^{*})-a_{i}^{l})=0,\
i=1,\cdots,k,$ (5.4)

then $X^{*}$ is a global minimizer of (ELS).

###### Proof.

For any $X\in E_{n,p,k}$, we have

$\displaystyle{\rm tr}(A_{0}X)$ $\displaystyle\geq\ $ $\displaystyle{\rm
tr}(A_{0}X)+\sum_{i=1}^{k}\left[\lambda_{i}^{+}({\rm
tr}(A_{i}X)-a_{i}^{u})-\lambda_{i}^{-}({\rm tr}(A_{i}X)-a_{i}^{l})\right]+{\rm
tr}(\Lambda(X^{T}X-I_{p}))$ (5.5) $\displaystyle=\ $ $\displaystyle{\rm
tr}(A_{0}X)+\sum_{i=1}^{k}\lambda_{i}{\rm tr}(A_{i}X)+{\rm
tr}(\Lambda(X^{T}X-I_{p}))-\sum_{i=1}^{k}(\lambda_{i}^{+}a_{i}^{u}-\lambda_{i}^{-}a_{i}^{l})$
(5.6) $\displaystyle\geq\ $ $\displaystyle{\rm
tr}(A_{0}X^{*})+\sum_{i=1}^{k}\lambda_{i}{\rm tr}(A_{i}X^{*})+{\rm
tr}(\Lambda(X^{*T}X^{*}-I_{p}))-\sum_{i=1}^{k}(\lambda_{i}^{+}a_{i}^{u}-\lambda_{i}^{-}a_{i}^{l})$
(5.7) $\displaystyle=\ $ $\displaystyle{\rm
tr}(A_{0}X^{*})+\sum_{i=1}^{k}\left[\lambda_{i}^{+}({\rm
tr}(A_{i}X^{*})-a_{i}^{u})-\lambda_{i}^{-}({\rm
tr}(A_{i}X^{*})-a_{i}^{l})\right]+{\rm tr}(\Lambda(X^{*T}X^{*}-I_{p}))$
$\displaystyle=\ $ $\displaystyle{\rm tr}(A_{0}X^{*}),$ (5.8)

where (5.5) holds due to the feasibility of $X$, the inequality (5.7) holds
since $X^{*}$ globally minimizes the unconstrained convex quadratic
programming (5.6) as $\Lambda\succeq 0$ and (5.2), the last equality (5.8)
holds due to (5.3) and (5.4). Therefore, $X^{*}$ is a global minimizer of
(ELS). ∎

Regrading (ELS) as a nonlinear programming problem in terms of the vector
variables in ${\mathbb{R}}^{np}$, based on the standard first- and second-
order necessary optimality conditions, we easily derive the following result.

######  Lemma 5.2.

Assume $p(p+1)/2+1\leq n-k$. Let $X^{*}\in E_{n,p,k}$ satisfy LICQ, and the
standard first- and second-order necessary optimality conditions. Then $X^{*}$
is a global minimizer of (ELS).

###### Proof.

Under the assumptions, there exists $\lambda_{i}\in{\mathbb{R}}$ for
$i=1,\cdots,k$ and $\Lambda\in{\mathbb{S}}^{p}$ such that the KKT condition
(5.2)-(5.4) holds. Let

$\mathcal{A}=\\{i:{\rm tr}(A_{i}X^{*})-a_{i}^{l}=0\ {\rm or}\ {\rm
tr}(A_{i}X^{*})-a_{i}^{u}=0\\}.$ (5.9)

Let $(\cdot)_{i}$ be the $i$-th column of matrix $(\cdot)$. According to the
second-order necessary optimality condition, for any
$V\in{\mathbb{R}}^{n\times p}$ such that

$\displaystyle{\rm tr}(A_{i}V)=0,\ i\in\mathcal{A},$ (5.10) $\displaystyle
V_{i}^{T}X^{*}_{j}+V_{j}^{T}X^{*}_{i}=0,\ 1\leq i\leq j\leq p,$ (5.11)

it must hold that

$\displaystyle{\rm vec}(V)^{T}(\Lambda\otimes I_{n}){\rm vec}(V)\geq 0.$
(5.12)

Note that $V\in{\mathbb{R}}^{n\times p}$ has $np$ components and (5.10)-(5.11)
contains at most $k+p(p+1)/2$ linear constraints. Under the assumption
$p(p+1)/2+1\leq n-k$, we have

$np-k-p(p+1)/2\geq np-(n-1).$

Therefore, the set of $V$ satisfying (5.10)-(5.11) is a subspace with
dimension no less than $np-(n-1).$ Then, according to (5.12), $\Lambda\otimes
I_{n}\in{\mathbb{R}}^{np\times np}$ has at most $n-1$ negative eigenvalues. On
the other hand, if $\Lambda$ has at least one negative eigenvalue, then
$\Lambda\otimes I_{n}$ has at least $n$ negative eigenvalues. Thus, it follows
that $\Lambda\succeq 0$. According to Lemma 5.1, $X^{*}$ is a global minimizer
of (ELS). The proof is complete. ∎

To our surprise, the above result in Lemma 5.2 can be greatly strengthened by
more careful analysis.

######  Theorem 5.1.

Suppose $p+1\leq n-k$. Let $X^{*}\in E_{n,p,k}$ satisfy LICQ, and the first-
and second-order necessary optimality conditions. Then, $X^{*}$ is a global
minimizer of (ELS).

###### Proof.

Let $X^{*}=[X_{1}^{*},\cdots,X_{p}^{*}]$. There are vectors
$X_{p+1}^{*},\cdots,X_{n}^{*}\in{\mathbb{R}}^{n}$ such that
$Q:=[X_{1}^{*},\cdots,$ $X_{p}^{*},X_{p+1}^{*},\cdots,X_{n}^{*}]$ is
orthogonal. Consider the orthogonal linear transformation $Z=Q^{T}X$. We have

$\displaystyle X^{T}X=I_{p}\Longleftrightarrow Z^{T}Z=I_{p},$
$\displaystyle{\rm tr}(A_{i}X)={\rm tr}(A_{i}QZ),\ i=0,\cdots,k.$

Define $\tilde{A}_{i}=A_{i}Q$ for $i=0,\cdots,k$. Then,
$Z^{*}=Q^{T}X^{*}=[e_{1},\cdots,e_{p}]$ is a global minimizer of the following
(ELS) problem:

$\begin{array}[]{cl}\min&{\rm tr}(\tilde{A}_{0}Z)\\\ {\rm s.t.}&Z\in\\{Z\in
St_{n,p}:~{}a_{i}^{l}\leq{\rm tr}(\tilde{A}_{i}Z)\leq
a_{i}^{u},~{}i=1,\cdots,k\\}.\end{array}$

Moreover, under any nonsingular linear transformation, LICQ and the first- and
second-order necessary optimality conditions remain invariant. Without loss of
generality, we simply assume that

$X^{*}=[e_{1},\cdots,e_{p}]$

is a global minimizer of the original (ELS).

Let $V\in{\mathbb{R}}^{n\times p}$ satisfy (5.10)-(5.11). Note that (5.11)
reduces to

$V_{i}^{T}e_{j}+V_{j}^{T}e_{i}=0,\ 1\leq i\leq j\leq p,$ (5.13)

where $V_{i(j)}$ is the $i(j)$-th column of $V$. Define
$U=V_{[1:p][1:p]}\in{\mathbb{R}}^{p\times p}$ and
$W=V_{[p+1:n][1:p]}\in{\mathbb{R}}^{(n-p)\times p}$. That is, $V=[U;W].$ The
equations (5.13) are equivalent to

$U_{ij}+U_{ji}=0,\ 1\leq i\leq j\leq p,$ (5.14)

i.e., $U$ is skew-symmetric.

In sum, by collecting all components $U_{ij}$ $(1\leq i<j\leq n)$ as
$u\in{\mathbb{R}}^{p(p-1)/2}$, we reorganize the second-order optimality
condition (5.10)-(5.12), and then rewrite it as

$[u^{T},\ {\rm vec}(W)^{T}]b_{i}=0,\ i\in\mathcal{A}\Longrightarrow[u^{T},\
{\rm vec}(W)^{T}]M\left[\begin{matrix}u\\\ {\rm
vec}(W)\end{matrix}\right]\succeq 0,$ (5.15)

where $\mathcal{A}$ is defined in (5.9), $b_{i}$ and
$M\in{\mathbb{S}}^{p(p-1)/2+(n-p)p}$ are the rearranged counterparts of
$A_{i}$ and $\Lambda\otimes I_{n}$, respectively.

Since $\mathcal{A}$ contains at most $k$ indices and we have assumed $k\leq
n-p-1$, it follows from (5.15) that $M$ has at most $n-p-1$ negative
eigenvalues. On the other hand, we can verify that

$M_{[p(p-1)/2+1:p(p-1)/2+(n-p)p]}={(\Lambda\otimes
I_{n})}_{[p^{2}+1:p^{2}+(n-p)p]}=\Lambda\otimes I_{n-p}.$ (5.16)

Thus, if $\Lambda$ has a negative eigenvalue, then $\Lambda\otimes I_{n-p}$
has at least $n-p$ negative eigenvalues and so does
$M_{[p(p-1)/2+1:p(p-1)/2+(n-p)p]}$ by (5.16). According to Cauchy’s interlace
theorem [17], $M$ has at least $n-p$ negative eigenvalues. The contradiction
implies that $\Lambda\succeq 0$. Then, according to Lemma 5.1, we complete the
proof. ∎

######  Remark 5.1.

According to Example 5.1, the condition $p+1\leq n-k$ assumed in Theorem 5.1
seems to be tight.

According to Theorem 5.1, we present a sufficient condition under which there
is no local non-global minimizer.

######  Corollary 5.1.

If $p+1\leq n-k$, and LICQ holds at any feasible point of (ELS), then any
local minimizer of (ELS) is globally optimal.

As a generalization of LICQ, the condition that the rank of the Jacobin matrix
of constraints is invariant for all points in $E_{n,p,k}$ is sufficient for
$E_{n,p,k}$ being an embedded submanifold in $St_{n,p}$ [18, Theorem 5.12].

######  Remark 5.2.

As shown in the following example, the rank of Jacobin matrix of constraints
may not be the same for different feasible points of the general (ELS). It
indicates that the intersection of a hyperplane and a Stiefel manifold could
not remain a manifold [1].

######  Example 5.2.

Consider an instance of (ELS) with $p=2$, $n=4$, and $k=1$:

$\min_{X\in{\mathbb{R}}^{4\times 2}}\\{0:~{}{\rm tr}(AX)=1,\ X^{T}X=I_{2}\\},$
(5.17)

where

$A=\begin{bmatrix}1&1&0&0\\\ 1&0&0&0\end{bmatrix}.$ (5.18)

Given the following two feasible points

$\bar{X}=\begin{bmatrix}\bar{X}_{1}&\bar{X}_{2}\end{bmatrix}=\begin{bmatrix}1&0\\\
0&1\\\ 0&0\\\ 0&0\end{bmatrix},\
\tilde{X}=\begin{bmatrix}\tilde{X}_{1}&\tilde{X}_{2}\end{bmatrix}=\begin{bmatrix}\frac{1}{2}&0\\\
\frac{1}{2}&\frac{\sqrt{2}}{2}\\\ 0&\frac{1}{2}\\\
\frac{\sqrt{2}}{2}&-\frac{1}{2}\end{bmatrix},$ (5.19)

the Jacobin matrices of the constraints are

$\displaystyle\bar{J}=\begin{bmatrix}\bar{X}_{1}&0&\bar{X}_{2}&{(A^{T})}_{1}\\\
0&\bar{X}_{2}&\bar{X}_{1}&{(A^{T})}_{2}\end{bmatrix}=\begin{bmatrix}1&0&0&1\\\
0&0&1&1\\\ 0&0&0&0\\\ 0&0&0&0\\\ 0&0&1&1\\\ 0&1&0&0\\\ 0&0&0&0\\\
0&0&0&0\end{bmatrix},$ (5.20) $\displaystyle{\rm and~{}}$
$\displaystyle\tilde{J}=\begin{bmatrix}\tilde{X}_{1}&0&\tilde{X}_{2}&{A^{T}}_{1}\\\
0&\tilde{X}_{2}&\tilde{X}_{1}&{A^{T}}_{2}\end{bmatrix}=\begin{bmatrix}\frac{1}{2}&0&0&1\\\
\frac{1}{2}&0&\frac{\sqrt{2}}{2}&1\\\ 0&0&\frac{1}{2}&0\\\
\frac{\sqrt{2}}{2}&0&-\frac{1}{2}&0\\\ 0&0&\frac{1}{2}&1\\\
0&\frac{\sqrt{2}}{2}&\frac{1}{2}&0\\\ 0&\frac{1}{2}&0&0\\\
0&-\frac{1}{2}&\frac{\sqrt{2}}{2}&0\end{bmatrix},$ (5.21)

respectively. One can easily verify that rank$(\bar{J})=3$ and
rank$(\tilde{J})=4$.

######  Remark 5.3.

There are two trivial cases of (ELS) where the feasible sets are manifolds:
(i) $k=0$; (ii) $p=1$. In case (i), the feasible set is $St_{n,p}$, and so
LICQ holds at any feasible point. In case (ii), the feasible set, if not
empty, is either a singleton or a sphere.

##  6 Conclusion

We have introduced a new optimization model (ELS) by adding $k$ linear
constraints in the classical linear optimization on Stiefel manifold
$St_{n,p}$. We provided some applications of (ELS), from which we showed that
(ELS) is NP-hard in general. Regrading (ELS) as a nonconvex quadratic matrix
programming and following Shapiro-Barvinok-Pataki theorem, one can verify that
(ELS) admits an exact SDP relaxation when $p(p+1)/2\leq n-k$. Our analysis
showed that this sufficient condition can be greatly strengthened to be $p\leq
n-k$, which includes the extreme case $n-p=0=k$. It should be noted that the
hidden convexity does not mean any local minimizer of (ELS) is globally
optimal. Regarding (ELS) as a smooth nonlinear programming problem, we showed
that, at a local minimizer satisfying LICQ, the standard first- and second-
order necessary optimality conditions are sufficient for the global optimality
when $p+1\leq n-k$. When and only when LICQ holds at all feasible points of
(ELS) is unknown. Another theoretical question is when and only when the
intersection of a hyperplane and $St_{n,p}$ remains a manifold.

## References

  * [1] P. Absil, R. Mahony, and R. Sepulchre, Optimization Algorithms on Matrix Manifolds, Princeton University Press, Princeton, 2007.
  * [2] Y.-H. Au-Yeung and Y. T. Poon., A remark on the convexity and positive definiteness concerning Hermitian matrices, Southeast Asian Bulletin of Mathematics, 3 (1979), pp. 85–92.
  * [3] A. Barvinok, Problems of distance geometry and convex properties of quadratic maps, Discrete and Computational Geometry, 13 (1995), pp. 189–202.
  * [4] A. Beck, Quadratic matrix programming, SIAM Journal on Optimization, 17 (2006), pp. 1224–1238.
  * [5] A. Beck, On the convexity of a class of quadratic mappings and its application to the problem of finding the smallest ball enclosing a given intersection of balls, Journal of Global Optimization, 39 (2007), pp. 113–126.
  * [6] A. Beck, Convexity properties associated with nonconvex quadratic matrix functions and applications to quadratic programming, Journal of Optimization Theory and Applications, 142 (2009), pp. 1–29.
  * [7] D. Bienstock and A. Michalka, Polynomial solvability of variants of the trust-region subproblem, in Proceedings of the twenty-fifth annual ACM-SIAM symposium on Discrete algorithms, 2014, pp. 380–390.
  * [8] G. Birkhoff, Tres observaciones sobre el algebre lineal, Univ. Nac. Tucumán Rev. Ser. A, 5 (1946).
  * [9] N. Boumal, V. Voroninski, and A. S. Bandeira, Deterministic guarantees for Burer-Monteiro factorizations of smooth semidefinite programs, Communications on Pure and Applied Mathematics, 73 (2020), pp. 581–608.
  * [10] S. Burer and R. Monteiro, A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization, Mathematical Programming, 95 (2003), pp. 329–357.
  * [11] R. Burkard and U. Derigs, Assignment and Matching Problems: Solution Methods with FORTRAN-Programs, Lecture Notes in Economics and Mathematical Systems, vol 184, Springer, Berlin, Heidelberg, 1980, ch. The Linear Bottleneck Assignment Problem.
  * [12] R. E. Burkard and E. Çela, Linear assignment problems and extensions, in Handbook of Combinatorial Optimization: Supplement Volume A, D.-Z. Du and P. M. Pardalos, eds., Springer US, Boston, MA, 1999, pp. 75–149.
  * [13] S. Chen, S. Ma, A. M. So, and T. Zhang, Proximal gradient method for nonsmooth optimization over the Stiefel manifold, SIAM Journal on Optimization, 30(1) (2020), pp. 210–239.
  * [14] L. Eldén and H. Park, A Procrustes problem on the Stiefel manifold, Numerische Mathematik, 82 (1999), pp. 599–619.
  * [15] G. Golub, Some modified eigenvalue problem, SIAM Review, 15 (1973), pp. 318–344.
  * [16] Y. Hsia and R.-L. Sheu, Trust region subproblem with a fixed number of additional linear inequality constraints has polynomial complexity.  Preprint, arXiv:1312.1398, 2013.
  * [17] S.-G. Hwang, Cauchy’s interlace theorem for eigenvalues of Hermitian matrices, The American Mathematical Monthly, 111 (2004), pp. 157–159.
  * [18] J. M. Lee, Introduction to Smooth Manifolds, Springer New York, NY, 2009\.
  * [19] Y. Li and W. Xie, On the exactness of Dantzig-Wolfe relaxation for rank constrained optimization problems.  arXiv:2210.16191, 2022.
  * [20] J. Manton, Optimization algorithms exploiting unitary constraints, IEEE Trans. Signal Process., 50(3) (2002), pp. 635–650.
  * [21] J. M. Martínez, Local minimizers of quadratic functions on Euclidean balls and spheres, SIAM Journal on Optimization, 4 (1994), pp. 159–176.
  * [22] G. Pataki, On the rank of extreme matrices in semidefinite programs and the multiplicity of optimal eigenvalues, Mathematics of Operations Research, 23 (1998), pp. 339–358.
  * [23] T. Rapcsák, On minimization of sums of heterogeneous quadratic functions on Stiefel manifolds, in From Local to Global Optimization, A. Migdalas, P. Pardalos, and P. Varbrand, eds., Kluwer Academic Publishers, Dordrecht, 2001.
  * [24] T. Rapcsák, On minimization on Stiefel manifolds, European Journal of Operational Research, 143 (2002), pp. 365–376.
  * [25] P. Schonemann, A generalized solution of the orthogonal Procrustes problem, Psychometrika, 31(1) (1966), pp. 1–10.
  * [26] A. Shapiro, Rank-reducibility of a symmetric matrix and sampling theory of minimum trace factor analysis, Psychometrika, 47 (1982), pp. 187–199.
  * [27] M. Song, H. Liu, J. Wang, and Y. Xia, On local minimizers of nonconvex quadratically constrained homogeneous quadratic optimization with at most two constraints, SIAM Journal on Optimization, (2022).  Accepted.
  * [28] M. Song, Y. Xia, and H. Liu, Local optimality conditions for a class of hidden convex optimization, INFORMS Journal on optimization, (2022).  Accepted.
  * [29] E. Stiefel, Richtungsfelder und fernparallelismus in $n$-dimensionalen mannigfaltigkeiten, Comm. Math. Helv., 4 (1935-1936), pp. 305–353.
  * [30] I. Waldspurger and A. Waters, Rank optimality for the Burer-Monteiro factorization, SIAM Journal on Optimization, 30 (2020), pp. 2577–2602.
  * [31] J. Wang and Y. Xia, Closing the gap between necessary and sufficient conditions for local nonglobal minimizer of trust region subproblem, SIAM Journal on Optimization, 30 (2020), pp. 1980–1995.
  * [32] S. Wang and Y. Xia, On the ball-constrained weighted maximin dispersion problem, SIAM Journal on Optimization, 26(3) (2016), pp. 1565–1588.
  * [33] H. Wolkowicz, A note on lack of strong duality for quadratic problems with orthogonal constraints, European Journal of Operational Research, 143 (2002), pp. 356–364.
  * [34] Y. Xia, A survey of hidden convex optimization, Journal of the Operations Research Society of China, 8 (2020), pp. 1–28.
  * [35] Y. Xia and Y. Han, Partial Lagrangian relaxation for the unbalanced orthogonal Procrustes problem, Mathematical Methods of Operations Research, 79 (2014), pp. 225–237.
  * [36] Y. Xue and V. K. N. Lau, Online orthogonal dictionary learning based on Frank-Wolfe method, IEEE Transactions on Neural Networks and Learning Systems, (2021), pp. 1–15.
  * [37] R. Yarlagadda and J. Hershey, Hadamard Matrix Analysis and Synthesis, Springer, Boston, MA, 1997.
  * [38] Y.-x. Yuan, Recent advances in trust region algorithms, Mathematical Programming, 151 (2015), pp. 249–281.
  * [39] F. Zohrizadeh, M. Kheirandishfard, F. Kamangar, and R. Madani, Non-smooth optimization over Stiefel manifolds with applications to dimensionality reduction and graph clustering, in Proceedings of the 28th International Joint Conference on Artificial Intelligence, 2019, pp. 1319–1326.


# Phase Coexistence Implications of Violating Newton’s Third Law

Yu-Jen Chiu  Department of Materials Science and Engineering, University of
California, Berkeley, California 94720, USA  Ahmad K. Omar<EMAIL_ADDRESS>Department of Materials Science and Engineering, University of California,
Berkeley, California 94720, USA  Materials Sciences Division, Lawrence
Berkeley National Laboratory, Berkeley, California 94720, USA

###### Abstract

Newton’s third law, actio = reactio, is a foundational statement of classical
mechanics. However, in natural and living systems, this law appears to be
routinely violated for constituents interacting in a nonequilibrium
environment. Here, we use computer simulations to explore the macroscopic
phase behavior implications of breaking microscopic interaction reciprocity
for a simple model system. We consider a binary mixture of attractive
particles and introduce a parameter which is a continuous measure of the
degree to which interaction reciprocity is broken. In the reciprocal limit,
the species are indistinguishable and the system phase separates into domains
with distinct densities and identical compositions. Increasing nonreciprocity
is found to drive the system to explore a rich assortment of phases, including
phases with strong composition asymmetries and three-phase coexistence. Many
of the states induced by these forces, including traveling crystals and
liquids, have no equilibrium analogue. By mapping the complete phase diagram
for this model system and characterizing these unique phases, our findings
offer a concrete path forward toward understanding how nonreciprocity shapes
the structures found in living systems and how this might be leveraged in
synthetic materials.

Introduction.– Living and natural systems across all length scales can appear
to violate Newton’s third law. A familiar example is the so-called “predator-
prey” interaction whereby one entity (the predator) feels an attractive force
towards the other while the other (the prey) is repelled. These effective
nonreciprocal interactions can emerge from a host of complex factors and can
have far-reaching implications on collective phenomena, phase transitions, and
pattern formation. In the living world, examples range from the interactions
amongst bacteria [1, 2] at the microscale to the dynamics of animal herds [3,
4] at the macroscale. Synthetic systems with nonreciprocal interactions have
also been increasingly realized, with nonreciprocity emerging in systems as
diverse as dusty plasmas [5, 6, 7, 8, 9, 10], colloidal suspensions [11, 12,
13, 14], and even solid metameterials [15, 16, 17]. The understanding of the
structure and phases of biological systems and leveraging nonreciprocal
interactions in synthetic materials will require a fundamental understanding
of both the origins and implications of these effective interactions.

Effective forces acting on particles of interest have long been known to shape
the structures and properties of condensed matter systems. These forces emerge
from the coarse-graining of degrees of freedom, such as other species present
within the system. When the coarse-grained degrees of freedom are in
equilibrium the structure of the effective forces is severely restricted: the
one-body forces (e.g., fluctuating Brownian and dissipative forces acting on
solute upon coarse-graining a molecular solvent) must satisfy the fluctuation-
dissipation theorem (FDT) and many-body forces (e.g., the pairwise depletion
force between larger solutes upon coarse-graining smaller solutes) must be
conservative and consistent with Newton’s third law [18, 19, 20, 21, 22].

Coarse-graining degrees of freedom that are out of equilibrium unlocks a new
range of possibilities. In a nonequilibrium environment, one-body forces need
not satisfy FDT with these “active” forces resulting in novel phase
transitions and collective phenomena [23, 24]. A driven environment may also
generate effective nonreciprocal many-body forces [25, 26, 27, 28]. The
consequences of violating interaction reciprocity on phase behavior remain
unclear (especially in comparison to those of one-body active forces) and have
motivated the development of a variety of theoretical perspectives [25, 27,
29].

In this Letter, we explore the consequences of violating Newton’s third law on
the phase behavior of a model system where noncreciprocity can be continuously
tuned, defining an additional axis on the phase diagram. Increasing
nonreciprocity not only leads to exotic states of coexistence, but may also
stabilize homogeneous states resembling active fluids. Our findings thus
reveal that nonreciprocal interactions may impact the phase behavior in new
and unexpected ways and constitute an important step forward toward the
development of a multicomponent nonequilibrium coexistence theory.

Figure 1: (a) Schematic of pairwise nonreciprocal interactions between species
$L$ and $G$. (b) Phase diagram in the density-nonreciprocity plane with
$\varepsilon/k_{B}T=2.0$ and $\chi=0.5$. (c) Representative snapshots of each
of the seven distinct regions of the phase diagram.

Model system.– Elucidating the influence of nonreciprocal interactions on
material phase behavior requires a model system in which the degree of
nonreciprocity can be continuously varied. We consider a binary mixture of
particles of species $L$ and $G$ with the pairwise force exerted on particle
$i$ by particle $j$ taking the following form [25, 30, 9, 31]:

$\mathbf{F}_{ij}(\mathbf{r})=\mathbf{F}_{ij}^{\rm
C}(\mathbf{r})\times\begin{cases}[1-\Delta(r)]\;&\text{, $ij$ $\in$ $LG$}\\\
[1+\Delta(r)]\;&\text{, $ij$ $\in$ $GL$}\\\ 1&\text{, $ij$ $\in$ $LL$ or
$GG$}\end{cases}$ (1)

where $\mathbf{r}$ is the interparticle distance vector (with magnitude $r$),
$\mathbf{F}_{ij}^{\rm C}(\mathbf{r})$ is a conservative interaction force
satisfying Newton’s third law $\left({\rm i.e.,}~{}\mathbf{F}_{ij}^{\rm
C}+\mathbf{F}_{ji}^{\rm C}=\mathbf{0}\right)$, and $\Delta(r)$ is a
nonreciprocity function. Interactions between particles of like species are
entirely reciprocal $\left(\mathbf{F}_{ij}=\mathbf{F}_{ij}^{\rm C}\right)$
while interspecies interactions exhibit both reciprocal $\mathbf{F}_{ij}^{\rm
R}=\mathbf{F}^{\rm C}_{ij}$ and nonreciprocal $\mathbf{F}_{ij}^{\rm
NR}=\pm\Delta\times\mathbf{F}^{\rm C}_{ij}$ contributions with
$\mathbf{F}_{ij}=\mathbf{F}_{ij}^{\rm R}+\mathbf{F}_{ij}^{\rm NR}$. Consider
$ij$ $\in$ $LG$: the sum of this pair interaction will result in a net force
(violating Newton’s third law) with
$\mathbf{F}_{ij}+\mathbf{F}_{ji}=2\mathbf{F}_{ij}^{\rm
NR}=2\Delta\mathbf{F}^{\rm C}_{ij}$. This resulting net force $\mathbf{F}^{\rm
NR}_{\rm net}=2\mathbf{F}^{\rm NR}_{ij}$ [see Fig. 1(a)] may be thought of as
an internally generated active force that entirely depends on particle
configurations and the nonreciprocity function $\Delta(r)$. For simplicity, we
consider a step nonreciprocity function $\Delta(r)=\Delta\Theta(r-d_{\rm
rec})$ such that interactions between particles are entirely reciprocal for
separations distances less than a reciprocity diameter, $d_{\rm rec}$, and
have a constant degree of nonreciprocity $\Delta$ for distances greater than
$d_{\rm rec}$.

Figure 1(a) illustrates the qualitatively distinct interaction regimes
controlled by $\Delta$. Departures from $\Delta=0$ break interaction
reciprocity, resulting in forces that are no longer conservative and a
nonequilibrium distribution of microstates. For $0<\Delta<1$, while the
magnitude of the force is different for each species, the forces continue to
oppose each other: repulsive (attractive) interactions will continue to be
repulsive (attractive). However, for $\Delta>1$, the nonreciprocity results in
forces that no longer oppose each other: a particle exerting a repulsive
(attractive) force will itself experience an attractive (repulsive) force [see
Fig. 1(a)]. This is the so-called “predator-prey” interaction that appears in
nature.

Our aim is to identify how nonreciprocity shapes the structural and dynamical
landscape of a simple model system. To this end, we take our conservative
force $\mathbf{F}^{C}$ to result from Lennard-Jones (LJ) interactions (with a
cutoff distance of $2.5\sigma$), introducing the LJ energy ($\varepsilon$) and
length ($\sigma$) scales. The reciprocal diameter is set to $d_{\rm
rec}=2^{1/6}\sigma$ such that all particle pairs experience strictly
reciprocal repulsion within separation distances of $d_{\rm rec}$. Particle
dynamics follow the overdamped Langevin equation which imparts an ideal
translational diffusivity $D_{T}=k_{B}T/\zeta$ (introducing our system
timescale $\tau=\sigma^{2}/D_{T}$) where $k_{B}T$ is the thermal energy and
$\zeta$ is the translational drag coefficient. Our system state is thus fully
described by four parameters: the degree of nonreciprocity $\Delta$, the
global volume fraction $\phi=\left(\rho_{L}+\rho_{G}\right)\pi(d_{\rm
rec})^{3}/6$ (where $\rho_{L}$ and $\rho_{G}$ are the number densities of
species $L$ and $G$), the global composition
$\chi=\rho_{G}/(\rho_{L}+\rho_{G})$, and the ratio of the interaction and
thermal energies $\varepsilon/k_{B}T$. In this work, we fix $\chi=0.5$ and
$\varepsilon/k_{B}T=2.0$, and systematically sweep $\phi$ and $\Delta$. All
simulations were conducted with a minimum of $50000$ particles using HOOMD-
blue [32].

Phase diagram.– Figure 1(b) displays the phase diagram obtained from extensive
computer simulations. For finite nonreciprocity, the system is out of
equilibrium, and, while there has been recent progress in the sampling of
driven systems [33, 34, 35, 36, 37, 38, 39, 40, 41], we address questions of
global stability by conducting a number of long-time simulations with distinct
initial configurations, as detailed in the Supplemental Material (SM) [42].

Stationary phase transitions.– In the reciprocal limit, $\Delta=0$, there is
no distinction between $L$ and $G$ particles. The stationary state corresponds
to that of a single-component attractive Lennard-Jones system at
$\varepsilon/k_{B}T=2.0$. While it has been established [43] that the
thermodynamic ground state for this system is crystal (fcc)-fluid coexistence,
observing this state requires an exceedingly rare fluctuation to generate the
critical nuclei necessary for crystal growth. A long-lived (metastable)
liquid-gas coexistence is instead observed for all densities examined in this
work with $\Delta=0$.

Figure 2: Nonreciprocity dependence of coexisting (a) densities and (b)
compositions for states of two-phase coexistence with $\phi=0.0875$.

Small departures from the reciprocal limit ($0<\Delta\leq 0.15$) continue to
result in liquid-gas coexistence with similar (slightly reduced) coexisting
densities [see Fig. 2(a)]. However, a composition asymmetry between the phases
is now generated as shown in Fig. 2(b). For small $\Delta$, it may be
permissible to assume that the system remains in a local equilibrium with
modified interaction potentials between particle pairs. However, a composition
distinction between phases in a symmetric mixture would require difference in
the modified $LL$ and $GG$ interaction energies. As nonreciprocity leaves
interactions between like species unaltered, systems with weak nonreciprocal
forces cannot be mapped to equilibrium with modified interaction energies.

An alternative intriguing possibility exists for this weak species segregation
in which, for small $\Delta$, rather than mapping the system to an equilibrium
system with modified interactions, we consider the effective temperatures of
each species to slightly deviate from the bath temperature and with
$\left(T_{L}\neq T_{G}\right)$ [44]. The species with the higher effective
temperature would be entropically driven to enrich the dilute phase [42].
However, there is little evidence to suggest that $T_{L}>T_{G}$ as the self
diffusion of $G$ is always greater than that of $L$ [see Fig. 3(b)].
Attempting to define the effective temperature with a Stokes-Einstein relation
(i.e., $k_{B}T_{L/G}=D^{\rm self}_{L/G}\zeta^{\rm eff}$ where $\zeta^{\rm
eff}$ is an effective translational drag coefficient assumed to be identical
for each species) would lead to the interpretation that $T_{G}>T_{L}(>T)$. We
thus conclude that even for the small values of nonreciprocity examined here,
effective equilibrium and effective temperature ideas are not sufficient to
describe the observed phase behavior.

Figure 3: Nonreciprocity dependence of the (a) correlation lengths (normalized
by the ideal gas value $\xi_{0}$ [42]) and (b) self diffusivities for
$\phi=0.2$ and $\Delta\geq 0.3$. Shaded regions highlight the active fluid
regime ($D^{\rm self}>D_{T}$). (c) Fluctuations of the local composition
environment of a tagged particle control nonreciprocal forcing.

Increasing the nonreciprocity beyond $\Delta>0.15$ results in a multitude of
states depending on the overall volume fraction. For a broad range of
intermediate densities $\left(0.15<\phi<0.35\right)$ the system is globally
homogeneous [45] (i.e., there is no phase separation) but forms clusters of
finite spatial extent. Each cluster is typically enriched in either $L$ and
$G$, and as a result, the system appears to be microphase separated,
particularly at smaller values of $\Delta$ [cf., Fig. 1(c)]. We characterize
the extent of species segregation by computing the partial static structure
factors, provided in the SM [42]. The correlation lengths for density
fluctuations of species $L$ and $G$, $\xi_{L}$ and $\xi_{G}$, are extracted
from the structure factor and shown in Fig. 3(a). With increasing $\Delta$ we
find a monotonic reduction in these correlation lengths: increasing
nonreciprocity homogenizes the system at these densities. We note that the
divergence of both $\xi_{L}$ and $\xi_{G}$ as $\Delta\rightarrow 0.2$ is
anticipated as the phase-separated regime is approached in this limit.

What’s more, the dramatic alteration to the system’s structural properties is
accompanied with a commensurate change in the species’ self diffusivity [see
Fig. 3(b)]. The diffusion constant in the absence of any interactions is the
ideal Stokes-Einstein translational diffusivity $D_{T}=k_{B}T/\zeta$ imparted
on the particles by the Langevin bath. In the case of passive systems with
purely reciprocal interactions, $D_{T}$ serves as an upper bound for the self-
diffusion constant with $D^{\rm self}\leq D_{T}$. Furthermore, if a Stokes-
Einstein relation holds, $D^{\rm self}=k_{B}T/\zeta^{\rm eff}$ where
$\zeta^{\rm eff}\geq\zeta$. Upon increasing $\Delta$, the self diffusivities
of both species increase significantly and even exceed the ideal diffusivity
for $\Delta\gtrsim 0.6$.

We can determine the precise origins of this enhanced diffusion by using
Green-Kubo relation for self-diffusion [46] and leveraging the substitution of
forces for velocities permitted by overdamped dynamics. Composition
fluctuations about a tagged particle [see Fig. 3(c)] generate instantaneous
nonreciprocal forces that act to increase the diffusivity (much like a one-
body active force) and the system resembles an active fluid. However, this
contribution is found to be relatively independent of $\Delta$, as detailed in
the SM [42]. The origin of the dramatic increase in particle diffusivity with
$\Delta$ is in fact found to be the result of the diminishing magnitude of a
number of force correlations that reduce particle mobility [42]. Increased
nonreciprocal forcing induces profound structural changes that have a clear
dynamical consequence: the breakup of segregated microphases reduces the
effective particle friction. The diminished friction with increasing $\Delta$
coupled with the direct contribution of nonreciprocal forcing to particle
diffusion are what allow the system to exceed the ideal diffusion.

Figure 4: Traveling crystal (a) long time-averaged $\phi$ and $q_{\rm 12}$
profiles (averaging time $\gg\tau_{\rm ts}$) with
$\left(\Delta=1.5,\phi=0.5\right)$ and (b) $\Delta$ dependence of the
traveling state dynamics for $\phi=0.5$.

At volume fractions between regions of global homogeneity [i.e.,
$0.005<\phi\leq 0.15$], nonreciprocity introduces a variety of possible
stationary states. Appreciable values of nonreciprocity
$\left(0.2\leq\Delta\leq 0.55\right)$ results in a state of three-phase
coexistence, shown in Fig. 1(b). The three coexisting fluids consist of two
high density spherical liquid phases [47], each nearly pure in species $L$ or
$G$, and a dilute phase with a relatively equal amount of each species. In
contrast to nonreciprocal two-phase coexistence, the observed densities and
compositions may more readily be mapped to an equilibrium analogue. In
equilibrium, reducing the magnitude of interspecies attraction, i.e.,
$\varepsilon_{LG}<\varepsilon$, can generate a state of three-phase
coexistence not dissimilar to what we observe here, as detailed in the SM
[42].

Further increasing $\Delta$, we return to two-fluid coexistence with densities
similar to those seen at low nonreciprocity [see Fig. 2(a)]. However, unlike
lower values of nonreciprocity, there is a sharp composition contrast between
the coexisting fluids as shown in Fig. 2(b). Intriguingly, the dense fluid
[again adopting a spherical morphology, see Fig. 1(c)] is nearly entirely
comprised of species $G$. Preparing a system seeded with a liquid droplet
purely of species $L$ in this two-phase region of the phase diagram again
results in the nucleation and growth of a liquid $G$ droplet – dense liquids
of $G$ appear to be the preferred state [42]. In single component
nonequilibrium systems, interfacial mechanics uniquely determine the
coexisting densities [48, 49]. While this remains to be generalized to
multicomponent systems, the interfacial nonreciprocal forces result in a pure
$G$ liquid phase coexisting with a dilute mixture as the mechanically stable
state of coexistence.

With increasing nonreciprocity, this region of two-phase coexistence with
strong species segregation is found to narrow while the stable homogeneous
region of the phase diagram expands [see Fig. 1(b)]. At the highest value of
nonreciprocity reported here, the phase diagram has four distinct regions.
Beginning at the lowest density, with increasing $\phi$ we move from the
homogeneous “ideal gas” regime to two-fluid coexistence, followed by a broad
regime of homogeneity. Further increasing the density beyond $\phi\geq 0.35$
results in a traveling state, a dynamical transition which we now discuss.

Traveling states.– At high $\phi$, performing a long-time average results in
density and composition profiles that resemble traditional states of two-phase
coexistence. Intriguingly, at these high concentrations and with increasing
nonreciprocity, the dense phase undergoes an ordering transition with the
emergence of an fcc crystal, quantified with the Steinhardt-Nelson-Ronchetti
order parameter measuring twelve-fold rotational symmetry, $q_{12}$ [50].
Typical density and $q_{12}$ profiles are shown in Fig. 4(a), with $z$ being
the direction normal to the interface. Did the mobility imparted by
nonreciprocal forces allow the system to surmount the nucleation barrier or is
nonreciprocity itself driving this ordering transition? We leave interrogating
these questions for a future study, but recent work on active order-disorder
transitions [51] suggests that nonequilibrium forcing can entirely reshape the
crystallization landscape.

While the long-time averages of these states appear consistent with a
stationary two-phase coexistence scenario, careful examination of the low-
density region reveals a kink in its density profile, suggesting that this is
not a typical state of coexistence between two homogeneous phases. Indeed,
observing the system dynamics reveals the emergence of a traveling state,
defined by ballistic center-of-mass motion (see Fig. 1(c) and SM [42]). This
density gradient generates an asymmetry that drives a persistent center-of-
mass velocity: the coexisting domains appear to chase each other.

These traveling states occur for both liquid-gas and crystal-fluid
coexistence, depending on the degree of nonreciprocity [see Fig. 1(b)]. Over a
time period $\tau_{\rm ts}$, the direction of the low-density phases’ density
gradient reverses, and with it so too does the center-of-mass velocity $v_{\rm
com}$. These dynamics are well-described by a simple oscillatory form ${v_{\rm
com}(t)=v_{\rm ts}\sin(2\pi t/\tau_{\rm ts})}$. Figure 4(c) reveals that while
the velocity increases nearly linearly with nonreciprocity, there is a marked
decrease in the period [52]. That the velocity scales nearly linearly with
$\Delta$ is perhaps indicative that the traveling speed simply scales with the
magnitude of the nonreciprocal forcing while the oscillation period is
dictated by the timescale for the reorganization of the density gradient. The
enhanced mobility with increasing $\Delta$ [cf., Fig. 3(b)] likely drives the
marked reduction in the gradient reorganization time, generating the observed
rapid oscillation frequency at high nonreciprocity.

The density gradient which drives the traveling state is also accompanied by a
composition gradient: the velocity direction is normal to the interface
enriched in $G$ (see Fig. 1(c) and SM movies). These traveling states thus
appear to be a version of the “predator-prey” interaction with the $G$
enriched regions chasing $G$ depleted domains. It is important to note that
traveling states do not require nonreciprocal interactions. Mixtures of active
and passive particles [53, 54, 55, 56], interacting with purely reciprocal
interactions, can also exhibit such states [54, 55, 27]. What determines the
stability of these dynamical states are the mutual diffusion coefficient
tensor $D_{\alpha\beta}$ the elements of which describe how a diffusive flux
of species $\alpha$ emerges from a density gradient in species $\beta$ [57].
However, while $D_{\alpha\beta}$, is crucial for determining stability
criteria, predicting the properties of the stable states (e.g., the density
and composition of phases) will require additional considerations.

Conclusions.– Nonreciprocity is a generic feature to living and natural
systems and is increasingly used as means to alter the structure and dynamics
of synthetic materials across a multitude of length scales. Here,
nonreciprocity is found to generate a wealth of stationary and dynamical phase
transitions, including three-fluid coexistence, liquid-liquid phase separation
with strong species segregation, and traveling crystal phases. The minimal
system presented here may thus serve as an ideal model for the development of
the nonequilibrium statistical mechanics of mixtures.

###### Acknowledgements.

Y.-J.C. acknowledges support from the UC Berkeley College of Engineering Jane
Lewis Fellowship. A.K.O. was supported by the Laboratory Directed Research and
Development Program of Lawrence Berkeley National Laboratory under U.S.
Department of Energy Contract No. DE-AC02-05CH11231 and the UC Berkeley
College of Engineering. This research used the Savio computational cluster
resource provided by the Berkeley Research Computing program.

## References

  * Long and Azam [2001] R. A. Long and F. Azam, Appl. Environ. Microb. 67, 4975 (2001).
  * Czirók  _et al._ [1996] A. Czirók, E. Ben-Jacob, I. Cohen, and T. Vicsek, Phys. Rev. E 54, 1791 (1996).
  * Vicsek and Zafeiris [2012] T. Vicsek and A. Zafeiris, Phys. Rep. 517, 71 (2012).
  * Strandburg-Peshkin  _et al._ [2013] A. Strandburg-Peshkin, C. R. Twomey, N. W. Bode, A. B. Kao, Y. Katz, C. C. Ioannou, S. B. Rosenthal, C. J. Torney, H. S. Wu, S. A. Levin, and I. D. Couzin, Curr. Biol. 23, R709 (2013).
  * Tsytovich [1997] V. N. Tsytovich, Phys.-Usp. 40, 53 (1997).
  * Chaudhuri  _et al._ [2011] M. Chaudhuri, A. V. Ivlev, S. A. Khrapak, H. M. Thomas, and G. E. Morfill, Soft Matter 7, 1287 (2011).
  * Morfill and Ivlev [2009] G. E. Morfill and A. V. Ivlev, Rev. Mod. Phys. 81, 1353 (2009).
  * Keh [2016] H. J. Keh, Curr. Opin. Colloid In. 24, 13 (2016).
  * Ivlev and Kompaneets [2017] A. V. Ivlev and R. Kompaneets, Phys. Rev. E 95, 053202 (2017).
  * Kompaneets  _et al._ [2016] R. Kompaneets, G. E. Morfill, and A. V. Ivlev, Phys. Rev. E 93, 063201 (2016).
  * Meredith  _et al._ [2020] C. H. Meredith, P. G. Moerman, J. Groenewold, Y.-J. Chiu, W. K. Kegel, A. van Blaaderen, and L. D. Zarzar, Nat. Chem. 12, 1136 (2020).
  * Sriram and Furst [2012] I. Sriram and E. M. Furst, Soft Matter 8, 3335 (2012).
  * Soto and Golestanian [2014] R. Soto and R. Golestanian, Phys. Rev. Lett. 112, 068301 (2014).
  * Saha  _et al._ [2019] S. Saha, S. Ramaswamy, and R. Golestanian, New J. Phys. 21, 063006 (2019).
  * Scheibner  _et al._ [2020] C. Scheibner, A. Souslov, D. Banerjee, P. Surówka, W. T. Irvine, and V. Vitelli, Nat. Phys. 16, 475 (2020).
  * Miri and Alù [2019] M. A. Miri and A. Alù, Science 363 (2019).
  * Nassar  _et al._ [2020] H. Nassar, B. Yousefzadeh, R. Fleury, M. Ruzzene, A. Alù, C. Daraio, A. N. Norris, G. Huang, and M. R. Haberman, Nat. Rev. Mater. 5, 667 (2020).
  * Israelachvili [1992] J. N. Israelachvili,  _Intermolecular and Surface Forces_ (Elsevier, Amsterdam, 1992).
  * Dijkstra  _et al._ [2000] M. Dijkstra, R. Van Roij, and R. Evans, J. Chem. Phys. 113, 4799 (2000).
  * Bolhuis  _et al._ [2001] P. G. Bolhuis, A. A. Louis, J. P. Hansen, and E. J. Meijer, J. Chem. Phys. 114, 4296 (2001).
  * Praprotnik  _et al._ [2008] M. Praprotnik, L. D. Site, and K. Kremer, Annu. Rev. Phys. Chem. 59, 545 (2008).
  * Mognetti  _et al._ [2009] B. M. Mognetti, P. Virnau, L. Yelash, W. Paul, K. Binder, M. Müller, and L. G. MacDowell, J. Chem. Phys. 130, 044101 (2009).
  * Cates and Tailleur [2015] M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys. 6, 219 (2015).
  * Bechinger  _et al._ [2016] C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Rev. Mod. Phys. 88, 045006 (2016).
  * Ivlev  _et al._ [2015] A. V. Ivlev, J. Bartnick, M. Heinen, C. R. Du, V. Nosenko, and H. Löwen, Phys. Rev. X 5, 011035 (2015).
  * Hayashi and Sasa [2006] K. Hayashi and S. I. Sasa, J. Phys. Condens. Matter. 18, 2825 (2006).
  * You  _et al._ [2020] Z. You, A. Baskaran, and M. C. Marchetti, Proc. Natl. Acad. Sci. USA 117, 19767 (2020).
  * Durve  _et al._ [2018] M. Durve, A. Saha, and A. Sayeed, Eur. Phys. J. E 41, 1 (2018).
  * Fruchart  _et al._ [2021] M. Fruchart, R. Hanai, P. B. Littlewood, and V. Vitelli, Nature 592, 363 (2021).
  * Bartnick  _et al._ [2015] J. Bartnick, M. Heinen, A. V. Ivlev, and H. Löwen, J. Phys. Condens. Matter 28, 025102 (2015).
  * Kryuchkov  _et al._ [2018] N. P. Kryuchkov, A. V. Ivlev, and S. O. Yurchenko, Soft Matter 14, 9720 (2018).
  * Anderson  _et al._ [2020] J. A. Anderson, J. Glaser, and S. C. Glotzer, Comput. Mater. Sci. 173, 109363 (2020).
  * Klymko  _et al._ [2018] K. Klymko, P. L. Geissler, J. P. Garrahan, and S. Whitelam, Phys. Rev. E 97, 032123 (2018).
  * Ray  _et al._ [2018] U. Ray, G. K.-L. Chan, and D. T. Limmer, J. Chem. Phys. 148, 124120 (2018).
  * Das and Limmer [2019] A. Das and D. T. Limmer, J. Chem. Phys. 151, 244105 (2019).
  * Whitelam  _et al._ [2020] S. Whitelam, D. Jacobson, and I. Tamblyn, J. Chem. Phys. 153, 044113 (2020).
  * Ray and Chan [2020] U. Ray and G. K.-L. Chan, J. Chem. Phys. 152, 104107 (2020).
  * Helms and Chan [2020] P. Helms and G. K.-L. Chan, Phys. Rev. Lett. 125, 140601 (2020).
  * Oakes  _et al._ [2020] T. H. E. Oakes, A. Moss, and J. P. Garrahan, Mach. Learn. Sci. Technol 1, 35004 (2020).
  * Rose  _et al._ [2021] D. C. Rose, J. F. Mair, and J. P. Garrahan, New J. Phys. 23, 013013 (2021).
  * Das  _et al._ [2022] A. Das, B. Kuznets-Speck, and D. T. Limmer, Phys. Rev. Lett. 128, 028005 (2022).
  * Note [1] See Supplemental Material, which includes Refs. [58, 59, 46, 60], for additional simulation and calculation details as well as simulation movies.
  * Schultz and Kofke [2018] A. J. Schultz and D. A. Kofke, J. Chem. Phys 149, 204508 (2018).
  * Note [2] The presence of multiple temperature scales necessarily implies that this system is also out of equilibrium. However, the statistical mechanics of systems with particles in contact with different heat baths is considerably simpler than a full treatment of nonreciprocity.
  * Note [3] This broad region of intermediate densities at which no phase transitions occur $\left(0.15<\phi<0.35\right)$ is not the only region of homogeneity. In the limit of vanishing density $\phi\rightarrow 0$, the system is indistinguishable from an ideal gas with temperature $k_{B}T$. We indeed find the system to be homogeneous at the lowest concentration shown in the phase diagram presented in Fig. 1(b).
  * Hargus  _et al._ [2021] C. Hargus, J. M. Epstein, and K. K. Mandadapu, Phys. Rev. Lett. 127, 178001 (2021).
  * Note [4] The low overall density results in these dense liquids occupying little volume, forming spherical droplets within the dilute fluid phase which comprises most of the system by volume.
  * Aifantis and Serrin [1983] E. C. Aifantis and J. B. Serrin, J. Colloid Interf. Sci. 96, 530 (1983).
  * Omar  _et al._ [2022] A. K. Omar, H. Row, S. A. Mallory, and J. F. Brady, arXiv:2211.12673 (2022).
  * Steinhardt  _et al._ [1983] P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. Rev. B 28, 784 (1983).
  * Omar  _et al._ [2021] A. K. Omar, K. Klymko, T. GrandPre, and P. L. Geissler, Phys. Rev. Lett. 126, 188002 (2021).
  * Note [5] We note that our simulation duration ($5000\tau$) limited the periods we could observe to those with $\tau_{\rm ts}<2500\tau$, which coincides with $\Delta\geq 1.2$.
  * Stenhammar  _et al._ [2015] J. Stenhammar, R. Wittkowski, D. Marenduzzo, and M. E. Cates, Phys. Rev. Lett. 114, 018301 (2015).
  * Wysocki  _et al._ [2016] A. Wysocki, R. G. Winkler, and G. Gompper, New J. Phys. 18, 123030 (2016).
  * Wittkowski  _et al._ [2017] R. Wittkowski, J. Stenhammar, and M. E. Cates, New J. Phys. 19, 105003 (2017).
  * Omar  _et al._ [2019] A. K. Omar, Y. Wu, Z.-G. Wang, and J. F. Brady, ACS Nano 13, 560 (2019).
  * Note [6] The diagonal elements of the mutual diffusion tensor are not to be confused with the self diffusion coefficients computed in this work. The mutual diffusion tensor describes the system response to imposed gradients, and its elements are only equal to self-diffusion coefficients in the dilute limit.
  * Martínez  _et al._ [2009] L. Martínez, R. Andrade, E. G. Birgin, and J. M. Martínez, J. Comput. Chem. 30, 2157 (2009).
  * Stukowski [2009] A. Stukowski, Model. Simul. Mat. Sci. En. 18, 015012 (2009).
  * Weeks  _et al._ [1971] J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys. 54, 5237 (1971).