# Post-Processing with Projection and Rescaling Algorithms for Semidefinite
Programming

Shin-ichi Kanoh and Akiko Yoshise  Graduate School of Systems and Information
Engineering, University of Tsukuba, Tsukuba, Ibaraki 305-8573, and Japan
Society for the Promotion of Science, 5-3-1 Kojimachi, Chiyoda-ku, Tokyo
102-0083, Japan. email<EMAIL_ADDRESS>Corresponding author. Faculty
of Engineering, Information and Systems, University of Tsukuba, Tsukuba,
Ibaraki 305-8573, Japan. email<EMAIL_ADDRESS>
###### Abstract

We propose the algorithm that solves the symmetric cone programs (SCPs) by
iteratively calling the projection and rescaling methods the algorithms for
solving exceptional cases of SCP. Although our algorithm can solve SCPs by
itself, we propose it intending to use it as a post-processing step for
interior point methods since it can solve the problems more efficiently by
using an approximate optimal (interior feasible) solution. We also conduct
numerical experiments to see the numerical performance of the proposed
algorithm when used as a post-processing step of the solvers implementing
interior point methods, using several instances where the symmetric cone is
given by a direct product of positive semidefinite cones. Numerical results
show that our algorithm can obtain approximate optimal solutions more
accurately than the solvers. When at least one of the primal and dual problems
did not have an interior feasible solution, the performance of our algorithm
was slightly reduced in terms of optimality. However, our algorithm stably
returned more accurate solutions than the solvers when the primal and dual
problems had interior feasible solutions.

##  1 Introduction

Let $\mathbb{E}$ be a real-valued vector space with an inner product
$\langle\cdot,\cdot\rangle$. Consider the following symmetric cone programs
(SCPs):

$\begin{array}[]{ccccll}{\rm(P)}&\displaystyle\inf_{x}&\langle c,x\rangle&{\rm
s.t.}&\mathcal{A}x=b,&x\in\mathcal{K},\\\
{\rm(D)}&\displaystyle\sup_{(z,y)}&b^{\top}y&{\rm
s.t.}&z=c-\mathcal{A}^{*}y,&(z,y)\in\mathcal{K}^{*}\times\mathbb{R}^{m},\end{array}$

where $\mathcal{K}\subseteq\mathbb{E}$ is a symmetric cone, $\mathcal{K}^{*}$
is a dual cone of $\mathcal{K}$, i.e.,
$\mathcal{K}^{*}=\\{s\in\mathbb{E}:\langle s,x\rangle\geq 0,\forall
x\in\mathcal{K}\\}$, $\mathcal{A}:\mathbb{E}\rightarrow\mathbb{R}^{m}$ is a
linear operator, $b\in\mathbb{R}^{m}$, $c\in\mathbb{E}$ and
$\mathcal{A}^{*}:\mathbb{R}^{m}\rightarrow\mathbb{E}$ is the adjoint operator
of $\mathcal{A}$, i.e., $\langle\mathcal{A}x,y\rangle=\left\langle
x,\mathcal{A}^{*}y\right\rangle$ for all $x\in\mathbb{E}$ and
$y\in\mathbb{R}^{m}$. Note that $\mathcal{K}^{*}=\mathcal{K}$ holds since
$\mathcal{K}$ is a symmetric cone. If we choose a positive semidefinite cone
as a conic constraint, SCP results in semidefinite programming (SDP). In this
study, we are mainly interested in the case where $\mathcal{K}$ is a Cartesian
product of $p$ simple positive semidefinite cones
$\mathbb{S}^{r_{1}}_{+},\dots,\mathbb{S}^{r_{p}}_{+}$, i.e.,
$\mathcal{K}=\mathbb{S}^{r_{1}}_{+}\times\dots\times\mathbb{S}^{r_{p}}_{+}$.

Interior point methods are among the most popular and practical algorithms to
solve SCPs. In addition to its practical performance, interior point methods
have the theoretical strength of a polynomial-time algorithm and have been
implemented in many solvers. Recently, projection and rescaling methods were
proposed as new polynomial-time algorithms for the special case of SCP by
Lourenço et al. [15], and Kanoh and Yoshise [12]. Their studies were motivated
by a work by Chubanov [5]. A similar type of algorithm, although not a
polynomial-time algorithm, was proposed by Peña and Soheili [20]. The
numerical result of [12] showed that projection and rescaling algorithms
solved ill-conditioned instances, i.e., instances with feasible solutions only
near the boundaries of the cone, more stably than the commercial solver Mosek.
It is well known that in iterations of interior point methods, the closer the
current solution moves to the optimal solution, the more difficult it becomes
to compute the search direction accurately. In other words, interior point
methods cannot work stably near the boundaries of the cone. Therefore, the
numerical results of [12] motivated us to use projection and rescaling
algorithms to solve general SCPs more accurately.

In this study, we propose the algorithms to find approximate optimal solutions
to (P) and (D) using projection and rescaling algorithms as an inner
procedure. Although the proposed methods can find approximate optimal
solutions to (P) and (D) by themselves, they can work more efficiently by
using an approximate optimal (interior feasible) solution. To take advantage
of this feature, the proposed algorithm is designed to be used as a post-
processing step for interior point methods. We also conduct numerical
experiments to compare the performance of our algorithm with the solvers.

The main contribution of this study is to provide comprehensive numerical
results showing whether the projection and rescaling algorithms can be used to
obtain accurate optimal solutions, which will lead to the development of
practical aspects of the projection and rescaling methods. We now mention two
studies [10, 11] focusing on accurately solving SDP. SDPA-GMP [10] is a very
accurate SDP solver that executes the primal-dual interior point method using
multiple-precision arithmetic libraries. Using multiple-precision arithmetic
enhances the accuracy of the solution but, at the same time, increases the
computational time. In [11], Henrion, Naldi and Din proposed an algorithm that
solves Linear matrix inequality (LMI) problems in exact arithmetic, which
means that their algorithm can be used to check the feasibility of dual SDPs.
Although their algorithm can accurately solve small LMI problems, it is
unsuitable for moderate or large problems due to exact computations. This
study differs from these previous studies in that it focuses on the projection
and rescaling methods. Furthermore, this study differs from these studies in
that it proposes an algorithm that can be used in a post-processing step. Our
algorithm is available on the website:
https://github.com/Shinichi-K-4649/Post-processing-using-projection-and-
rescaling-methods.git

We developed our algorithm for SDP. The theoretical foundations of our
algorithm discussed in Section 3 hold for SCP.

###  1.1 Motivation

The primary motivation for focusing on solving SDP accurately is to address
the practical barriers that prevent the implementation of the Facial Reduction
Algorithm (FRA) [4, 3, 19, 32]. The FRA is a regularization technique to
address SDP that does not satisfy Slater’s condition and was developed
originally by Browein and Wolkowicz [4, 3].

We say that (P) satisfies Slater’s condition or (P) is strongly feasible if
there exists a vector $x\in{\rm int}\mathcal{K}$ satisfying $\mathcal{A}x=b$,
where ${\rm int}\mathcal{K}$ is an interior of $\mathcal{K}$. Similarly,
Slater’s condition holds for (D) if there exists $z\in{\rm
int}\mathcal{K}^{*}$ satisfying $c-\mathcal{A}^{*}y=z$. If (P) and/or (D) do
not have the interior feasible solutions, optimal solutions might not exist,
and a duality gap might exist [28, 17, 30], which prevents interior point
methods from working stably.

The FRA transforms SDP with no interior feasible solution into SDP satisfying
Slater’s condition or detects its infeasibility, without changing the optimal
value by iteratively projecting the cone into a specific subspace. Let us
briefly see how the FRA works for (P) here. If (P) has no interior feasible
solution, then there exists a vector $f\in\mathbb{R}^{m}$ such that
$b^{\top}f\geq 0$ and $-\mathcal{A}^{*}f\in\mathcal{K}^{*}\setminus\\{0\\}$,
see Proposition 2.4 in Section 2.1. If $b^{\top}f=0$ holds, $f$ is called a
reducing direction for (P). Otherwise, $f$ is an infeasibility certificate for
(P). The FRA first finds such a vector. If a reducing direction is obtained,
then the FRA generates a new problem with a feasible region equivalent to the
feasible region of (P) by replacing the conic constraint $\mathcal{K}$ with
$\mathcal{K}\cap\\{-\mathcal{A}^{*}f\\}^{\bot}$, where
$\\{-\mathcal{A}^{*}f\\}^{\bot}$ is the space of elements orthogonal to
$-\mathcal{A}^{*}f$. If a reducing direction exists for a new problem, the FRA
tries to find it and generates a new reduced problem again. The FRA repeats
the above process until an infeasibility certificate for a reduced problem or
a strongly feasible problem whose feasible region is equivalent to the
feasible region of (P). Thus, one iteration of the FRA is equivalent to
solving a certain SDP to find a reducing direction or an infeasibility
certificate. For example, implementation at the initial iteration of the FRA
for (P) corresponds to solving the following auxiliary system.

$\begin{array}[]{cccclll}\mbox{(AUX-P)}&{\rm find}\ \ f\in\mathbb{R}^{m}&{\rm
s.t.}&b^{\top}f\geq
0,&-\mathcal{A}^{*}f\in\mathcal{K}^{*}\setminus\\{0\\}.\end{array}$

Therefore, what is required to implement the FRA is to solve the auxiliary
systems exactly. However, this requirement is inaccessible for the current
solvers because the auxiliary systems rarely satisfy Slater’s condition.
Despite such a barrier to implementing the FRA, several studies have conducted
the facial reduction scheme for some problems that arise in applications [36,
33, 13, 31]. These studies directly obtain reducing directions using the
structure of the problem without solving the auxiliary systems numerically.

While these approaches, recently, an implementation of the FRA for any SDP has
been studied [22, 37, 21, 16]. Permenter and Parrilo [22], and Zhu, Pataki and
Tran-Dinh [37] proposed a practical facial reduction scheme. The common idea
of these two studies is that instead of solving an auxiliary system, we solve
another problem whose feasible region is contained in the feasible region of
the original auxiliary system. Although their algorithms might not recover the
interior feasible solutions or detect infeasibility, we can easily implement
them. The algorithms of [22] and [37] can be implemented using only LP solvers
and eigenvalue decomposition, respectively.

Permenter, Friberg and Andersen [21], and Lourenço, Muramatsu and Tsuchiya
[16] proposed an algorithm solving for arbitrary SDP based on a facial
reduction scheme. Permenter, Friberg and Andersen [21] showed that reducing
directions for (P) and (D) can be obtained from relative interior feasible
solutions to a self-dual homogenous model of (P) and (D). In addition, based
on this result, they proposed a theoretical algorithm for solving arbitrary
SDP. We remark that their algorithm finds reducing directions only when
needed, unlike the conventional FRA. Even if (P) and/or (D) do not satisfy
Slater’s condition, their algorithm directly obtains the optimal solution as
long as a complementary solution for (P) and (D) exists. However, their
algorithm requires an oracle that returns relative interior feasible solutions
to the self-dual homogenous model of (P) and (D). Although such solutions are
theoretically obtained using an interior point method tracking a central path,
the numerical experiments of [21] showed practical barriers. On the other
hand, the algorithm of [16] requires only an oracle that returns the optimal
solution to (P) and (D) satisfying Slater’s condition, which is a milder
requirement than [21]. In [16], the authors showed that solving the auxiliary
systems can be replaced by solving a certain SDP with primal and dual interior
feasible points. This study shows that the first step toward fully
implementing the FRA is to solve SDP with primal and dual interior feasible
points as accurately as possible.

One might think this goal is already achievable with the current SDP solvers.
However, the accuracy of the solution that satisfies the termination
tolerances specified by the solvers will not be sufficient for a stable
execution of the FRA. Furthermore, setting the value of the termination
tolerances too small can lead to numerical instability of interior point
methods. Thus, it is worth studying how to solve SDP more accurately than ever
before.

To see practical barriers to implementing the FRA, let us consider a simple
SDP:

######  Example 1.1.

$\begin{array}[]{llc}{\rm(Ex1.1)}\ \ \displaystyle\inf_{x}\
\left\langle\begin{pmatrix}1&0&0\\\ 0&1&0\\\
0&0&1\end{pmatrix},x\right\rangle&{\rm
s.t}&\left\langle\begin{pmatrix}1&1&0\\\ 1&0&0\\\
0&0&0\end{pmatrix},x\right\rangle=1,\\\ &&\left\langle\begin{pmatrix}0&0&0\\\
0&1&0\\\ 0&0&0\end{pmatrix},x\right\rangle=0,\\\
&&\left\langle\begin{pmatrix}0&0&1\\\ 0&0&0\\\
1&0&2\end{pmatrix},x\right\rangle=0,\\\ &&x\in\mathbb{S}^{3}_{+}.\end{array}$

The optimal value of (Ex1.1) is $1$. This problem is feasible but not strongly
feasible. Indeed, any feasible solution $x\in\mathbb{S}^{3}_{+}$ to (Ex1.1)
can be represented as

$x=\begin{pmatrix}1&0&-t\\\ 0&0&0\\\ -t&0&t\end{pmatrix},$

where $t$ is a real value such that $0\leq t\leq 1$. Thus, using a reducing
direction, we can reformulate (Ex1.1) over a lower dimensional positive
semidefinite cone. Any reducing direction $f\in\mathbb{R}^{3}$ for (Ex1.1) can
be represented as $(0,-k,0)^{\top}$, where $k>0$. Indeed, the following holds
for such a vector $f$.

$b^{\top}f=1\times 0+0\times(-k)+0\times 0=0,\ {\rm
and}-\mathcal{A}^{*}f=\begin{pmatrix}0&0&0\\\ 0&k&0\\\
0&0&0\end{pmatrix}\in\mathbb{S}^{3}_{+}.$

Since $-\mathcal{A}^{*}f\in\mathbb{S}^{3}_{+}$, the region
$\mathbb{S}^{3}_{+}\cap\\{-\mathcal{A}^{*}f\\}^{\bot}$ can be represented as
$\mathbb{S}^{3}_{+}\cap\\{-\mathcal{A}^{*}f\\}^{\bot}=U\mathbb{S}^{3-r}_{+}U^{\top}$,
where $r$ is a rank of $-\mathcal{A}^{*}f$ and $U\in\mathbb{R}^{3\times 3-r}$
is a matrix whose columns are the eigenvectors corresponding to the zero
eigenvalues of $-\mathcal{A}^{*}f$. Suppose that $U$ is given by

$U=\begin{pmatrix}1&0\\\ 0&0\\\ 0&1\end{pmatrix}.$

Then, we can generate a reduced problem with interior feasible points as
follows:

$\begin{array}[]{llc}\displaystyle\inf_{\bar{x}}\
\left\langle\begin{pmatrix}1&0\\\ 0&1\end{pmatrix},\bar{x}\right\rangle&{\rm
s.t}&\left\langle\begin{pmatrix}1&0\\\
0&0\end{pmatrix},\bar{x}\right\rangle=1,\\\
&&\left\langle\begin{pmatrix}0&0\\\
0&0\end{pmatrix},\bar{x}\right\rangle=0,\\\
&&\left\langle\begin{pmatrix}0&1\\\
1&2\end{pmatrix},\bar{x}\right\rangle=0,\\\
&&\bar{x}\in\mathbb{S}^{2}_{+}.\end{array}$

Let us apply the FRA to this problem using the commercial solver Mosek. First,
we computed the reducing direction with Mosek by using the formulation
proposed in [16] (see Section 4.3.1). Then, the following vectors were
obtained.

$f_{\rm Mosek}=\begin{pmatrix}\mbox{3.46e-09}\\\ \mbox{-4.00e-00}\\\
\mbox{-1.50e-08}\end{pmatrix}\ {\rm and}\ -\mathcal{A}^{*}f_{\rm
Mosek}=\begin{pmatrix}\mbox{-3.46e-09}&\mbox{-3.46e-09}&\mbox{1.50e-08}\\\
\mbox{-3.46e-09}&\mbox{4.00e-00}&0\\\
\mbox{1.50e-08}&0&\mbox{3.00e-08}\end{pmatrix}.$

The eigenvalue decomposition of $-\mathcal{A}^{*}f_{\rm Mosek}$ was computed
as follows:

$-\mathcal{A}^{*}f_{\rm Mosek}=P_{\rm
Mosek}\begin{pmatrix}\mbox{-9.20e-09}&0&0\\\ 0&\mbox{3.57e-08}&0\\\
0&0&\mbox{4.00e-00}\end{pmatrix}P_{\rm Mosek}^{\top},$

where $P_{\rm Mosek}\in\mathbb{R}^{3\times 3}$ is a matrix whose columns are
the eigenvectors corresponding to the eigenvalues of $-\mathcal{A}^{*}f_{\rm
Mosek}$. Next, we generated a reduced problem. To identify
$\mathbb{S}^{3}_{+}\cap\\{-\mathcal{A}^{*}f_{\rm Mosek}\\}^{\bot}$, for each
eigenvalue of $-\mathcal{A}^{*}f_{\rm Mosek}$, we must correctly determine
whether it is zero or not. In this example, we estimated eigenvalues of
$-\mathcal{A}^{*}f_{\rm Mosek}$ are zero if they are less than a threshold
$T_{\lambda}$. Table 1 shows that the computed optimal value of the reduced
problem is very sensitive to the value of $T_{\lambda}$. If the threshold
$T_{\lambda}$ is set to 1e-7, the rank of $-\mathcal{A}^{*}f_{\rm Mosek}$ is
1, yielding a reduced problem with an optimal value very close to 1. On the
other hand, if the threshold $T_{\lambda}$ is set to 1e-8, the rank of
$-\mathcal{A}^{*}f_{\rm Mosek}$ is 2. Then, solving this reduced problem with
Mosek, we obtained the infeasibility certificate.

Table 1: Comparison of the reduced problems Thresholds | $\mathbb{S}^{3}_{+}\cap\\{-\mathcal{A}^{*}f_{\rm Mosek}\\}^{\bot}$ | Computed optimal value of the reduced problem
---|---|---
$T_{\lambda}=$ 1e-7 | $U\mathbb{S}^{2}_{+}U^{\top}$ | 1
$T_{\lambda}=$ 1e-8 | $U\mathbb{S}^{1}_{+}U^{\top}$ | infeasible

Let us now compute the reducing direction using our algorithm proposed in
Section 4 and the approximate optimal solution returned by Mosek. The computed
reducing direction was as follows.

$f_{\rm Pro}=\begin{pmatrix}\mbox{7.34e-14}\\\ \mbox{-4.00e-00}\\\
\mbox{-4.90e-14}\end{pmatrix}\ {\rm and}\ -\mathcal{A}^{*}f_{\rm
Pro}=\begin{pmatrix}\mbox{-7.34e-14}&\mbox{-7.34e-14}&\mbox{4.90e-14}\\\
\mbox{-7.34e-14}&\mbox{4.00e-00}&0\\\
\mbox{4.90e-14}&0&\mbox{9.80e-14}\end{pmatrix}.$

The eigenvalue decomposition of $-\mathcal{A}^{*}f_{\rm Pro}$ was computed as
follows:

$-\mathcal{A}^{*}f_{\rm Pro}=P_{\rm Pro}\begin{pmatrix}\mbox{-8.67e-14}&0&0\\\
0&\mbox{1.11e-13}&0\\\ 0&0&\mbox{4.00e-00}\end{pmatrix}P_{\rm Pro}^{\top},$

where $P_{\rm Pro}\in\mathbb{R}^{3\times 3}$ is a matrix whose columns are the
eigenvectors corresponding to the eigenvalues of $-\mathcal{A}^{*}f_{\rm
Pro}$. The proposed method found a more accurate solution to the formulation
proposed in [16], resulting in a higher accuracy of the computed reducing
direction. The behavior of the FRA is more stable if the reduced problem is
generated using the vector $f_{\rm Pro}$ rather than the vector $f_{\rm
Mosek}$. From this example, we can see that a more accurate method of solving
SDP is very effective in obtaining accurate reducing directions and that
accurate reducing directions are essential for a stable execution of the FRA.

###  1.2 Structure of this paper

The remainder of this paper is organized as follows. Section 2 is devoted to
some preliminaries and notations. We review the feasibility statuses of SCP
and the projection and rescaling algorithm proposed by Kanoh and Yoshise [12].
Section 3 presents the theoretical foundations of our algorithm and explains
our algorithm. We also describe a practical version of our algorithm that
employs some implementation strategies. Section 4 presents numerical results
showing that our algorithm stably obtains optimal solutions with higher
accuracy than the solvers. The conclusions are summarized in Section 5.

##  2 Preliminaries and Notations

This section describes the preliminaries and notations used throughout this
paper. Section 2.1 provides an overview of SCP, focusing on the feasibility
statuses. Section 2.2 briefly introduces Euclidean Jordan algebras. Although
our algorithm uploaded on the website is an algorithm for solving SDP, its
theoretical foundations can be extended to SCP without problems. The algorithm
for SCP can be concisely described using Euclidean Jordan algebras. Thus, we
summarize the minimum knowledge of Euclidean Jordan algebras required to
understand this study. Section 2.3 briefly describes the projection and
rescaling algorithm proposed by [12], which is used in our algorithm. Our
algorithm employed the algorithm of [12] because the numerical results in [12]
showed that their algorithm was superior to the other methods in terms of
computation time.

###  2.1 Feasibility classes of SCP

Let $\theta_{p}\in\mathbb{R}\cup\\{\pm\infty\\}$ and
$\theta_{d}\in\mathbb{R}\cup\\{\pm\infty\\}$ be the primal and dual optimal
values, respectively. If (P) is infeasible, $\theta_{p}$ takes $+\infty$ and
if (D) is infeasible, $\theta_{d}$ takes $-\infty$. By
$C^{\bot}\subseteq\mathbb{E}$, we denote the space of elements orthogonal to a
set $C\subseteq\mathbb{E}$. Here, we review four different mutually exclusive
feasibility classes of SCP. These feasibility classes are defined in the field
of conic linear programming. Therefore, in this section, we purposely leave
$\mathcal{K}^{*}$ as it is. First, let us look at the definition of strong
feasibility.

######  Definition 2.1 (Strongly feasible).

We say that (P) is strongly feasible (or satisfies Slater’s condition) if
there exists $x\in\mbox{int }\mathcal{K}$ satisfying $\mathcal{A}x=b$.
Similarly, we say that (D) is strongly feasible (or satisfies Slater’s
condition) if there exists $y\in\mathbb{R}^{m}$ satisfying
$c-\mathcal{A}^{*}y\in\mbox{int }\mathcal{K}^{*}$.

If primal and dual problems are strongly feasible, $\theta_{p}=\theta_{d}$ and
the existence of optimal solutions to both problems are guaranteed.

######  Proposition 2.2.

  1. 1.

If (P) is strongly feasible, then $\theta_{p}=\theta_{d}$. In addition, if (D)
is feasible, then (D) has an optimal solution.

  2. 2.

If (D) is strongly feasible, then $\theta_{p}=\theta_{d}$. In addition, if (P)
is feasible, then (P) has an optimal solution.

###### Proof.

See Theorems 3.2.6 and 3.2.8 in [23]. ∎

We next see the characterization of strong infeasibility.

######  Definition 2.3 (Strongly infeasible).

We say that (P) is strongly infeasible if there exists $y\in\mathbb{R}^{m}$
such that $-\mathcal{A}^{*}y\in\mathcal{K}^{*}$ and $b^{\top}y>0$. Similarly,
we say that (D) is strongly infeasible if there exists $x\in\mathcal{K}$ such
that $\mathcal{A}x=0$ and $\langle c,x\rangle<0$.

A vector $y\in\mathbb{R}^{m}$ satisfying $-\mathcal{A}^{*}y\in\mathcal{K}^{*}$
and $b^{\top}y>0$ is called an improving ray of (D) because $y$ makes (D)
unbounded, i.e., $\theta_{d}=+\infty$, if (D) is feasible. We also call
$x\in\mathcal{K}$ an improving ray of (P) if $\mathcal{A}x=0$ and $\langle
c,x\rangle<0$. If there exists an improving ray of (P) (or (D)), then we can
see that the hyperplane $\\{\mathcal{A}^{*}y\\}^{\bot}$ (or $\\{x\\}^{\bot}$)
strictly separates the affine space of (P) (or (D)) and $\mathcal{K}$ (or
$\mathcal{K}^{*}$), which implies the infeasibility of (P) (or (D)).

To define the remaining feasibility class, we introduce the following
proposition.

######  Proposition 2.4.

  1. 1.

(P) is not strongly feasible if and only if there exists $y\in\mathbb{R}^{m}$
such that $-\mathcal{A}^{*}y\in\mathcal{K}^{*}$, $\mathcal{A}^{*}y\neq 0$ and
$b^{\top}y\geq 0$.

  2. 2.

(D) is not strongly feasible if and only if there exists a nonzero
$x\in\mathcal{K}$ such that $\mathcal{A}x=0$ and $\langle c,x\rangle\leq 0$.

###### Proof.

See Lemma 3.2 in [32] or Theorems 3.1.3 and 3.1.4 in [6]. ∎

Here, we define the following feasibility class.

######  Definition 2.5 (Weak status).

We say that (P) is in weak status if there exists $y\in\mathbb{R}^{m}$ such
that $-\mathcal{A}^{*}y\in\mathcal{K}^{*}$, $\mathcal{A}^{*}y\neq 0$ and
$b^{\top}y=0$. Similarly, we say that (D) is in weak status if there exists a
nonzero $x\in\mathcal{K}$ such that $\mathcal{A}x=0$ and $\langle
c,x\rangle=0$.

We sometimes divide this class into two classes and call them weakly feasible
and weakly infeasible, respectively. It is understood that a problem is weakly
feasible if it is feasible and does not have interior feasible solutions.
Similarly, it is understood that a problem is weakly infeasible if it is
infeasible and its dual does not have an improving ray. See [17] for more
details.

A vector $y\in\mathbb{R}^{m}$ (or a vector $-\mathcal{A}^{*}y\in\mathbb{E}$)
satisfying $-\mathcal{A}^{*}y\in\mathcal{K}^{*}$, $\mathcal{A}^{*}y\neq 0$ and
$b^{\top}y=0$ is called a reducing direction for (P). We also call a nonzero
$x\in\mathcal{K}$ a reducing direction for (D) if $\mathcal{A}x=0$ and
$\langle c,x\rangle=0$. As discussed in Section 1.1, reducing directions are
used to regularize problems in the FRA.

###  2.2 Euclidean Jordan algebras and Symmetric cones

Let $\mathbb{E}$ be a real-valued vector space equipped with an inner product
$\langle\cdot,\cdot\rangle$ and a bilinear operation $\circ$ :
$\mathbb{E}\times\mathbb{E}\rightarrow\mathbb{E}$, and $e$ be the identity
element, i.e.,$x\circ e=e\circ x=x$ holds for any $x\in\mathbb{E}$.
$(\mathbb{E},\circ)$ is called a Euclidean Jordan algebra if it satisfies

$\displaystyle x\circ y=y\circ x,\ \ x\circ(x^{2}\circ y)=x^{2}\circ(x\circ
y),\ \ \langle x\circ y,z\rangle=\langle y,x\circ z\rangle$

for all $x,y,z\in\mathbb{E}$ and $x^{2}:=x\circ x$. We denote $y\in\mathbb{E}$
as $x^{-1}$ if $y$ satisfies $x\circ y=e$. $c\in\mathbb{E}$ is called an
idempotent if it satisfies $c\circ c=c$, and an idempotent $c$ is called
primitive if it can not be written as a sum of two or more nonzero
idempotents. A set of primitive idempotents $c_{1},c_{2},\ldots c_{k}$ is
called a Jordan frame if $c_{1},\ldots c_{k}$ satisfy

$c_{i}\circ c_{j}=0\ (i\neq j),\ \ c_{i}\circ c_{i}=c_{i}\ (i=1,\dots,k),\ \
\sum_{i=1}^{k}c_{i}=e.$

For $x\in\mathbb{E}$, the degree of $x$ is the smallest integer $d$ such that
the set $\\{e,x,x^{2},\ldots,x^{d}\\}$ is linearly independent. The rank of
$\mathbb{E}$ is the maximum integer $r$ of the degree of $x$ over all
$x\in\mathbb{E}$. The following properties are known.

######  Proposition 2.6 (Spectral theorem (cf. Theorem III.1.2 of [7])).

Let $(\mathbb{E},\circ)$ be a Euclidean Jordan algebra having rank $r$. For
any $x\in\mathbb{E}$, there exist real numbers $\lambda_{1},\dots,\lambda_{r}$
and a Jordan frame $c_{1},\dots,c_{r}$ for which the following holds:

$x=\sum_{i=1}^{r}\lambda_{i}c_{i}.$

The numbers $\lambda_{1},\dots,\lambda_{r}$ are uniquely determined
eigenvalues of $x$ (with their multiplicities). Furthermore, ${\rm
trace}(x):=\sum_{i=1}^{r}\lambda_{i}$, $\det(x):=\prod_{i=1}^{r}\lambda_{i}$.

For any $x,y\in\mathbb{E}$, we define the inner product
$\langle\cdot,\cdot\rangle$ and the norm $\|\cdot\|$ as $\langle
x,y\rangle:={\rm trace}(x\circ y)$ and $\|x\|:=\sqrt{\langle x,x\rangle}$,
respectively. For any $x\in\mathbb{E}$ having spectral decomposition
$x=\sum_{i=1}^{r}\lambda_{i}c_{i}$ as in Proposition 2.6, we also define
$\|x\|_{\infty}:=\max\\{|\lambda_{1}|,\dots,|\lambda_{r}|\\}$.

It is known that the set of squares $\mathcal{K}=\\{x^{2}:x\in\mathbb{E}\\}$
is the symmetric cone of $\mathbb{E}$ (cf. Theorems III.2.1 and III.3.1 of
[7]). The following properties can be derived from the results in [7], as in
Corollary 2.3 of [35]:

######  Proposition 2.7.

Let $x\in\mathbb{E}$ and let $\sum_{j=1}^{r}\lambda_{j}c_{j}$ be a
decomposition of $x$ given by Propositoin 2.6. Then

(i)

$x\in\mathcal{K}$ if and only if $\lambda_{j}\geq 0\ (j=1,2,\ldots,r)$,

(ii)

$x\in{\rm int}\hskip 2.13394pt\mathcal{K}$ if and only if $\lambda_{j}>0\
(j=1,2,\ldots,r)$.

A Euclidean Jordan algebra $(\mathbb{E},\circ)$ is called simple if it cannot
be written as any Cartesian product of non-zero Euclidean Jordan algebras. If
the Euclidean Jordan algebra $(\mathbb{E},\circ)$ associated with a symmetric
cone $\mathcal{K}$ is simple, then we say that $\mathcal{K}$ is simple. In
this study, we will consider that $\mathcal{K}$ is given by a Cartesian
product of $p$ simple positive semidefinite cones $\mathbb{S}^{r_{\ell}}_{+}$,
$\mathcal{K}:=\mathbb{S}^{r_{1}}_{+}\times\dots\times\mathbb{S}^{r_{p}}_{+}$,
whose rank and identity element are $r_{\ell}$ and $e_{\ell}$
$(\ell=1,\ldots,p)$. The rank $r$ and the identity element of $\mathcal{K}$
are given by

$r=\sum_{\ell=1}^{p}r_{\ell},\ \ e=(e_{1},\dots,e_{p}).$

Furthermore, for any $x\in\mathbb{E}_{1}\times\dots\times\mathbb{E}_{p}$,
$\det(x)=\prod_{\ell=1}^{p}\prod_{i=1}^{r_{\ell}}\lambda_{i}$.

In what follows, $x_{\ell}$ stands for the $\ell$-th block element of
$x\in\mathbb{E}$, i.e.,
$x=(x_{1},\dots,x_{p})\in\mathbb{E}_{1}\times\dots\times\mathbb{E}_{p}$. For
each $\ell=1,\cdots,p$, we define
$\lambda_{\min}(x_{\ell}):=\min\\{\lambda_{1},\cdots,\lambda_{r_{\ell}}\\}$
and
$\lambda_{\max}(x_{\ell}):=\max\\{\lambda_{1},\cdots,\lambda_{r_{\ell}}\\}$
where $\lambda_{1},\cdots,\lambda_{r_{\ell}}$ are eigenvalues of $x_{\ell}$.
The minimum and maximum eigenvalues of $x\in\mathcal{K}$ are given by
$\lambda_{\min}(x)=\min\\{\lambda_{\min}(x_{1}),\cdots,\lambda_{\min}(x_{p})\\}$
and
$\lambda_{\max}(x)=\max\\{\lambda_{\max}(x_{1}),\cdots,\lambda_{\max}(x_{p})\\}$,
respectively.

Next, we consider the quadratic representation $Q_{v}(x)$ defined by
$Q_{v}(x):=2v\circ(v\circ x)-v^{2}\circ x$. For the Euclidean Jordan algebra
$(\mathbb{E},\circ)$ such as
$\mathbb{E}=\mathbb{E}_{1}\times\dots\times\mathbb{E}_{p}$, the quadratic
representation $Q_{v}(x)$ of $x\in\mathbb{E}$ is denoted by
$Q_{v}(x)=\left(Q_{v_{1}}(x_{1}),\dots,Q_{v_{p}}(x_{p})\right)$. Letting
$\mathcal{I}_{\ell}$ be the identity operator of the Euclidean Jordan algebra
$(\mathbb{E}_{\ell},\circ_{\ell})$ associated with the cone
$\mathcal{K}_{\ell}$, we have $Q_{e_{\ell}}=\mathcal{I}_{\ell}$ for
$\ell=1,\ldots,p$. The following properties can also be retrieved from the
results in [7] as in Proposition 3 of [15]:

######  Proposition 2.8.

For any $v\in{\rm int}\mathcal{K}$, $Q_{v}(\mathcal{K})=\mathcal{K}$.

More detailed descriptions, including concrete examples of symmetric cone
optimization, can be found in, e.g., [7, 8, 25]. Here, we will explain the
bilinear operation, the identity element, the inner product, the eigenvalues,
the primitive idempotents, and the quadratic representation of the cone when
the cone is a positive semidefinite cone.

######  Example 2.9 ($\mathcal{K}$ is the semidefinite cone
$\mathbb{S}^{n}_{+}$).

Let $\mathbb{S}^{n}$ be the set of symmetric matrices of $n\times n$. The
semidefinite cone $\mathbb{S}^{n}_{+}$ is given by
$\mathbb{S}^{n}_{+}=\\{X\in\mathbb{S}^{n}:X\succeq O\\}$. For any symmetric
matrices $X,Y\in\mathbb{S}^{n}$, define the bilinear operation $\circ$ and
inner product as $X\circ Y=\frac{XY+YX}{2}$ and $\langle
X,Y\rangle=\mbox{tr}(XY)=\sum_{i=1}^{n}\sum_{j=1}^{n}X_{ij}Y_{ij}$,
respectively. For any $X\in\mathbb{S}^{n}$, perform the eigenvalue
decomposition and let $u_{1},\dots,u_{n}$ be the corresponding normalized
eigenvectors for the eigenvalues $\lambda_{1},\dots,\lambda_{n}$:
$X=\sum_{i=1}^{n}\lambda_{i}u_{i}u_{i}^{T}$. The eigenvalues of $X$ in the
Jordan algebra are $\lambda_{1},\dots,\lambda_{n}$ and the primitive
idempotents are $c_{1}=u_{1}u_{1}^{T},\dots,c_{n}=u_{n}u_{n}^{T}$, which
implies that the rank of the semidefinite cone $\mathbb{S}^{n}_{+}$ is $r=n$.
The identity element is the identity matrix $I$. The quadratic representation
of $V\in\mathbb{S}^{n}$ is given by $Q_{V}(X)=VXV$.

###  2.3 Projection and rescaling algorithm

For two sets $C_{1}$ and $C_{2}$, we denote by $\mbox{FP}(C_{1},C_{2})$ the
feasibility problem

$\begin{array}[]{ll}\mbox{find}&x\in C_{1}\cap C_{2}.\end{array}$

The projection and rescaling algorithms [20, 15, 12] solve
$\mbox{FP}(\mathcal{L},\mbox{int }\mathcal{K})$, where
$\mathcal{K}\subseteq\mathbb{E}$ is a symmetric cone and
$\mathcal{L}\subseteq\mathbb{E}$ is a linear subspace. The feasibility of
$\mbox{FP}(\mathcal{L},\mbox{int }\mathcal{K})$ is closely related to the
feasibility of $\mbox{FP}(\mathcal{L}^{\bot},\mathcal{K}\setminus\\{0\\})$,
where $\mathcal{L}^{\bot}$ is the orthogonal complement of $\mathcal{L}$. The
proof of Proposition 2.10 is straightforward using Theorem 20.2 of [24] and
therefore omitted.

######  Proposition 2.10.

$\mbox{FP}(\mathcal{L},\mbox{int }\mathcal{K})$ is infeasible if and only if
$\mbox{FP}(\mathcal{L}^{\bot},\mathcal{K}\setminus\\{0\\})$ is feasible.

For any $k>0$ and feasible solution $x$ of $\mbox{FP}(\mathcal{L},\mbox{int
}\mathcal{K})$, since $kx$ is a feasible solution of
$\mbox{FP}(\mathcal{L},\mbox{int }\mathcal{K})$, it makes sense to consider
the positive scaled version of $\mbox{FP}(\mathcal{L},\mbox{int
}\mathcal{K})$. In [12], they consider the following feasibility problem. We
will denote it by $\mbox{FP}_{S_{\infty}}(\mathcal{L},\mbox{int }\mathcal{K})$
in this study.

$\begin{array}[]{ll}\mbox{find}&x\in\mathcal{L}\cap\mbox{int }\mathcal{K}\ \ \
\mbox{s.t.}\ \ \|x\|_{\infty}\leq 1.\end{array}$

######  Proposition 2.11.

$\mbox{FP}(\mathcal{L},\mbox{int }\mathcal{K})$ is feasible if and only if
$\mbox{FP}_{S_{\infty}}(\mathcal{L},\mbox{int }\mathcal{K})$ is feasible.

###### Proof.

If $\mbox{FP}_{S_{\infty}}(\mathcal{L},\mbox{int }\mathcal{K})$ is feasible,
it is clear that $\mbox{FP}(\mathcal{L},\mbox{int }\mathcal{K})$ is feasible.

Let $x$ be a feasible solution of $\mbox{FP}(\mathcal{L},\mbox{int
}\mathcal{K})$ and let $\lambda_{\max}(x)$ and $\lambda_{\min}(x)$ be a
maximum eigenvalue and a minimum eigenvalue of $x$, respectively. Since
$x\in\mbox{int }\mathcal{K}$, we can see that
$\lambda_{\max}(x)\geq\lambda_{\min}(x)>0$ and hence
$\frac{1}{\lambda_{\max}(x)}x\in\mbox{int }\mathcal{K}$. In addition,
$\|\frac{1}{\lambda_{\max}(x)}x\|_{\infty}=1$ holds. Thus,
$\frac{1}{\lambda_{\max}(x)}x$ is a feasible solution of
$\mbox{FP}_{S_{\infty}}(\mathcal{L},\mbox{int }\mathcal{K})$. ∎

The projection and rescaling algorithms consist of two ingredients: the “main
algorithm” and the “basic procedure.” The structure of the method is as
follows: In the outer iteration, the main algorithm calls the basic procedure
with $\mathcal{L}$ and $\mathcal{K}$. The basic procedure proposed in [12]
generates a sequence in $\mathbb{E}$ using projection to $\mathcal{L}$ and
terminates in a finite number of iterations returning one of the following:
(i). a solution of problem $\mbox{FP}(\mathcal{L},\mbox{int }\mathcal{K})$,
(ii). a solution of problem
$\mbox{FP}(\mathcal{L}^{\bot},\mathcal{K}\setminus\\{0\\})$, or (iii). a cut
of $\mbox{FP}_{S_{\infty}}(\mathcal{L},\mbox{int }\mathcal{K})$, i.e., a
Jordan frame $\\{c_{1},c_{2},\dots,c_{r}\\}$ such that $\langle
c_{i},x\rangle\leq\xi$ holds for any feasible solution $x$ of problem
$\mbox{FP}_{S_{\infty}}(\mathcal{L},\mbox{int }\mathcal{K})$ and for some
$i\in\\{1,2,\dots,r\\}$, where $r$ is a rank of $\mathcal{K}$ and $\xi$ is a
parameter specified by the user such that $0<\xi<1$. If the result (i) or (ii)
is returned by the basic procedure, then the feasibility of problem
$\mbox{FP}(\mathcal{L},\mbox{int }\mathcal{K})$ can be determined, and the
main procedure stops. If the result (iii) is returned, then the main procedure
scales the problem $\mbox{FP}_{S_{\infty}}(\mathcal{L},\mbox{int
}\mathcal{K})$ as $\mbox{FP}_{S_{\infty}}(Q_{v}(\mathcal{L}),\mbox{int
}\mathcal{K})$, where $v=\frac{1}{\sqrt{\xi}}\sum_{h\in H}c_{h}+\sum_{h\notin
H}c_{h}$ and $H=\\{i:\langle c_{i},x\rangle\leq\xi\\}$. Then, the main
procedure calls the basic procedure with $Q_{v}(\mathcal{L})$ and
$\mathcal{K}$. Noting that $v\in{\rm int}\mathcal{K}$ and Proposition 2.8, the
feasibility of $\mbox{FP}(\mathcal{L},\mbox{int }\mathcal{K})$ can be checked
by solving $\mbox{FP}_{S_{\infty}}(Q_{v}(\mathcal{L}),\mbox{int
}\mathcal{K})$. Thus, the projection and rescaling algorithm of [12] checks
the feasibility of $\mbox{FP}(\mathcal{L},\mbox{int }\mathcal{K})$ by
repeating the above procedures.

Their algorithm has a feature that the main algorithm works while keeping
information about the minimum eigenvalues of any feasible solution of
$\mbox{FP}_{S_{\infty}}(\mathcal{L},\mbox{int }\mathcal{K})$. That is, their
method can determine whether there exists a feasible solution of
$\mbox{FP}_{S_{\infty}}(\mathcal{L},\mbox{int }\mathcal{K})$ whose minimum
eigenvalue is greater than $\varepsilon$, where $\varepsilon>0$ is a parameter
specified by the user. In [12], a feasible solution of
$\mbox{FP}_{S_{\infty}}(\mathcal{L},\mbox{int }\mathcal{K})$ whose minimum
eigenvalue is greater than or equal to $\varepsilon$ is called an
$\varepsilon$-feasible solution of
$\mbox{FP}_{S_{\infty}}(\mathcal{L},\mbox{int }\mathcal{K})$.

Algorithm 1 The projection and rescaling algorithm of [12]

1: Input: $\mathcal{L}$, $\mathcal{K}$, $\varepsilon>0$ and a constant $\xi$
such that $0<\xi<1$.

2: Output: A solution to $\mbox{FP}(\mathcal{L},\mbox{int }\mathcal{K})$ or
$\mbox{FP}(\mathcal{L}^{\bot},\mathcal{K}\setminus\\{0\\})$ or a certificate
that there is no $\varepsilon$ feasible solution to problem
$\mbox{FP}_{S_{\infty}}(\mathcal{L},\mbox{int }\mathcal{K})$.

3: initialization: $k\leftarrow 1,\ \mathcal{L}^{k}\leftarrow\mathcal{L}$

4: Call the basic procedure with $\mathcal{L}^{k}$, $\mathcal{K}$ and $\xi$.

5: if a solution to $\mbox{FP}(\mathcal{L}^{k},\mbox{int }\mathcal{K})$ is
obtained then

6: Rescale the obtained solution to the solution of
$\mbox{FP}(\mathcal{L},\mbox{int }\mathcal{K})$

7: return the solution to $\mbox{FP}(\mathcal{L},\mbox{int }\mathcal{K})$

8: else if a solution to
$\mbox{FP}({\mathcal{L}^{k}}^{\bot},\mathcal{K}\setminus\\{0\\})$ is obtained
then

9: Rescale the obtained solution to the solution of
$\mbox{FP}(\mathcal{L}^{\bot},\mathcal{K}\setminus\\{0\\})$

10: return the solution to
$\mbox{FP}(\mathcal{L}^{\bot},\mathcal{K}\setminus\\{0\\})$

11: else

12: Compute an upper bound for the minimum eigenvalue of any feasible solution
of $\mbox{FP}_{S_{\infty}}(\mathcal{L},\mbox{int }\mathcal{K})$

13: if the computed upper bound is less than $\varepsilon$ then

14: stop Algorithm 1 (There is no $\varepsilon$ feasible solution to
$\mbox{FP}_{S_{\infty}}(\mathcal{L},\mbox{int }\mathcal{K})$.)

15: end if

16: Compute the vector $v$ used to scale the problem

17: end if

18: Scale the linear subspace, i.e., $\mathcal{L}^{k+1}\leftarrow
Q_{v}(\mathcal{L}^{k})$

19: $k\leftarrow k+1$. Go back to line 4.

##  3 Proposed algorithms

Since projection and rescaling algorithms can only solve the special case of
SCP, i.e., $\mbox{FP}(\mathcal{L},\mbox{int }\mathcal{K})$ for a linear
subspace $\mathcal{L}\subseteq\mathbb{E}$, we consider how to use projection
and rescaling algorithms to obtain an approximate optimal solution to (P) or
(D). In Section 3.1, we introduce two types of formulations to which
projection and rescaling methods can be applied. These formulations require a
real value $\theta\in\mathbb{R}$, and their feasibilities depend on the value
of $\theta$. In Section 3.1.1, we introduce the formulation that gives the
interior feasible solution $x$ of (P) such that $\langle c,x\rangle<\theta$
when (P) is strongly feasible and $\theta>\theta_{p}$. In Section 3.1.2, we
introduce the formulation that gives the interior feasible solution $(y,z)$ of
(D) such that $b^{\top}y>\theta$ when (D) is strongly feasible and
$\theta<\theta_{d}$. Then, we present a basic idea of our algorithm in Section
3.2. We also discuss some implementation strategies of our algorithm in
Section 3.3. We then develop a practical version of our algorithm in Section
3.4.

###  3.1 Theoretical foundations

####  3.1.1 Primal model

Let $\theta\in\mathbb{R}$. For $\theta$, let us define the linear operator
$\mathcal{A}(\theta)$ as follows:

$\mathcal{A}(\theta)=\begin{pmatrix}\mathcal{A}&-b&\bm{0}\\\
c^{\top}&-\theta&1\end{pmatrix}.$

We define the symmetric cone
$\bar{\mathcal{K}}=\mathcal{K}\times\mathbb{R}^{2}_{+}$ and consider the
feasibility problem $\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\mbox{int
}\bar{\mathcal{K}})$.

######  Proposition 3.1.

Suppose that $\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\mbox{int
}\bar{\mathcal{K}})$ is feasible and $(x,\tau,\rho)$ is a feasible solution of
it. Then, $\frac{1}{\tau}x$ is an interior feasible solution to (P), and
$\langle c,\frac{1}{\tau}x\rangle<\theta$ holds.

###### Proof.

Since $(x,\tau,\rho)$ is a feasible solution of
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\mbox{int }\bar{\mathcal{K}})$,
$\mathcal{A}x-\tau b=0$ and $\langle c,x\rangle-\tau\theta+\rho=0$ hold.
Noting that $\tau>0$, we have $\mathcal{A}\frac{1}{\tau}x-b=0$ and $\langle
c,\frac{1}{\tau}x\rangle-\theta+\frac{\rho}{\tau}=0$ and hence,
$\frac{1}{\tau}x$ is an interior feasible solution of (P) such that $\langle
c,\frac{1}{\tau}x\rangle<\theta$ because $x\in\mbox{int }\mathcal{K}$ and
$\rho>0$. ∎

From Proposition 3.1, if we obtain the feasible solution of
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\mbox{int }\bar{\mathcal{K}})$, then
we can construct the interior feasible solution to (P) whose objective value
is smaller than $\theta$. The following proposition gives us a necessary and
sufficient condition for $\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\mbox{int
}\bar{\mathcal{K}})$ to be feasible.

######  Proposition 3.2.

$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\mbox{int }\bar{\mathcal{K}})$ is
feasible if and only if (P) is strongly feasible and $\theta_{p}<\theta$.

###### Proof.

If $\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\mbox{int }\bar{\mathcal{K}})$ is
feasible, then there exists a feasible solution $(x,\tau,\rho)$ to
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\mbox{int }\bar{\mathcal{K}})$. For
$x$ and $\tau$, $\frac{1}{\tau}x$ is an interior feasible solution to (P)
satisfying $\langle c,\frac{1}{\tau}x\rangle<\theta$ from Proposition 3.1,
which implies that (P) has an interior feasible solution and
$\theta_{p}\leq\langle c,\frac{1}{\tau}x\rangle<\theta$.

Conversely, if (P) is stronlgy feasible and $\theta_{p}<\theta$, then there
exists an interior feasible solution $x$ to (P) such that $\theta_{p}<\langle
c,x\rangle<\theta$. We can easily see that $(x,1,\theta-\langle c,x\rangle)$
is a feasible solution for $\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\mbox{int
}\bar{\mathcal{K}})$. ∎

Combining Proposition 3.2 with Proposition 2.10, we have a necessary and
sufficient condition for alternative problem
$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta)}^{*},\bar{\mathcal{K}}\setminus\\{0\\})$
of $\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\mbox{int }\bar{\mathcal{K}})$ to
be feasible.

######  Corollary 3.3.

$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta)}^{*},\bar{\mathcal{K}}\setminus\\{0\\})$
is feasible if and only if (P) is not strongly feasible or
$\theta_{p}\geq\theta$.

While feasible solutions of $\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\mbox{int
}\bar{\mathcal{K}})$ give us interior feasible solutions to (P) whose
objective value is smaller than $\theta$, feasible solutions of
$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta)}^{*},\bar{\mathcal{K}}\setminus\\{0\\})$
give us information about the feasibility of (P) or a feasible solution to (D)
whose objective value is greater than or equal to $\theta$.

######  Proposition 3.4.

Suppose that
$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta)}^{*},\bar{\mathcal{K}}\setminus\\{0\\})$
is feasible and $(z,\omega,\kappa)$ is a feasible solution of
$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta)}^{*},\bar{\mathcal{K}}\setminus\\{0\\})$.
Then, there exists $y\in\mathbb{R}^{m}$ such that $z=\mathcal{A}^{*}y+\kappa
c$ and $\omega=-b^{\top}y-\kappa\theta$. For $(y,\kappa)$, one of the
following three cases holds:

  1. 1.

$\kappa>0$ meaning that
$(-\frac{1}{\kappa}y,c-\mathcal{A}^{*}(\frac{-1}{\kappa}y))$ is a feasible
solution to (D) and its objective value is greater than or equal to $\theta$,

  2. 2.

$\kappa=0$ and $\omega>0$ meaning that $-y$ is an improving ray of (D), i.e.,
$-\mathcal{A}^{*}(-y)\in\mathcal{K}$ and $b^{\top}(-y)>0$, or

  3. 3.

$\kappa=\omega=0$ meaning that $-y$ is a reucing direction for (P), i.e.,
$-\mathcal{A}^{*}(-y)\in\mathcal{K}\setminus\\{0\\}$ and $b^{\top}(-y)=0$.

###### Proof.

Since $(z,\omega,\kappa)\in\mbox{range}{\mathcal{A}(\theta)}^{*}$, there
exists $(y,\gamma)\in\mathbb{R}^{m+1}$ such that

$\begin{pmatrix}z\\\ \omega\\\
\kappa\end{pmatrix}=\begin{pmatrix}\mathcal{A}^{*}&c\\\ -b^{\top}&-\theta\\\
\bm{0}^{\top}&1\end{pmatrix}\begin{pmatrix}y\\\ \gamma\end{pmatrix}.$

From this equation, we can easily see that $\kappa=\gamma$, and hence
$z=\mathcal{A}^{*}y+\kappa c$ and $\omega=-b^{\top}y-\kappa\theta$ hold for
$(y,\kappa)$.
(1): If $\kappa>0$, then
$-\mathcal{A}^{*}(\frac{-1}{\kappa}y)+c\in\mathcal{K}$ and
$b^{\top}(\frac{-1}{\kappa}y)-\theta\geq 0$ hold since
$(z,\omega,\kappa)\in\bar{\mathcal{K}}\setminus\\{0\\}$.
(2) & (3): If $\kappa=0$, then $z=-\mathcal{A}^{*}(-y)\in\mathcal{K}$ and
$\omega=b^{\top}(-y)\geq 0$ hold for $y$ since
$(z,\omega,\kappa)\in\bar{\mathcal{K}}\setminus\\{0\\}$.

If $\omega>0$, $y$ satisfies $-\mathcal{A}^{*}(-y)\in\mathcal{K}$ and
$b^{\top}(-y)>0$ and hence, $-y$ is an improving ray of (D).

If $\omega=0$, we can easily see that $b^{\top}(-y)=0$. In addition,
$-\mathcal{A}^{*}(-y)\in\mathcal{K}\setminus\\{0\\}$ holds since
$(z,\omega,\kappa)\neq(0,0,0)$. Thus, $-y$ is a reducing direction for (P). ∎

Proposition 3.4 ensures that feasible solutions of
$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta)}^{*},\bar{\mathcal{K}}\setminus\\{0\\})$
give feasible solutions to (D) if (P) is strongly feasible and (D) is
feasible. By adding one more assumption, we can guarantee that the feasible
solution of
$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta)}^{*},\bar{\mathcal{K}}\setminus\\{0\\})$
gives the optimal solution of (D).

######  Corollary 3.5.

Suppose that (P) is strongly feasible, (D) is feasible, and
$\theta=\theta_{p}$. Then,
$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta)}^{*},\bar{\mathcal{K}}\setminus\\{0\\})$
is feasible, i.e., there exists $y\in\mathbb{R}^{m}$ such that
$z=\mathcal{A}^{*}y+\kappa c$ and $\omega=-b^{\top}y-\kappa\theta$ for any
feasible solution $(z,\omega,\kappa)$ of
$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta)}^{*},\bar{\mathcal{K}}\setminus\\{0\\})$.
Furthermore, $(-\frac{1}{\kappa}y,\frac{1}{\kappa}z)$ is an optimal solution
for (D).

###### Proof.

Since $\theta=\theta_{p}$,
$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta)}^{*},\bar{\mathcal{K}}\setminus\\{0\\})$
is feasible by Corollary 3.3. In addition, for any feasible solution
$(z,\omega,\kappa)$ for
$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta)}^{*},\bar{\mathcal{K}}\setminus\\{0\\})$,
there exists $y\in\mathbb{R}^{m}$ satisfying $z=\mathcal{A}^{*}y+\kappa c$ and
$\omega=-b^{\top}y-\kappa\theta$ by Proposition 3.4.

Noting that (P) is strongly feasible and (D) is feasible, we can see that
$\kappa>0$ holds for any feasible solution $(z,\omega,\kappa)$ for
$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta)}^{*},\bar{\mathcal{K}}\setminus\\{0\\})$
from Proposition 3.4. Thus,
$(-\frac{1}{\kappa}y,c-\mathcal{A}^{*}(\frac{-1}{\kappa}y))$ is a feasible
solution for (D) and $b^{\top}(\frac{-1}{\kappa}y)\geq\theta_{p}$ holds.

Since (P) is strongly feasible and (D) is feasible, by Proposition 2.2, (D)
has an optimal solution and $\theta_{p}=\theta_{d}$ and hence,
$(-\frac{1}{\kappa}y,c-\mathcal{A}^{*}(\frac{-1}{\kappa}y))$ is an optimal
solution for (D). ∎

####  3.1.2 Dual model

Let $\theta\in\mathbb{R}$. We define the linear operator $\mathcal{A}(\theta)$
and the symmetric cone $\bar{\mathcal{K}}$ in the same way as defined in
Section 3.1.1. In this section, we consider the feasibility problem
$\mbox{FP}(\mbox{range}\mathcal{A}(\theta)^{*},\mbox{int }\bar{\mathcal{K}})$.

######  Proposition 3.6.

Suppose that $\mbox{FP}(\mbox{range}\mathcal{A}(\theta)^{*},\mbox{int
}\bar{\mathcal{K}})$ is feasible and $(z,\omega,\kappa)$ is a feasible
solution of it. Then, there exists $y\in\mathbb{R}^{m}$ such that
$z=\mathcal{A}^{*}y+\kappa c$ and $\omega=-b^{\top}y-\kappa\theta$. In
addition, $(-\frac{1}{\kappa}y,c-\mathcal{A}^{*}(\frac{-1}{\kappa}y))$ is an
interior feasible solution to (D), and $b^{\top}(\frac{-1}{\kappa}y)>\theta$
holds.

###### Proof.

Since this proposition can be proved in the same way as the proof of
Proposition 3.1, we omit the proof. ∎

Similar to Proposition 3.1, Proposition 3.6 implies that we can construct the
interior feasible solution of (D) whose objective value is greater than
$\theta$ using the feasible solutions of
$\mbox{FP}(\mbox{range}\mathcal{A}(\theta)^{*},\mbox{int }\bar{\mathcal{K}})$
The following proposition gives us a necessary and sufficient condition for
$\mbox{FP}(\mbox{range}\mathcal{A}(\theta)^{*},\mbox{int }\bar{\mathcal{K}})$
to be feasible.

######  Proposition 3.7.

$\mbox{FP}(\mbox{range}\mathcal{A}(\theta)^{*},\mbox{int }\bar{\mathcal{K}})$
is feasible if and only if (D) is strongly feasible and $\theta_{d}>\theta$.

###### Proof.

Since this proposition can be proved in the same way as the proof of
Proposition 3.2, we omit the proof. ∎

Combining Proposition 3.7 with Proposition 2.10, we have a necessary and
sufficient condition for alternative problem
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\bar{\mathcal{K}}\setminus\\{0\\})$
of $\mbox{FP}(\mbox{range}\mathcal{A}(\theta)^{*},\mbox{int
}\bar{\mathcal{K}})$ to be feasible.

######  Corollary 3.8.

$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\bar{\mathcal{K}}\setminus\\{0\\})$
is feasible if and only if (D) is not strongly feasible or
$\theta_{d}\leq\theta$.

While feasible solutions for
$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta)}^{*},\mbox{int
}\bar{\mathcal{K}})$ give us interior feasible solutions to (D) whose
objective value is greater than $\theta$, feasible solutions for
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\bar{\mathcal{K}}\setminus\\{0\\})$
give us information about the feasibility of (D) or a feasible solution to (P)
whose objective value is smaller than or equal to $\theta$.

######  Proposition 3.9.

Suppose that
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\bar{\mathcal{K}}\setminus\\{0\\})$
is feasible and $(x,\tau,\rho)$ is a feasible solution of
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\bar{\mathcal{K}}\setminus\\{0\\})$.
Then, for $(x,\tau,\rho)$, one of the following three cases holds:

  1. 1.

$\tau>0$ meaning that $\frac{1}{\tau}x$ is a feasible solution to (P) and its
objective value is smaller than or equal to $\theta$,

  2. 2.

$\tau=0$ and $\rho>0$, meaning that $x$ is an improving ray of (P), i.e.,
$x\in\mathcal{K}$, $\mathcal{A}x=0$ and $\langle c,x\rangle<0$, or

  3. 3.

$\tau=\rho=0$ meaning that $x$ is a reucing direction for (D), i.e.,
$x\in\mathcal{K}\setminus\\{0\\}$, $\mathcal{A}x=0$ and $\langle
c,x\rangle=0$.

###### Proof.

Since this proposition can be proved in the same way as the proof of
Proposition 3.4, we omit the proof. ∎

Proposition 3.9 ensures that feasible solutions for
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\bar{\mathcal{K}}\setminus\\{0\\})$
provide feasible solutions for (P) if (D) is strongly feasible and (P) is
feasible. Similar to Corollary 3.5, by adding one more assumption, we can
guarantee that the feasible solution of
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\bar{\mathcal{K}}\setminus\\{0\\})$
gives the optimal solution of (P).

######  Corollary 3.10.

Suppose that (D) is strongly feasible, (P) is feasible, and
$\theta=\theta_{d}$. Then,
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\bar{\mathcal{K}}\setminus\\{0\\})$
is feasible. In addition, for any feasible solution $(x,\tau,\rho)$ for
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\bar{\mathcal{K}}\setminus\\{0\\})$,
$\frac{1}{\tau}x$ is an optimal solution for (P).

###### Proof.

Since this proposition can be proved in the same way as the proof of Corollary
3.5, we omit the proof. ∎

###  3.2 Basic concept of the proposed algorithm

To briefly illustrate the concept of the proposed method, suppose that (P) is
strongly feasible. Then, the feasibility of
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\mbox{int }\bar{\mathcal{K}})$ only
depends on the value of $\theta$ by Proposition 3.2. If $\theta>\theta_{p}$,
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\mbox{int }\bar{\mathcal{K}})$ is
feasible and its solution gives an interior feasible solution $x$ to (P) such
that $\langle c,x\rangle<\theta$ by Proposition 3.1. If
$\theta\leq\theta_{p}$, $\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\mbox{int
}\bar{\mathcal{K}})$ is infeasible and the infeasibility certificates for
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\mbox{int }\bar{\mathcal{K}})$, i.e.,
feasible solutions for
$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta)}^{*},\bar{\mathcal{K}}\setminus\\{0\\})$,
give the dual feasible solution $y$ such that $b^{\top}y\geq\theta$ by
Proposition 3.4. Therefore, the closer the value of $\theta$ is to the optimal
value of (P) and (D), the more accurate the approximate optimal solution of
(P) and (D) obtained by solving
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\mbox{int }\bar{\mathcal{K}})$. In
addition, we can know whether $\theta>\theta_{p}$ or $\theta\leq\theta_{p}$ by
solving $\mbox{FP}(\mbox{ker}\mathcal{A}(\theta),\mbox{int
}\bar{\mathcal{K}})$.

We present our algorithm for finding an approximate optimal interior feasible
solution to (P).

Algorithm 2 Projection and rescaling algorithm with the primal model

1: Input: $\mathcal{A}$, $b$, $c$, $\mathcal{K}$, and $\theta_{acc}>0$.

2: Output: A feasible solution to (P) or a vector that determines the
feasibility status of (P).

3: initialization: $k\leftarrow 0,\ LB\leftarrow-\infty,UB\leftarrow\infty$,
$\bar{\mathcal{K}}\leftarrow\mathcal{K}\times\mathbb{R}^{2}_{+}$

4: Choose $\theta^{k}\in(LB,UB)$ and construct
$\mathcal{A}(\theta^{k})=\begin{pmatrix}\mathcal{A}&-b&\bm{0}\\\
c^{\top}&-\theta^{k}&1\end{pmatrix}.$

5: Let $\varepsilon$ be a sufficiently small positive value and $\xi$ be a
constant such that $0<\xi<1$.

6: while $UB-LB>\theta_{acc}$ do

7: Call the projection and rescaling algorithm of [12] with
$\mbox{ker}\mathcal{A}(\theta^{k})$, $\bar{\mathcal{K}}$, $\varepsilon$, and
$\xi$.

8: if a solution $(x,\tau,\rho)$ to
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta^{k}),\mbox{int }\bar{\mathcal{K}})$ is
obtained then

9: $x_{tmp}\leftarrow\frac{1}{\tau}x$, $UB\leftarrow\theta^{k}$

10: else if a solution $(z,\omega,\kappa)$ to
$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta^{k})}^{*},\bar{\mathcal{K}}\setminus\\{0\\})$
is obtained then

11: Compute $(y,\gamma)\in\mathbb{R}^{m+1}$ such that
$z=\mathcal{A}^{*}y+\gamma c$, $\omega=-b^{\top}y-\gamma\theta$ and
$\kappa=\gamma$.

12: if $-y$ is an improving ray of (D) or a reducing direction for (P) then

13: stop Algorithm 2 and return $(-y,\mathcal{A}^{*}y)$

14: end if

15: $y_{tmp}\leftarrow\frac{-1}{\kappa}y$, $z_{tmp}\leftarrow
c-\mathcal{A}^{*}y_{tmp}$, $LB\leftarrow\theta^{k}$

16: else

17: $LB\leftarrow\theta^{k}$

18: end if

19: Choose $\theta^{k+1}\in(LB,UB)$, $k\leftarrow k+1$

20: end while

21: return $x_{tmp}$

Algorithm 2 works as follows. First, Algorithm 2 chooses the input value
$\theta^{k}\in(LB,UB)$ and execute a projection and rescaling algorithm with
the corresponding problem
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta^{k}),\mbox{int }\bar{\mathcal{K}})$.
Next, Algorithm 2 performs the operations according to the returned result, as
follows:

  1. 1.

If a solution to (P) is obtained, the current primal solution $x_{tmp}$ and
$UB$ are updated.

  2. 2.

If a solution to (D) is obtained, the current dual solution $y_{tmp}$ and $LB$
are updated.

  3. 3.

If a reducing direction for (P) or an improving ray of (D) is obtained, then
Algorithm 2 terminates.

  4. 4.

If a projection and rescaling algorithm determines that the minimum eigenvalue
of any feasible solution of the input problem, i.e.,
$\mbox{FP}_{S_{\infty}}(\mbox{ker }\mathcal{A}(\theta^{k}),\mbox{int
}\bar{\mathcal{K}})$, is less than $\varepsilon$, then $LB$ is updated.

Algorithm 2 repeats the above operations until $UB-LB>\theta_{acc}$ holds.
$UB$ and $LB$ play the role of upper and lower bounds on $\theta_{p}$,
respectively. Thus, if $\theta_{acc}$ is sufficiently small, the output
$x_{tmp}$ from Algorithm 2 can be seen as an approximate optimal solution to
(P).

Here, we note that Algorithm 2 updates $LB$ on line 17. Suppose that Algorithm
2 reaches line 17 at the $k$-th iteration. In this case, it can be inferred
that the projection and rescaling algorithm terminates with the same
termination criteria with $\mbox{FP}_{S_{\infty}}(\mbox{ker
}\mathcal{A}(\theta),\mbox{int }\bar{\mathcal{K}})$ for any $\theta$ such that
$\theta<\theta^{k}$, unless $\theta^{k}$ is not too large compared to the
maximum value of the objective function in (P). Noting that Algorithm 2 is
used as a post-processing step, it is reasonable to update $LB$ as on line 17.
We also note that Algorithm 2 can find an approximate optimal solution to (D).
Suppose that Algorithm 2 reaches line 15 at the $k$-th iteration and
$\theta^{k}$ is very close to the optimal value of (D). Then,
$(y_{tmp},z_{tmp})$ satisfies $b^{\top}y_{tmp}\geq\theta^{k}$ and
$c-\mathcal{A}^{*}y_{tmp}=z_{tmp}\in\mathcal{K}$ by Proposition 3.4. However,
even if $\mbox{FP}(\mbox{ker}\mathcal{A}(\theta^{k}),\mbox{int
}\bar{\mathcal{K}})$ is infeasible, the projection and rescaling algorithm
does not necessarily return a solution to its alternative problem. Thus,
Algorithm 2 is just a method to find the approximate optimal interior feasible
solution of (P).

From the contents of Section 3.1.2, the algorithm for finding an approximate
optimal interior feasible solution to (D) can be considered similarly. (See
Algorithm 3.)

Algorithm 3 Projection and rescaling algorithm with the dual model

1: Input: $\mathcal{A}$, $b$, $c$, $\mathcal{K}$, and $\theta_{acc}>0$.

2: Output: A feasible solution to (D) or a vector that determines the
feasibility status of (D).

3: Same as lines 3-5 of Algorithm 2

4: while $UB-LB>\theta_{acc}$ do

5: Call the projection and rescaling algorithm of [12] with
$\mbox{range}{\mathcal{A}(\theta^{k})}^{*}$, $\bar{\mathcal{K}}$,
$\varepsilon$, and $\xi$.

6: if a solution $(z,\omega,\kappa)$ to
$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta^{k})}^{*},\mbox{int
}\bar{\mathcal{K}})$ is obtained then

7: Compute $(y,\gamma)\in\mathbb{R}^{m+1}$ such that
$z=\mathcal{A}^{*}y+\gamma c$, $\omega=-b^{\top}y-\gamma\theta$ and
$\kappa=\gamma$.

8: $y_{tmp}\leftarrow\frac{-1}{\kappa}y$, $z_{tmp}\leftarrow
c-\mathcal{A}^{*}y_{tmp}$, $LB\leftarrow\theta^{k}$

9: else if a solution $(x,\tau,\rho)$ to $\mbox{FP}(\mbox{ker
}\mathcal{A}(\theta^{k}),\bar{\mathcal{K}}\setminus\\{0\\})$ is obtained then

10: if $x$ is an improving ray of (P) or a reducing direction for (D) then

11: stop Algorithm 3 and return $x$

12: end if

13: $x_{tmp}\leftarrow\frac{1}{\tau}x$, $UB\leftarrow\theta^{k}$

14: else

15: $UB\leftarrow\theta^{k}$

16: end if

17: Choose $\theta^{k+1}\in(LB,UB)$, $k\leftarrow k+1$

18: end while

19: return $y_{tmp}$

###  3.3 Implementation strategies

Algorithms 2 and 3 can operate more efficiently by fully utilizing the
information obtained at each iteration. In this section, we modify Algorithm 2
to make it more practical. These modifications can also be employed in
Algorithm 3.

####  3.3.1 Updating LB using a dual feasible solution

Algorithm 2 keeps $LB$ as the lower bound on the $\theta_{p}$. If a dual
feasible solution $(y_{tmp},z_{tmp})$ is obtained at the $k$-th iteration of
Algorithm 2, we update $LB$ with $\theta^{k}$. Here, we note that
$b^{\top}y_{tmp}\geq\theta^{k}$ by Proposition 3.4 and all dual feasible
solutions give the lower bound on the $\theta_{p}$. Therefore, $LB$ can be
updated with $b^{\top}y_{tmp}$ instead of $\theta^{k}$, which reduces the
number of iterations of Algorithm 2. However, Algorithm 2 obtains approximate
feasible solutions in practice. That is, $z_{tmp}\notin\mathcal{K}$ or
$b^{\top}y_{tmp}<\theta^{k}$ might hold. Considering that $(y_{tmp},z_{tmp})$
might be an approximate dual feasible solution, we update $LB$ as follows:

$LB:=\begin{cases}\max\\{\theta^{k},b^{\top}y_{tmp}\\}&z_{tmp}\in\mathcal{K}\\\
\theta^{k}&\mbox{otherwise}\end{cases}.$ (1)

####  3.3.2 Keeping and updating a current dual feasible solution

The modification in this section has the same spirit as Section 3.3.1. That
is, we update $LB$ using dual feasible solutions. In the previous section, we
propose to update $LB$ as in (1) only if $(y_{tmp},z_{tmp})$ is a feasible
solution for (D). However, even if $(y_{tmp},z_{tmp})$ is an approximate
feasible solution, it can be used to update $LB$ as long as we have a dual
feasible solution $(\bar{y},\bar{z})$. By considering the linear combination
of $y_{tmp}$ and $\bar{y}$, $(y_{new},z_{new})$ might be obtained such that
$z_{new}=c-\mathcal{A}^{*}y_{new}\in\mathcal{K}$ and
$b^{\top}y_{new}\geq\max\\{\theta^{k},b^{\top}\bar{y}\\}$, which leads to
reducing the number of iterations of Algorithm 2. In addition, this
modification allows Algorithm 2 to return the current dual feasible solution
$(\bar{y},\bar{z})$.

Therefore, the following operations are added to Algorithm 2.

  1. 1.

Initialize $\bar{y}$ as $\bar{y}\leftarrow\emptyset$. ( If a dual feasible
solution $y$ is known in advance, initialize $\bar{y}$ and $LB$ as
$\bar{y}\leftarrow y$ and $LB\leftarrow b^{\top}y$, respectively.)

  2. 2.

Suppose that $(y_{tmp},z_{tmp})$ is obtained at the $k$-th iteration. Then, we
perform the following operations.

     * •

If $z_{tmp}\in\mathcal{K}$ and $\bar{y}$ is empty, then $\bar{y}=y_{tmp}$.

     * •

If $\bar{y}$ is not empty, compute $y_{new}$ such that $b^{\top}y_{new}\geq
b^{\top}\bar{y}$ and $c-\mathcal{A}^{*}y_{new}\in\mathcal{K}$, and then update
$LB$ and $\bar{y}$ as $LB=\max\\{\theta^{k},b^{\top}y_{new}\\}$ and
$\bar{y}=y_{new}$, respectively.

  3. 3.

Return $\bar{y}$ on line 21 if $\bar{y}$ is not empty.

In Algorithm 2, $\bar{y}$ plays the role of the current dual feasible
solution. If we have a dual feasible solution before running Algorithm 2, it
can be used to initialize $\bar{y}$. Since Algorithm 2 is used as a post-
processing step for the interior point method, $\bar{y}$ can be initialized
with the output of the interior point method. The method for computing
$y_{new}$ is described in the Appendix A.

####  3.3.3 Use of scaling information from previous iterations

Algorithm 2 calls the projection and rescaling algorithm of [12] iteratively.
It is natural to consider ways to reduce the computational time of the
projection and rescaling algorithm at the $k+1$-th iteration using the
information obtained by the $k$-th iteration. Recall that the projection and
rescaling algorithm of [12] solves the feasibility problem
$\mbox{FP}_{S_{\infty}}(\mathcal{L},\mbox{int }\mathcal{K})$ by repeating two
steps: (i). find a cut for $\mbox{FP}_{S_{\infty}}(\mathcal{L},\mbox{int
}\mathcal{K})$, (ii). scale the problem to an isomorphic problem equivalent to
$\mbox{FP}_{S_{\infty}}(\mathcal{L},\mbox{int }\mathcal{K})$ such that the
region narrowed by the cut is efficiently explored. Therefore, in Algorithm 2,
if the cuts obtained by a projection and rescaling algorithm at the $i(<k)$-th
iteration hold for any feasible solution of the feasibility problem considered
at the $k$-th iteration, such cuts can be used to reduce the execution time of
the projection and rescaling algorithm at the $k$-th iteration. Proposition
3.11 provides the sufficient condition for a cut to
$\mbox{FP}_{S_{\infty}}(\mbox{ker}\mathcal{A}(\bar{\theta}),\mbox{int
}\bar{\mathcal{K}})$ to be valid for
$\mbox{FP}_{S_{\infty}}(\mbox{ker}\mathcal{A}(\theta),\mbox{int
}\bar{\mathcal{K}})$ for two real values $\bar{\theta}\in\mathbb{R}$ and
$\theta\in\mathbb{R}$.

######  Proposition 3.11.

Suppose that (P) is strongly feasible and that $\bar{\theta}\in\mathbb{R}$
satisfies $\bar{\theta}>\theta_{p}$ and $\bar{\theta}-\theta_{p}\leq 1$. Then,
$\mbox{FP}_{S_{\infty}}(\mbox{ker}\mathcal{A}(\theta),\mbox{int
}\bar{\mathcal{K}})$ is feasible for any $\theta\in\mathbb{R}$ such that
$\theta_{p}<\theta<\bar{\theta}$. Furthermore, if $\langle
v,\bar{x}\rangle\leq\xi_{x}$, $\bar{\tau}\leq\xi_{\tau}$ and
$\bar{\rho}\leq\xi_{\rho}$ hold for any feasible solution
$(\bar{x},\bar{\tau},\bar{\rho})$ of
$\mbox{FP}_{S_{\infty}}(\mbox{ker}\mathcal{A}(\bar{\theta}),\mbox{int
}\bar{\mathcal{K}})$ and for some $\xi_{x}<1,\xi_{\tau}<1,\xi_{\rho}<1$ and
$v\in\mathcal{K}$, then, for any $\theta\in\mathbb{R}$ such that
$\theta_{p}<\theta<\bar{\theta}$ and for any feasible solution $(x,\tau,\rho)$
of $\mbox{FP}_{S_{\infty}}(\mbox{ker}\mathcal{A}(\theta),\mbox{int
}\bar{\mathcal{K}})$, $\langle v,x\rangle\leq\xi_{x}$, $\tau\leq\xi_{\tau}$
and $\rho\leq\xi_{\rho}$ hold.

###### Proof.

For any $\theta\in\mathbb{R}$ such that $\theta_{p}<\theta<\bar{\theta}$,
$\mbox{FP}_{S_{\infty}}(\mbox{ker}\mathcal{A}(\theta),\mbox{int
}\bar{\mathcal{K}})$ is feasible from Proposition 3.2. Let $(x,\tau,\rho)$ be
a feasible solution for
$\mbox{FP}_{S_{\infty}}(\mbox{ker}\mathcal{A}(\theta),\mbox{int
}\bar{\mathcal{K}})$. Then, we have

$\langle c,x\rangle-\tau\theta+\rho=0.$ (2)

Noting that $\tau>0$, we have $\langle
c,\frac{1}{\tau}x\rangle+\frac{1}{\tau}\rho=\theta$. Since $\rho>0$ and
$\frac{1}{\tau}x$ is an interior feasible solution for (P),
$\theta_{p}<\langle c,\frac{1}{\tau}x\rangle<\theta<\bar{\theta}$. Thus, we
find

$\bar{\theta}-\langle c,\frac{1}{\tau}x\rangle<\bar{\theta}-\theta_{p}.$ (3)

Let $\alpha=\bar{\theta}-\theta$. By substituting $\theta=\bar{\theta}-\alpha$
into (2), we have

$\langle c,x\rangle-\tau\bar{\theta}+\tau\alpha+\rho=0.$ (4)

Since $\bar{\theta}-\theta_{p}\leq 1$ holds from the assumption, we find

$\displaystyle\tau\alpha+\rho$ $\displaystyle=\tau\bar{\theta}-\langle
c,x\rangle\ \ \ \mbox{(by (\ref{proof: eq3}))}$
$\displaystyle<\tau(\bar{\theta}-\theta_{p})\ \ \ \mbox{(by (\ref{proof:
eq2}))}$ $\displaystyle\leq\tau\leq 1.$

Therefore, $\tau\alpha+\rho<1$ holds, and hence $(x,\tau,\tau\alpha+\rho)$ is
a feasible solution for
$\mbox{FP}_{S_{\infty}}(\mbox{ker}\mathcal{A}(\bar{\theta}),\mbox{int
}\bar{\mathcal{K}})$. From the assumption, $\langle v,x\rangle\leq\xi_{x}$,
$\tau\leq\xi_{\tau}$ and $\tau\alpha+\rho\leq\xi_{\rho}$ hold, which completes
the proof. ∎

Based on Proposition 3.11, we add the following operations to Algorithm 2.

  1. 1.

Initialize $\bar{v}$ as $\bar{v}\leftarrow(e,1,1)$.

  2. 2.

Call the projection and rescaling algorithm with
$Q_{\bar{v}}\left(\mbox{ker}\mathcal{A}(\theta^{k})\right)$ and
$\bar{\mathcal{K}}$ to solve
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta^{k}),\mbox{int }\bar{\mathcal{K}})$.

     * •

Suppose that the projection and rescaling algorithm returns the solution of
the input problem by finding the solution of the scaled problem
$\mbox{FP}_{S_{\infty}}(Q_{v}\left(\mbox{ker}\mathcal{A}(\theta^{k})\right),\mbox{int
}\bar{\mathcal{K}})$, where $v=(v_{1},v_{2},v_{3})\in\mbox{int
}\bar{\mathcal{K}}$.

       * –

After updating $UB$ as $UB\leftarrow\theta^{k}$, if $UB-LB\leq 1$ holds, then
preserve the scaling information as $\bar{v}\leftarrow v$.

####  3.3.4 Use of approximate optimal solutions

Let us consider the indicator $\delta_{\infty}(\mathcal{L}\cap\mbox{int
}\mathcal{K}):=\max_{x}\left\\{{\rm det}(x)\mid x\in\mathcal{L}\cap\mbox{int
}\mathcal{K},\|x\|_{\infty}=1\right\\}$, which is equivalent to the indicators
used in [20] and Section 6.2 of [12]. If $\mathcal{L}\cap\mbox{int
}\mathcal{K}\neq\emptyset$, then $\delta_{\infty}(\mathcal{L}\cap\mbox{int
}\mathcal{K})\in(0,1]$ holds, and if $e\in\mathcal{L}\cap\mbox{int
}\mathcal{K}$, then $\delta_{\infty}(\mathcal{L}\cap\mbox{int }\mathcal{K})=1$
holds. The larger the value of the indicator
$\delta_{\infty}(\mathcal{L}\cap\mbox{int }\mathcal{K})$, the sooner
projection and rescaling algorithms will find a solution for
$\mbox{FP}(\mathcal{L},\mbox{int }\mathcal{K})$. Therefore, scaling with
$v\in\mbox{int }\bar{\mathcal{K}}$ such that
$\delta_{\infty}(Q_{v}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int
}\bar{\mathcal{K}})\geq\delta_{\infty}(\mbox{ker}\mathcal{A}(\theta)\cap\mbox{int
}\bar{\mathcal{K}})$ holds can reduce the computational time of Algorithm 2.

Let us introduce Proposition 3.12, which gives us the basic idea of obtaining
$v\in\mbox{int }\bar{\mathcal{K}}$ such that
$\delta_{\infty}(Q_{v}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int
}\bar{\mathcal{K}})\geq\delta_{\infty}(\mbox{ker}\mathcal{A}(\theta)\cap\mbox{int
}\bar{\mathcal{K}})$.

######  Proposition 3.12.

Suppose that (P) is strongly feasible and for a given $\theta\in\mathbb{R}$
and an interior feasible solution $x\in\mathbb{E}$ for (P),
$\max\\{\theta_{p},\theta-1\\}<\langle c,x\rangle<\theta$ holds. Define
$v\in\mathbb{E}\times\mathbb{R}^{2}$ as $v:=(x^{-\frac{1}{2}},1,1)$. Then, the
following inequality holds:

$\delta_{\infty}(Q_{v}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int
}\bar{\mathcal{K}})\geq\theta-\langle c,x\rangle.$

###### Proof.

It is obvious that $(x,1,\theta-\langle
c,x\rangle)\in\mbox{ker}\mathcal{A}(\theta)\cap\mbox{int }\bar{\mathcal{K}}$.
Noting that $Q_{v}\left((x,1,\theta-\langle
c,x\rangle)\right)=(e,1,\theta-\langle c,x\rangle)\in
Q_{v}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int }\bar{\mathcal{K}}$ and
$\theta-\langle c,x\rangle<1$, we find
$\delta_{\infty}(Q_{v}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int
}\bar{\mathcal{K}})\geq\det(e)\times 1\times(\theta-\langle
c,x\rangle)=\theta-\langle c,x\rangle$. ∎

The next corollary follows similarly to Proposition 3.12.

######  Corollary 3.13.

Suppose that (P) is strongly feasible and for a given $\theta\in\mathbb{R}$,
$\theta_{p}<\theta-1$ holds. For any interior feasible solution
$x\in\mathbb{E}$ for (P) such that $\langle c,x\rangle=\theta-1$, define
$v\in\mathbb{E}\times\mathbb{R}^{2}$ as $v:=(x^{-\frac{1}{2}},1,1)$. Then,
$\delta_{\infty}(Q_{v}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int
}\bar{\mathcal{K}})=1$ holds.

Proposition 3.12 and Corollary 3.13 raise interest in whether the following
assumption holds or not.

######  Assumption 3.14.

Suppose that (P) is strongly feasible. Let $x\in\mathbb{E}$ and
$s\in\mathbb{E}$ be interior feasible solutions for (P) such that
$\max\\{\theta_{p},\theta-1\\}<\langle c,x\rangle<\langle c,s\rangle<\theta$
for a given $\theta\in\mathbb{R}$. Define
$v_{x}\in\mathbb{E}\times\mathbb{R}^{2}$ and
$v_{s}\in\mathbb{E}\times\mathbb{R}^{2}$ as $v_{x}:=(x^{-\frac{1}{2}},1,1)$
and $v_{s}:=(s^{-\frac{1}{2}},1,1)$, respectively. Then, the following
relation holds.

$\delta_{\infty}(Q_{v_{x}}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int
}\bar{\mathcal{K}})>\delta_{\infty}(Q_{v_{s}}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int
}\bar{\mathcal{K}}).$

Unfortunately, the authors could not prove whether Assumption 3.14 holds. What
prevents us from proving this assumption is revealed by the following
proposition.

######  Proposition 3.15.

Suppose that (P) is strongly feasible and that for a given
$\theta\in\mathbb{R}$ and an interior feasible solution $s\in\mathbb{E}$ for
(P), $\max\\{\theta_{p},\theta-1\\}<\langle c,s\rangle<\theta$ holds. Let
$v_{s}:=(s^{-\frac{1}{2}},1,1)$ and $(s_{1},s_{2},s_{3})$ be the point giving
the maximum value of
$\delta_{\infty}(Q_{v_{s}}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int
}\bar{\mathcal{K}})$. If $s_{1}=e$, then for any interior feasible solution
$x\in\mathbb{E}$ for (P) such that $\max\\{\theta_{p},\theta-1\\}<\langle
c,x\rangle<\langle c,s\rangle$,

$\delta_{\infty}(Q_{v_{x}}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int
}\bar{\mathcal{K}})>\delta_{\infty}(Q_{v_{s}}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int
}\bar{\mathcal{K}})$

holds, where $v_{x}:=(x^{-\frac{1}{2}},1,1)$.

###### Proof.

It is obvious that $(e,1,\theta-\langle c,s\rangle)\in
Q_{v_{s}}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int }\bar{\mathcal{K}}$ and
$\left(Q_{s^{\frac{1}{2}}}(s_{1}),s_{2},s_{3}\right)\in\mbox{ker}\mathcal{A}(\theta)\cap\mbox{int
}\bar{\mathcal{K}}$ holds.

If $s_{1}=e$, then $Q_{s^{\frac{1}{2}}}(s_{1})=s$ holds. In addition, we can
easily see that $s_{2}=1$ and $s_{3}=\theta-\langle c,s\rangle$ since $s$ is a
feasible solution for (P) and
$(s,s_{2},s_{3})\in\mbox{ker}\mathcal{A}(\theta)\cap\mbox{int
}\bar{\mathcal{K}}$. Thus, we have

$\delta_{\infty}(Q_{v_{s}}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int
}\bar{\mathcal{K}})=\det(e)\times 1\times(\theta-\langle
c,s\rangle)=\theta-\langle c,s\rangle.$

By Proposition 3.12, for any interior feasible solution $x$ for (P) such that
$\max\\{\theta_{p},\theta-1\\}<\langle c,x\rangle<\langle c,s\rangle$, we find
$\delta_{\infty}(Q_{v_{x}}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int
}\bar{\mathcal{K}})\geq\theta-\langle c,x\rangle$, and hence we have
$\delta_{\infty}(Q_{v_{x}}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int
}\bar{\mathcal{K}})>\delta_{\infty}(Q_{v_{s}}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int
}\bar{\mathcal{K}})$. ∎

As in the proof of Proposition 3.15, Assumption 3.14 holds in the case
$s_{1}=e$. However, proving whether Assumption 3.14 holds when $s_{1}\neq e$
is difficult. In this case, the specific value or upper bound of
$\delta_{\infty}(Q_{v_{s}}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int
}\bar{\mathcal{K}})$ is challenging to obtain. Even if $(s_{1},s_{2},s_{3})$
were obtained, it is unclear how the existence of $(s_{1},s_{2},s_{3})$ such
that $s_{1}\neq e$ affects the value of
$\delta_{\infty}(Q_{v_{x}}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int
}\bar{\mathcal{K}})$ for any interior feasible solution $x$ satisfying
$\max\\{\theta_{p},\theta-1\\}<\langle c,x\rangle<\langle c,s\rangle$. These
obstacles prevent us from proving Assumption 3.14.

So far, we have focused on scaling with an interior feasible solution
$x\in\mathbb{E}$ such that $\max\\{\theta_{p},\theta-1\\}<\langle
c,x\rangle<\theta$. The critical concern for our algorithm is whether scaling
with such $x$ can reduce the execution time of Algorithm 2. In other words,
our algorithm needs to determine whether the following relation holds for any
interior feasible solution $x$ satisfying
$\max\\{\theta_{p},\theta-1\\}<\langle c,x\rangle<\theta$ and
$v=(x^{-\frac{1}{2}},1,1)$.

$\delta_{\infty}(Q_{v}(\mbox{ker}\mathcal{A}(\theta))\cap\mbox{int
}\bar{\mathcal{K}})\geq\delta_{\infty}(\mbox{ker}\mathcal{A}(\theta)\cap\mbox{int
}\bar{\mathcal{K}}).$

Unfortunately, the above relation does not always hold. The easiest
counterexample is when $e$ is a feasible solution for (P) and $\langle
c,e\rangle=\theta-1$, i.e.,
$\delta_{\infty}(\mbox{ker}\mathcal{A}(\theta)\cap\mbox{int
}\bar{\mathcal{K}})=1$. Even if $e$ is not a feasible solution for (P), for
the same reason that it is challenging to prove Assumption 3.14, it is also
difficult to ascertain whether this relation holds. It isn’t easy to see
whether this relation holds when $x$ is an approximate feasible solution.
However, the authors believe that Assumption 3.14 holds and scaling with an
arbitrary interior feasible solution $x$ such that
$\max\\{\theta_{p},\theta-1\\}<\langle c,x\rangle<\theta$ can reduce the
computational time of Algorithm 2 based on the observation of a simple example
discussed in Appendix B.

The observations in Appendix B imply that even if $x$ is an approximate
solution for (P), we can expect that the scaling with such $x$ reduces the
computational time of Algorithm 2, as long as $\langle
c,x\rangle\simeq\max\\{\theta_{p},\theta-1\\}$ and $e$ is not an approximate
or feasible solution for (P) such that $\langle
c,e\rangle\simeq\max\\{\theta_{p},\theta-1\\}$. Noting that Algorithm 2 will
be used as a post-processing step of interior point methods, there is no
problem in supposing that we have an approximate optimal solution
$x\in\mbox{int }\mathcal{K}$ and $\theta$ such that $\theta\simeq\theta_{p}$,
i.e., $\max\\{\theta_{p},\theta-1\\}=\theta_{p}$, before running Algorithm 2.

Thus, the following operations are added to Algorithm 2.

  1. 1.

Initialize $\bar{v}$ as $\bar{v}\leftarrow(e,1,1)$. ( If an approximate primal
solution $x\in\mbox{int }\mathcal{K}$ is known in advance, initialize
$\bar{v}$ as $\bar{v}\leftarrow(x^{-\frac{1}{2}},1,1)$.)

  2. 2.

Call the projection and rescaling algorithm with
$Q_{\bar{v}}\left(\mbox{ker}\mathcal{A}(\theta^{k})\right)$ and
$\bar{\mathcal{K}}$ to solve
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta^{k}),\mbox{int }\bar{\mathcal{K}})$.

     * •

Suppose that the projection and rescaling algorithm returns the solution of
the input problem, and we obtain the solution $(x,\tau,\rho)$ for
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta^{k}),\mbox{int }\bar{\mathcal{K}})$.

       * –

After updating $UB$ as $UB\leftarrow\theta^{k}$, if $UB-LB>1$ holds, then
update $\bar{v}$ as $\bar{v}\leftarrow((\frac{1}{\tau}x)^{-\frac{1}{2}},1,1)$.

Note that Algorithm 2 with this modification does not update $\bar{v}$ in the
item (2) above if $UB-LB\leq 1$. This is because the technique in Section
3.3.3 updates $\bar{v}$ when $UB-LB\leq 1$. This modification is not
theoretically guaranteed to work. Our numerical experiments in Section 4
confirm whether these techniques work well.

####  3.3.5 Extracting a highly accurate solution

Algorithm 2 will obtain approximate solutions many times before they
terminate. Such solutions can be used to make the accuracy of the outputs from
Algorithm 2 robust.

The following operations are added to Algorithm 2. Since our algorithm is
intended to be used as a post-processing step, it is assumed that an
approximate optimal solution $(x^{0},y^{0},z^{0})$ is known in advance.

  1. 1.

Initialize $Sol_{\rm P}$ and $Sol_{\rm D}$ as $Sol_{\rm P}\leftarrow x^{0}$
and $Sol_{\rm D}\leftarrow y^{0}$, respectively.

  2. 2.

If $x_{tmp}$ is obtained at the $k$-th iteration, then we add $x_{tmp}$ to
$Sol_{\rm P}$.

  3. 3.

If $y_{tmp}$ is obtained at the $k$-th iteration, then we add $y_{tmp}$ to
$Sol_{\rm D}$.

  4. 4.

Choose $y^{*}$ as in (6) and then compute $x^{*}$ as in (5) using $y^{*}$
before line 21

  5. 5.

Return $x^{*}$ instead of $x_{tmp}$.

Since the set $Sol_{\rm P}$ will contain many approximate primal solutions,
Algorithm 2 can extract the highly accurate approximate optimal solution
$x^{*}$ from this set as in

$x^{*}\leftarrow\mathop{\rm arg~{}min}\limits_{x\in Sol_{\rm
P}}f(x,y^{*},c-\mathcal{A}^{*}y^{*}),$ (5)

where

$f(x,y,z)=\frac{\|\mathcal{A}x-b\|_{2}}{1+\displaystyle\max_{i=1,\dots,m}|b_{i}|}+\max\left\\{0,\frac{-\lambda_{\min}(x)}{1+\displaystyle\max_{i=1,\dots,m}|b_{i}|}\right\\}+\frac{|\langle
c,x\rangle-b^{\top}y|}{1+|\langle c,x\rangle|+|b^{\top}y|}+\frac{|\langle
x,z\rangle|}{1+|\langle c,x\rangle|+|b^{\top}y|}.$

Each term of the function $f(x,y,z)$ is based on the DIMACS Error [18], which
is a measure of accuracy as an optimal solution for (P) and (D). Note that a
dual solution $(y,z)$ is required to extract $x^{*}$ as in (5). If the dual
optimal solution were known, $f(x,y,z)$ could be used to accurately evaluate
the accuracy of $x$ as the primal optimal solution, but such cases would be
sporadic. Thus, we choose $y^{*}$ from $Sol_{\rm D}$ and use it as the
approximate optimal dual solution to extract $x^{*}$. Algorithm 2 chooses
$y^{*}$ as in

$y^{*}\leftarrow\mathop{\rm arg~{}max}\limits_{y\in Sol}b^{\top}y,$ (6)

where $Sol:=\left\\{y\in Sol_{\rm
D}:\lambda_{\min}(c-\mathcal{A}^{*}y)\geq\min\\{\lambda_{\min}(c-\mathcal{A}^{*}y^{0}),0\\}\right\\}$.
With the modification of Algorithm 2 proposed in Section 3.3.2, we can use a
current dual solution $(\bar{y},c-\mathcal{A}^{*}\bar{y})$ to extract $x^{*}$
as long as $\bar{y}$ is not empty.

###  3.4 Practical versions of Algorithms 2 and 3

We now describe Algorithms 2 and 3 employing the modifications proposed in the
previous section. Both algorithms are designed to use the approximate optimal
solutions $(x^{0},y^{0},z^{0})$ of (P) and (D). Algorithms 4 and 5 are
practical versions of Algorithms 2 and 3, respectively. These algorithms
terminate when $UB-LB\leq\theta_{acc}$ holds or when a vector that determines
the feasibility status of (P) or (D) is found. We consider that a reducing
direction for (P) is obtained when Algorithm 4 finds a vector $(y,\gamma)$
such that

$|\gamma|\leq\mbox{1e-12},\ |-b^{\top}y|\leq\mbox{1e-12},\
\lambda_{\min}(\mathcal{A}^{*}y)\geq\mbox{-1e-12},\ \mbox{and}\
\|\mathcal{A}^{*}y\|>\mbox{1e-12}.$

In addition, we consider that an improving ray of (D) is obtained when
Algorithm 4 finds a vector $(y,\gamma)$ such that

$|\gamma|\leq\mbox{1e-12},\ -b^{\top}y>\mbox{1e-12},\ \mbox{and}\
\frac{1}{-b^{\top}y}\lambda_{\min}(\mathcal{A}^{*}y)\geq\mbox{-1e-12}.$

Similarly, when Algorithm 5 finds a vector $(x,\tau,\rho)$ satisfying

$|\tau|\leq\mbox{1e-12},\ |\langle c,x\rangle|\leq\mbox{1e-12},\
\lambda_{\min}(x)\geq\mbox{-1e-12},\ \mbox{and}\ \|x\|>\mbox{1e-12},$

we consider that a reducing direction for (D) is obtained, and when Algorithm
5 finds a vector satisfying

$|\tau|\leq\mbox{1e-12},\ \langle c,x\rangle<\mbox{-1e-12},\ \mbox{and}\
\frac{1}{-\langle c,x\rangle}\lambda_{\min}(x)\geq\mbox{-1e-12},$

we consider that an improving ray of (P) is obtained.

Algorithm 4 Practical version of Algorithm 2

1: Input: $\mathcal{A}$, $b$, $c$, $\mathcal{K}$, $(x^{0},y^{0},z^{0})$ and
$\theta_{acc}>0$.

2: Output: A primal solution $x$ (and a dual solution $y$) or a vector that
determines the feasibility status of (P).

3: initialization: $k\leftarrow 0,\
LB\leftarrow-\infty,UB\leftarrow\infty,\bar{\mathcal{K}}\leftarrow\mathcal{K}\times\mathbb{R}^{2}_{+},\bar{y}\leftarrow\emptyset,\bar{v}\leftarrow(e,1,1),Sol_{\rm
P}\leftarrow x^{0},Sol_{\rm D}\leftarrow y^{0}$ ( If
$c-\mathcal{A}^{*}y^{0}\in\mathcal{K}$, $\bar{y}\leftarrow y^{0}$,
$LB\leftarrow b^{\top}\bar{y}$. If $x^{0}\in{\rm int}\mathcal{K}$,
$\bar{x}\leftarrow x^{0}$, $\bar{v}\leftarrow(\bar{x}^{-\frac{1}{2}},1,1)$. )

4: Choose $\theta^{k}\in(LB,UB)$ and construct
$\mathcal{A}(\theta^{k})=\begin{pmatrix}\mathcal{A}&-b&\bm{0}\\\
c^{\top}&-\theta^{k}&1\end{pmatrix}.$

5: Let $\varepsilon$ be a sufficiently small positive value and $\xi$ be a
constant such that $0<\xi<1$.

6: while $UB-LB>\theta_{acc}$ do

7: Call the projection and rescaling algorithm of [12] with
$Q_{\bar{v}}(\mbox{ker}\mathcal{A}(\theta^{k}))$, $\bar{\mathcal{K}}$,
$\varepsilon$, and $\xi$ to solve
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta^{k}),\mbox{int }\bar{\mathcal{K}})$
and then obtain the scaling information $v=(v_{1},v_{2},v_{3})$.

8: if a solution $(x,\tau,\rho)$ to
$\mbox{FP}(\mbox{ker}\mathcal{A}(\theta^{k}),\mbox{int }\bar{\mathcal{K}})$ is
obtained then

9: $x_{tmp}\leftarrow\frac{1}{\tau}x$, $Sol_{\rm P}\leftarrow Sol_{\rm
P}\cup\\{x_{tmp}\\}$, $UB\leftarrow\theta^{k}$

10: if $UB-LB\leq 1$ then

11: $\bar{v}\leftarrow(v_{1},v_{2},v_{3})$

12: else

13: if $x_{tmp}\in{\rm int}\mathcal{K}$ then

14: $\bar{v}\leftarrow(x_{tmp}^{-\frac{1}{2}},1,1)$

15: end if

16: end if

17: else if a solution $(z,\omega,\kappa)$ to
$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta^{k})}^{*},\bar{\mathcal{K}}\setminus\\{0\\})$
is obtained then

18: Compute $(y,\gamma)\in\mathbb{R}^{m+1}$ such that
$z=\mathcal{A}^{*}y+\gamma c$, $\omega=-b^{\top}y-\gamma\theta$ and
$\kappa=\gamma$.

19: if $-y$ is an improving ray of (D) or a reducing direction for (P) then

20: stop Algorithm 4 and return $(-y,\mathcal{A}^{*}y)$

21: else

22: $y_{tmp}\leftarrow\frac{-1}{\kappa}y$, $z_{tmp}\leftarrow
c-\mathcal{A}^{*}y_{tmp}$, $Sol_{\rm D}\leftarrow Sol_{\rm
D}\cup\\{y_{tmp}\\}$, $LB\leftarrow\theta^{k}$

23: if $\bar{y}\neq\emptyset$ then

24: Compute $y_{new}$ such that $b^{\top}y_{new}\geq b^{\top}\bar{y}$ and
$c-\mathcal{A}^{*}y_{new}\in\mathcal{K}$ using $\bar{y}$ and $y_{tmp}$.

25: $\bar{y}\leftarrow y_{new}$,
$LB\leftarrow\max\\{\theta^{k},b^{\top}\bar{y}\\}$

26: else

27: if $z_{tmp}\in\mathcal{K}$ then

28: $\bar{y}\leftarrow y_{tmp}$,
$LB\leftarrow\max\\{\theta^{k},b^{\top}\bar{y}\\}$

29: end if

30: end if

31: end if

32: else

33: $LB\leftarrow\theta^{k}$

34: end if

35: Choose $\theta^{k+1}\in(LB,UB)$, $k\leftarrow k+1$

36: end while

37: Choose $y^{*}$ from $Sol_{\rm D}$ as in (6) // If $\bar{y}=\emptyset$,
$y^{*}\leftarrow\bar{y}$.

38: Choose $x^{*}$ from $Sol_{\rm P}$ as in (5)

39: return $x^{*}$ (and $y^{*}$)

Algorithm 5 Practical version of Algorithm 3

1: Input: $\mathcal{A}$, $b$, $c$, $\mathcal{K}$, $(x^{0},y^{0},z^{0})$ and
$\theta_{acc}>0$.

2: Output: A dual solution $y$ (and a primal solution $x$) or a vector that
determines the feasibility status of (D).

3: initialization: $k\leftarrow 0,\
LB\leftarrow-\infty,UB\leftarrow\infty,\bar{\mathcal{K}}\leftarrow\mathcal{K}\times\mathbb{R}^{2}_{+},\bar{y}\leftarrow\emptyset,\bar{v}\leftarrow(e,1,1),Sol_{\rm
P}\leftarrow x^{0},Sol_{\rm D}\leftarrow y^{0}$ ( If
$c-\mathcal{A}^{*}y^{0}\in\mathcal{K}$, $\bar{y}\leftarrow y^{0}$,
$LB\leftarrow b^{\top}\bar{y}$. If $z^{0}\in{\rm int}\mathcal{K}$,
$\bar{z}\leftarrow z^{0}$, $\bar{v}\leftarrow(\bar{z}^{-\frac{1}{2}},1,1)$. )

4: Choose $\theta^{k}\in(LB,UB)$ and construct
$\mathcal{A}(\theta^{k})=\begin{pmatrix}\mathcal{A}&-b&\bm{0}\\\
c^{\top}&-\theta^{k}&1\end{pmatrix}.$

5: Let $\varepsilon$ be a sufficiently small positive value and $\xi$ be a
constant such that $0<\xi<1$.

6: while $UB-LB>\theta_{acc}$ do

7: Call the projection and rescaling algorithm of [12] with
$Q_{\bar{v}}(\mbox{range}{\mathcal{A}(\theta^{k})}^{*})$, $\bar{\mathcal{K}}$,
$\varepsilon$, and $\xi$ to solve
$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta^{k})}^{*},\mbox{int
}\bar{\mathcal{K}})$ and then obtain the scaling information
$v=(v_{1},v_{2},v_{3})$.

8: if a solution $(z,\omega,\kappa)$ to
$\mbox{FP}(\mbox{range}{\mathcal{A}(\theta^{k})}^{*},\mbox{int
}\bar{\mathcal{K}})$ is obtained then

9: Compute $(y,\gamma)\in\mathbb{R}^{m+1}$ such that
$z=\mathcal{A}^{*}y+\gamma c$, $\omega=-b^{\top}y-\gamma\theta$ and
$\kappa=\gamma$.

10: $y_{tmp}\leftarrow\frac{-1}{\kappa}y$, $z_{tmp}\leftarrow
c-\mathcal{A}^{*}y_{tmp}$, $Sol_{\rm D}\leftarrow Sol_{\rm
D}\cup\\{y_{tmp}\\}$, $LB\leftarrow\theta^{k}$

11: if $\bar{y}\neq\emptyset$ then

12: Compute $y_{new}$ such that $b^{\top}y_{new}\geq b^{\top}\bar{y}$ and
$c-\mathcal{A}^{*}y_{new}\in\mathcal{K}$ using $\bar{y}$ and $y_{tmp}$.

13: $\bar{y}\leftarrow y_{new}$,
$LB\leftarrow\max\\{\theta^{k},b^{\top}\bar{y}\\}$

14: else

15: if $z_{tmp}\in\mathcal{K}$ then

16: $\bar{y}\leftarrow y_{tmp}$,
$LB\leftarrow\max\\{\theta^{k},b^{\top}\bar{y}\\}$

17: end if

18: end if

19: if $UB-LB\leq 1$ then

20: $\bar{v}\leftarrow(v_{1},v_{2},v_{3})$

21: else

22: if $z_{tmp}\in{\rm int}\mathcal{K}$ then

23: $\bar{v}\leftarrow(z_{tmp}^{-\frac{1}{2}},1,1)$

24: end if

25: end if

26: else if a solution $(x,\tau,\rho)$ to $\mbox{FP}(\mbox{ker
}\mathcal{A}(\theta^{k}),\bar{\mathcal{K}}\setminus\\{0\\})$ is obtained then

27: if $x$ is an improving ray of (P) or a reducing direction for (D) then

28: stop Algorithm 5 and return $x$

29: end if

30: $x_{tmp}\leftarrow\frac{1}{\tau}x$, $Sol\leftarrow Sol\cup\\{x_{tmp}\\}$,
$UB\leftarrow\theta^{k}$

31: else

32: $UB\leftarrow\theta^{k}$

33: end if

34: Choose $\theta^{k+1}\in(LB,UB)$, $k\leftarrow k+1$

35: end while

36: Choose $y^{*}$ from $Sol_{\rm D}$ as in (6) // If $\bar{y}=\emptyset$,
$y^{*}\leftarrow\bar{y}$.

37: Choose $x^{*}$ from $Sol_{\rm P}$ as in (5)

38: return $y^{*}$ (and $x^{*}$)

##  4 Numerical results

In this section, we apply our algorithm to some instances to show the
numerical performance of the proposed methods. Figure 1 shows the flow of our
numerical experiments. First, we solve the instances with SDP solvers to
obtain approximate optimal solutions $(x^{0},y^{0},z^{0})$ to (P) and (D).
SDPA [34], SDPT3 [29], and Mosek [1] were used in our experiment. Next, we
call Algorithm 6 with $(x^{0},y^{0},z^{0})$. Algorithm 6 calls Algorithm 5 and
Algorithm 4 in that order if $\lambda_{\min}(z^{0})>$-1e-12. Algorithm 5 is
called first to take advantage of the modification in Section 3.3.2. Since
Algorithm 5 can find an approximate optimal interior feasible solution to (D),
a highly accurate optimal dual solution $y$ will likely be obtained after
running Algorithm 5. Such a vector $y$ can reduce the execution time of
Algorithm 4 via the modification proposed in Section 3.3.2. In our experiment,
SDPA and SPDT3 returned $z^{0}$ such that $\lambda_{\min}(z^{0})>$-1e-14 for
all instances. However, Mosek returned $z^{0}$ such that
$\lambda_{\min}(z^{0})<$-1e-10 for some instances in the experiment of Section
4.2. Algorithm 6 processed these problems in the order of Algorithm 4 and
Algorithm 5. This is because we preferred to use the modifications proposed in
Section 3.3.4 rather than find a highly accurate optimal dual solution before
executing Algorithm 4. (The next section will show the effect of the
modifications proposed in Section 3.3.4.) Even if
$\lambda_{\min}(z^{0})<$-1e-12, we can generate an interior point of the cone
$\mathcal{K}$ by adding a slight perturbation to $z^{0}$. However, the
interior points obtained in this way are not likely to be accurate dual
feasible solutions. On the other hand, Algorithm 4 can take advantage of the
modifications proposed in Section 3.3.4 using $x^{0}$. If Algorithm 4 can find
a dual feasible solution, then Algorithm 5, called after Algorithm 4, can also
use the modifications proposed in Section 3.3.4 from the first iteration. We
note that Mosek did not yield $(x^{0},y^{0},z^{0})$ such that
$\lambda_{\min}(z^{0})<$ -1e-12 and $\lambda_{\min}(x^{0})<$ -1e-12 for all
instances in our experiment. Thus, Algorithm 6 did not reach step 30 in our
experiment. As long as $(x^{0},y^{0},z^{0})$ is an approximate interior
feasible solution, step 30 of Algorithm 6 is not reached.

We set the upper limit for the execution time of Algorithms 4 and 5 to 30
minutes and $\theta_{acc}$ = 1e-12. In addition, we added practical
termination conditions to Algorithms 4 and 5 to account for numerical errors.
(See Appendix C.1.) For detailed settings of our algorithm, see Appendix C.

Figure 1: Experiment flow Algorithm 6 Proposed Algorithm

1: Input: $\mathcal{A}$, $b$, $c$, $\mathcal{K}$, $(x^{0},y^{0},z^{0})$ and
$\theta_{acc}>0$.

2: Output: Solutions to (P) and (D) or a vector that determines the
feasibility status of (P) or (D).

3: if $\lambda_{\min}(z^{0})\geq$ $-1$e-12 then

4: $(x^{*},y^{*},z^{*})\leftarrow(x^{0},y^{0},z^{0})$

5: if $z^{*}\notin\mbox{int}\mathcal{K}$ then

6: $z^{*}\leftarrow z^{*}+(-\lambda_{\min}(z^{*})+$1e-15$)e$

7: end if

8: Call Algorithm 5 with $\mathcal{A}$, $b$, $c$, $\mathcal{K}$,
$(x^{*},y^{*},z^{*})$ and $\theta_{acc}>0$.

9: if a solution $(x,y)$ to (P) and (D) is obtained then

10: $x^{*}\leftarrow x$, $y^{*}\leftarrow y$, $z^{*}\leftarrow
c-\mathcal{A}^{*}y$

11: else if $x$, a reducing direction for (D) or an improving ray of (P), is
obtained then

12: stop Algorithm 6 and return $x$

13: end if

14: Call Algorithm 4 with $\mathcal{A}$, $b$, $c$, $\mathcal{K}$,
$(x^{*},y^{*},z^{*})$ and $\theta_{acc}>0$.

15: if a solution $(x,y)$ to (P) and (D) is obtained then

16: $x^{*}\leftarrow x$, $y*\leftarrow y$, $z^{*}\leftarrow
c-\mathcal{A}^{*}y$

17: else if $(y,z)$, a reducing direction for (P) or an improving ray of (D),
is obtained then

18: stop Algorithm 6 and return $(y,z)$

19: end if

20: else if $\lambda_{\min}(x^{0})\geq-1$e-12 then

21: $(x^{*},y^{*},z^{*})\leftarrow(x^{0},y^{0},z^{0})$

22: if $x^{*}\notin\mbox{int}\mathcal{K}$ then

23: $x^{*}\leftarrow x^{*}+(-\lambda_{\min}(x^{*})+$1e-15$)e$

24: end if

25: Call Algorithm 4 with $\mathcal{A}$, $b$, $c$, $\mathcal{K}$,
$(x^{*},y^{*},z^{*})$ and $\theta_{acc}>0$.

26: Same as lines 15 - 19

27: Call Algorithm 5 with $\mathcal{A}$, $b$, $c$, $\mathcal{K}$,
$(x^{*},y^{*},z^{*})$ and $\theta_{acc}>0$.

28: Same as lines 9 - 13

29: else

30: $(x^{*},y^{*},z^{*})\leftarrow(x^{0},y^{0},z^{0})$

31: Same as lines 8 - 19

32: end if

33: return $(x^{*},y^{*},z^{*})$

We conducted three types of numerical experiments. In Section 4.1, we verified
whether the modification proposed in Section 3.3.4 contributes to reducing the
execution time of Algorithm 6. In Section 4.2, we checked whether Algorithm 6
can obtain approximate optimal solutions with higher accuracy than the SDP
solvers. In Section 4.3, we tested whether Algorithm 6 detects the feasibility
status of SDP more accurately than the SDP solvers. All executions were
performed using MATLAB R2022a on an Intel(R) Core(TM) i9-10980XE CPU @ 3.00GHz
machine with 128GB of RAM.

In Section 4.1 and Section 4.2, we measured the accuracy of the output
solution $(x,y,z)$ to (P) and (D) using the DIMACS errors [18]. Letting $d$ be
the dimension of the Euclidean space $\mathbb{E}$ corresponding to
$\mathcal{K}$, the DIMACS errors consists of six measures as follows:

$\begin{array}[]{ll}{\rm
err_{1}}(x,y,z)=\cfrac{\|\mathcal{A}x-b\|_{2}}{1+\displaystyle\max_{i=1,\dots,m}|b_{i}|},&{\rm
err_{2}}(x,y,z)=\max\left\\{0,\cfrac{-\lambda_{\min}(x)}{1+\displaystyle\max_{i=1,\dots,m}|b_{i}|}\right\\},\\\
{\rm
err_{3}}(x,y,z)=\cfrac{\|c-\mathcal{A}^{*}y-z\|_{J}}{1+\displaystyle\max_{i=1,\dots,d}|c_{i}|},&{\rm
err_{4}}(x,y,z)=\max\left\\{0,\cfrac{-\lambda_{\min}(z)}{1+\displaystyle\max_{i=1,\dots,d}|c_{i}|}\right\\},\\\
{\rm err_{5}}(x,y,z)=\cfrac{\langle c,x\rangle-b^{\top}y}{1+|\langle
c,x\rangle|+|b^{\top}y|},&{\rm err_{6}}(x,y,z)=\cfrac{\langle
x,z\rangle}{1+|\langle c,x\rangle|+|b^{\top}y|}.\end{array}$

Primal feasibility is measured by ${\rm err_{1}}(x,y,z)$ and ${\rm
err_{2}}(x,y,z)$ and dual feasibility is measured by ${\rm err_{3}}(x,y,z)$
and ${\rm err_{4}}(x,y,z)$. Optimality is measured by ${\rm err_{5}}(x,y,z)$
and ${\rm err_{6}}(x,y,z)$. We note that the values of ${\rm err_{5}}(x,y,z)$
and ${\rm err_{6}}(x,y,z)$ might be negative as long as $(x,y,z)$ is an
approximate solution. If ${\rm err_{5}}(x,y,z)=-\delta$ and $\delta>0$, then
${\rm err_{5}}(x,y,z)$ corresponds to a “worse” solution than if ${\rm
err_{5}}(x,y,z)=\delta$.

###  4.1 The effectiveness of the modification in Section 3.3.4

The instances used in this experiment are picked from SDPLIB [2]. The
feasibility status of SDPLIB instances can be inferred based on Freund et
al.’s study [9]. Since this experiment was only to evaluate the modification
in Section 3.3.4, we chose only small size instances where both the primal and
dual problems are strongly feasible, i.e., control1, control2, control3,
truss1, truss3, and truss4. We first solved these instances with Mosek to
obtain an approximate optimal interior feasible solution
$(x^{0},y^{0},z^{0})$. However, Mosek returned dual solutions $(y^{0},z^{0})$
such that $c-\mathcal{A}^{*}y^{0}\notin\mathcal{K}$ and
$z^{0}\notin\mathcal{K}$ to truss1, truss3, and truss4. Thus, we only checked
for six instances how the computation time of Algorithm 4 changed with and
without the modification proposed in Section 3.3.4. In this section, we
referred to Algorithm 4 without the modification proposed in Section 3.3.4 as
the “Naive” method.

Table 2 summarizes the results of our experiments. The “times(s)” column shows
the average CPU time of the corresponding method, and the other columns show
the average value of the DIMACS errors. Table 2 shows that the computational
time of Algorithm 4 is significantly smaller than that of the “Naive” method.
Algorithm 4 scales the problem with an approximate optimal interior feasible
solution before executing the projection and rescaling algorithm. Such scaling
reduced the number of iterations of the projection and rescaling algorithm in
this experiment. Thus, the computational time of Algorithm 4 is much smaller
than that of the “Naive” method.

Table 2: Experimental results testing the effectiveness of the modification proposed in Section 3.3.4 Method | ${\rm err_{1}}$ | ${\rm err_{2}}$ | ${\rm err_{3}}$ | ${\rm err_{4}}$ | ${\rm err_{5}}$ | ${\rm err_{6}}$ | time(s)
---|---|---|---|---|---|---|---
Naive | 7.15e-14 | 0 | 0 | 0 | 3.89e-10 | 3.89e-10 | 13.36
Algorithm 4 | 1.50e-14 | 4.81e-17 | 0 | 0 | 1.81e-13 | 1.81e-13 | 2.66

###  4.2 The effectiveness of our post-processing algorithm

The outline of our numerical experiment in this section is as follows:

  1. 1.

Solve the instances with the SDP solvers using the default settings to obtain
approximate optimal solutions $(x_{\rm def},y_{\rm def},z_{\rm def})$.

  2. 2.

Solve the instances using $(x_{\rm def},y_{\rm def},z_{\rm def})$ and
Algorithm 6 to obtain approximate optimal solutions $(x_{\rm pro},y_{\rm
pro},z_{\rm pro})$.

  3. 3.

Solve the instances with the SDP solvers using the tight settings to obtain
approximate optimal solutions $(x_{\rm tight},y_{\rm tight},z_{\rm tight})$.

  4. 4.

Compare the above three results regarding accuracy and computational time.

In this experiment, we used SDPA, SDPT3, and Mosek. These solvers have
parameters that allow us to tune the feasibility and the optimality
tolerances. For example, SDPA terminates if an approximate optimal solution
$(x,y,z)$ is obtained for parameters $\epsilon_{\rm primal}$, $\epsilon_{\rm
dual}$ and $\epsilon_{\rm gap}$ such that

$\begin{array}[]{lll}\displaystyle\max_{i=1,\dots,m}|(\mathcal{A}x-b)_{i}|\leq\epsilon_{\rm
primal},&\displaystyle\max_{i=1,\dots,d}|(c-\mathcal{A}^{*}y-z)_{i}|\leq\epsilon_{\rm
dual},\ \mbox{and}&\cfrac{|\langle
c,x\rangle-b^{\top}y|}{\max\left\\{\left(|\langle
c,x\rangle|+|b^{\top}y|\right)/2,1\right\\}}\leq\epsilon_{\rm gap}\end{array}$

are satisfied. That is, the parameters $\epsilon_{\rm primal}$, $\epsilon_{\rm
dual}$ and $\epsilon_{\rm gap}$ are the tolerances for primal feasibility,
dual feasibility, and optimality measures, respectively. The definitions of
feasibility and optimality measures are slightly different for each solver,
but the other solvers have parameters that play the same role. Table 3
summarizes the values of the parameters used in our experiment. The “default
tolerances” column shows the default values of the parameters for each solver.
The “tight tolerances” column shows the values of the parameters used as the
tight setting in our experiment.

Table 3: Solver setting in our numerical experiments | default tolerances | tight tolerances
---|---|---
Solver | $\epsilon_{\rm primal}$ | $\epsilon_{\rm dual}$ | $\epsilon_{\rm gap}$ | $\epsilon_{\rm primal}$ | $\epsilon_{\rm dual}$ | $\epsilon_{\rm gap}$
Mosek | 1.0e-08 | 1.0e-12
SDPT3 | 1.0e-08
SDPA | 1.0e-07

The instances of SDPLIB were used in this experiment. Note that the instances
for which the SDP solver using the default setting returned infeasibility
criteria were excluded. In addition, we excluded some instances due to memory
limitations and then conducted numerical experiments with 57 instances in
total. To appropriately observe the results, we classified these instances
into two groups, the well-posed and the ill-posed groups, based on Freund et
al.’s study [9]. The well-posed group includes 34 instances where both the
primal and dual problems are expected to be strongly feasible. On the other
hand, the ill-posed group includes 23 instances where either the primal or
dual problem is not expected to be strongly feasible. According to [9], SDPLIB
does not include the instances whose primal problem is strongly feasible, but
the dual is not. Thus, the ill-posed group includes only instances where the
dual problem is expected to be strongly feasible, but the primal problem is
not.

The numerical results for the well-posed group and ill-posed group are
summarized in Figures 2-4 and 5-7, respectively. Figures 2-7 summarize the
values of DIMACS errors for $(x_{\rm def},y_{\rm def},z_{\rm def})$, $(x_{\rm
tight},y_{\rm tight},z_{\rm tight})$ and $(x_{\rm pro},y_{\rm pro},z_{\rm
pro})$. In this experiment, no output $(x_{\rm pro},y_{\rm pro},z_{\rm pro})$
was obtained such that ${\rm err_{2}}(x_{\rm pro},y_{\rm pro},z_{\rm
pro})>$1e-11 or ${\rm err_{4}}(x_{\rm pro},y_{\rm pro},z_{\rm pro})>0$.
Therefore, we omitted figures summarizing the values of ${\rm err_{2}}$ and
${\rm err_{4}}$. In Figures 2-7, the black solid line, the blue circle and the
inverted orange triangle show the values of the DIMACS errors logarithmized by
the base $10$ for $(x_{\rm def},y_{\rm def},z_{\rm def})$, $(x_{\rm
tight},y_{\rm tight},z_{\rm tight})$ and $(x_{\rm pro},y_{\rm pro},z_{\rm
pro})$, respectively. It is not plotted if the corresponding DIMACS error
value is $0$. We note that the values of $\log_{10}|{\rm err_{5}}(x,y,z)|$ and
$\log_{10}|{\rm err_{6}}(x,y,z)|$ are plotted because the values of ${\rm
err_{5}}(x,y,z)$ and ${\rm err_{6}}(x,y,z)$ can be negative. In addition, we
plot the points where the corresponding DIMACS error values are negative by
filling them in. In each graph of Figures 2-7, the horizontal axis indicates
the instances, sorted so that the corresponding DIMACS error values for
$(x_{\rm def},y_{\rm def},z_{\rm def})$ are in ascending order.

First, let us compare the results for the well-posed group. From Figures 2-4,
we can observe that Algorithm 6 returned an approximate optimal solution
$(x_{\rm pro},y_{\rm pro},z_{\rm pro})$ that was superior to the $(x_{\rm
def},y_{\rm def},z_{\rm def})$ and $(x_{\rm tight},y_{\rm tight},z_{\rm
tight})$ in terms of feasibility and optimality for almost all instances.
Algorithm 6 returns $(x_{\rm pro},y_{\rm pro},z_{\rm pro})$ such that $z_{\rm
pro}=c-\mathcal{A}^{*}y_{\rm pro}$, thus the value of ${\rm err_{3}}(x_{\rm
pro},y_{\rm pro},z_{\rm pro})$ should be $0$ for all instances. However, in
(b) of Figures 2-4, four inverted triangles representing the results for four
instances, theta1, theta2, theta3, and theta4, were plotted. This is probably
due to slight numerical errors in computing $c-\mathcal{A}^{*}y_{\rm pro}$ for
these instances. Note that $(x_{\rm pro},y_{\rm pro},z_{\rm pro})$ satisfies
$\langle c,x_{\rm pro}\rangle–b^{\top}y_{\rm pro}>0$ for almost all instances.

We can also see that the solvers with the tight setting tend to return more
accurate approximate optimal solutions $(x_{\rm tight},y_{\rm tight},z_{\rm
tight})$ than $(x_{\rm def},y_{\rm def},z_{\rm def})$. However, Mosek and
SDPA, with the tight setting, returned strange outputs for some instances.
Mosek, with the default setting, determined the primal and dual problems are
expected to be feasible and returned an approximate optimal solution for all
instances. The error message “rescode = 100006” was returned for only one
instance. This error message means that “the optimizer is terminated due to
slow progress”. On the other hand, when using Mosek with the tight settings,
it returned the error message “rescode = 100006” for 21 instances. Among these
instances, 4 instances were determined that their primal problems were
expected to be infeasible, and 13 instances were not determined their
feasibility. For the aforementioned 4 instances and 13 instances, Mosek, with
the tight setting, provided incorrect infeasibility certificates and
inaccurate dual feasible solutions, respectively. SDPA, with the default
setting, determined the primal and/or dual problems are feasible and returned
an approximate optimal solution for all instances. However, when using SDPA
with tight settings, it was determined that at least one of the primal or dual
problems was expected to be infeasible for 13 instances, and incorrect primal
solutions were returned for 11 instances of them. A possible reason for these
strange outputs is that tightening the tolerances might cause numerical
instability. Figure 4 shows that SDPT3 with the tight setting worked stably
and obtained more accurate approximate optimal solutions $(x_{\rm
tight},y_{\rm tight},z_{\rm tight})$ than $(x_{\rm def},y_{\rm def},z_{\rm
def})$ for almost all instances. We also solved all instances using SDPT3 with
the tolerances $\epsilon_{\rm primal}$, $\epsilon_{\rm dual}$ and
$\epsilon_{\rm gap}$ set to 1e-13 and obtained the same results as Figure 4.

Next, let us compare the results for the ill-posed group. Since Mosek, using
the default setting, returned the reducing direction for hinf12, Figure 5
excludes the plot representing the result for this instance. On the other
hand, Algorithm 6 did not find a reducing direction for (P) in this
experiment. Thus, Algorithm 6 terminated by returning an approximate optimal
solution for all instances. From Figures 5-7, we can observe that Algorithm 6
returned an approximate optimal solution $(x_{\rm pro},y_{\rm pro},z_{\rm
pro})$ that was superior to the $(x_{\rm def},y_{\rm def},z_{\rm def})$ and
$(x_{\rm tight},y_{\rm tight},z_{\rm tight})$ in terms of feasibility for
almost all instances. In addition, we can see that the value of $|{\rm
err_{5}}(x_{\rm pro},y_{\rm pro},z_{\rm pro})|$ or $|{\rm err_{6}}(x_{\rm
pro},y_{\rm pro},z_{\rm pro})|$ was larger than the value of $|{\rm
err_{5}}(x_{\rm def},y_{\rm def},z_{\rm def})|$ or $|{\rm err_{6}}(x_{\rm
def},y_{\rm def},z_{\rm def})|$ for some instances. However, this observation
does not imply that $(x_{\rm pro},y_{\rm pro},z_{\rm pro})$ was inferior to
$(x_{\rm def},y_{\rm def},z_{\rm def})$ and $(x_{\rm tight},y_{\rm
tight},z_{\rm tight})$ in terms of optimality.

Tables 4-6 compare the dual objective value of $(x_{\rm pro},y_{\rm
pro},z_{\rm pro})$ with the primal or dual objective values of $(x_{\rm
def},y_{\rm def},z_{\rm def})$ and $(x_{\rm tight},y_{\rm tight},z_{\rm
tight})$. From Tables 4-6, we can see that

$b^{\top}y_{\rm pro}>\max\\{b^{\top}y_{\rm def},b^{\top}y_{\rm tight}\\}$

holds for almost all instances of the ill-posed group. Noting that
$c-\mathcal{A}^{*}y_{\rm pro}=z_{\rm pro}\in\mathcal{K}$ holds for all
instances, $(y_{\rm pro},z_{\rm pro})$ can be regarded as a more accurate
optimal dual solution than $(y_{\rm def},z_{\rm def})$ and $(y_{\rm
tight},z_{\rm tight})$. Therefore, the reason why the value of $|{\rm
err_{5}}(x_{\rm pro},y_{\rm pro},z_{\rm pro})|$ or $|{\rm err_{6}}(x_{\rm
pro},y_{\rm pro},z_{\rm pro})|$ was sometimes greater than the corresponding
DIMACS error values for $(x_{\rm def},y_{\rm def},z_{\rm def})$ and $(x_{\rm
tight},y_{\rm tight},z_{\rm tight})$ is that Algorithm 6 did not obtain a
sufficient accurate primal solution $x_{\rm pro}$ such that $\langle c,x_{\rm
pro}\rangle\simeq b^{\top}y_{\rm pro}$ for some instances. Here, we note that
the primal problem of all instances in the ill-posed group is not expected to
be strongly feasible, which might prevent the projection and rescaling methods
in Algorithm 6 from working stably and obtaining accurate optimal primal
solutions.

Table 4: Comparison of the objective values of $(x_{\rm pro},y_{\rm pro},z_{\rm pro})$, $(x_{\rm def},y_{\rm def},z_{\rm def})$ and $(x_{\rm tight},y_{\rm tight},z_{\rm tight})$ for the ill-posed group : Mosek instance | $b^{\top}y_{\rm pro}-\langle c,x_{\rm def}\rangle$ | $b^{\top}y_{\rm pro}-b^{\top}y_{\rm def}$ | $b^{\top}y_{\rm pro}-\langle c,x_{\rm tight}\rangle$ | $b^{\top}y_{\rm pro}-b^{\top}y_{\rm tight}$
---|---|---|---|---
gpp100 | 7.38e-06 | 7.37e-06 | 7.38e-06 | 7.37e-06
gpp124-1 | 4.88e-06 | 4.86e-06 | 1.24e-06 | 1.23e-06
gpp124-2 | 5.75e-05 | 5.81e-05 | 5.75e-05 | 5.81e-05
gpp124-3 | 4.76e-05 | 4.94e-05 | 4.76e-05 | 4.94e-05
gpp124-4 | 5.10e-05 | 5.09e-05 | 5.10e-05 | 5.09e-05
hinf1 | 6.65e-07 | 6.01e-07 | 1.30e-06 | 1.18e-06
hinf3 | 3.79e-05 | 3.59e-05 | 3.85e-05 | 3.64e-05
hinf4 | 6.85e-05 | 6.61e-05 | 1.62e-05 | 1.53e-05
hinf5 | 8.33e+00 | 8.33e+00 | 8.33e+00 | 8.33e+00
hinf6 | 1.86e-03 | 1.75e-03 | 7.51e-03 | 7.21e-03
hinf7 | 7.17e+00 | 7.17e+00 | 7.17e+00 | 7.17e+00
hinf8 | 4.63e-03 | 4.46e-03 | 4.63e-03 | 4.46e-03
hinf10 | 7.17e-05 | 4.00e-05 | 3.09e-04 | 2.43e-04
hinf11 | 2.23e-04 | 1.93e-04 | 1.60e-03 | 1.48e-03
hinf13 | 3.40e-04 | 2.59e-04 | 3.40e-04 | 2.59e-04
hinf14 | 5.49e-04 | 5.25e-04 | 5.58e-04 | 5.34e-04
hinf15 | 8.12e-04 | 6.42e-04 | 2.24e-03 | 2.08e-03
qap5 | 1.29e-06 | 1.22e-06 | 1.90e-09 | 1.79e-09
qap6 | 4.38e-04 | 4.36e-04 | 2.68e-03 | 2.61e-03
qap7 | 4.81e-05 | 4.43e-05 | 8.71e-04 | 8.53e-04
qap8 | 5.26e-03 | 5.25e-03 | 4.50e-03 | 4.49e-03
qap9 | 1.14e-03 | 1.13e-03 | 4.89e-03 | 4.84e-03
qap10 | 1.56e-02 | 1.55e-02 | 1.56e-02 | 1.55e-02
Table 5: Comparison of the objective values of $(x_{\rm pro},y_{\rm pro},z_{\rm pro})$, $(x_{\rm def},y_{\rm def},z_{\rm def})$ and $(x_{\rm tight},y_{\rm tight},z_{\rm tight})$ for the ill-posed group : SDPA instance | $b^{\top}y_{\rm pro}-\langle c,x_{\rm def}\rangle$ | $b^{\top}y_{\rm pro}-b^{\top}y_{\rm def}$ | $b^{\top}y_{\rm pro}-\langle c,x_{\rm tight}\rangle$ | $b^{\top}y_{\rm pro}-b^{\top}y_{\rm tight}$
---|---|---|---|---
gpp100 | -1.88e-07 | 2.07e-07 | -1.55e-03 | 3.00e-05
gpp124-1 | -1.29e-07 | 5.94e-08 | -1.04e-03 | 1.28e-05
gpp124-2 | -4.51e-06 | 1.03e-07 | -1.21e-03 | 9.34e-06
gpp124-3 | 5.07e-05 | 6.27e-05 | -2.05e-03 | 1.04e-04
gpp124-4 | -3.14e-05 | 1.88e-06 | -4.54e-03 | 1.03e-04
hinf1 | 1.18e-05 | 1.27e-05 | -2.06e-04 | 1.09e-04
hinf3 | -5.02e-03 | 9.78e-03 | 4.19e-03 | 5.67e-03
hinf4 | 3.46e-04 | 3.46e-04 | 5.96e-04 | 4.22e-04
hinf5 | 2.05e+01 | 2.00e+01 | 2.06e+01 | 2.01e+01
hinf6 | 4.65e-02 | 2.41e-02 | 2.86e-01 | 2.40e-01
hinf7 | -3.08e+00 | 4.77e+00 | -3.08e+00 | 4.77e+00
hinf8 | 1.69e-01 | 1.69e-01 | 1.68e-01 | 1.68e-01
hinf10 | 1.07e-01 | 5.37e-02 | 9.72e+00 | 6.91e+00
hinf11 | 1.25e-01 | 6.30e-02 | 3.11e+00 | 2.58e+00
hinf12 | 2.63e+01 | 3.21e+00 | 4.93e+01 | 4.17e+01
hinf13 | 4.39e+00 | 2.88e+00 | 4.39e+00 | 2.88e+00
hinf14 | -2.90e-03 | 3.14e-03 | -2.03e-02 | 1.84e-02
hinf15 | 3.63e+00 | 2.64e+00 | 3.63e+00 | 2.64e+00
qap5 | -2.24e-02 | 3.36e-04 | -1.60e-03 | 1.08e-04
qap6 | -5.04e-02 | 1.21e-02 | -1.93e-02 | 5.86e-03
qap7 | -1.01e-01 | 2.34e-02 | -4.57e-02 | 1.58e-02
qap8 | -1.53e-01 | 2.11e-02 | -4.95e-01 | 1.62e-01
qap9 | -3.95e-01 | 8.01e-02 | -3.01e-01 | 5.72e-02
qap10 | -7.21e-01 | 8.16e-02 | -6.57e-01 | 7.75e-02
Table 6: Comparison of the objective values of $(x_{\rm pro},y_{\rm pro},z_{\rm pro})$, $(x_{\rm def},y_{\rm def},z_{\rm def})$ and $(x_{\rm tight},y_{\rm tight},z_{\rm tight})$ for the ill-posed group : SDPT3 instance | $b^{\top}y_{\rm pro}-\langle c,x_{\rm def}\rangle$ | $b^{\top}y_{\rm pro}-b^{\top}y_{\rm def}$ | $b^{\top}y_{\rm pro}-\langle c,x_{\rm tight}\rangle$ | $b^{\top}y_{\rm pro}-b^{\top}y_{\rm tight}$
---|---|---|---|---
gpp100 | 2.88e-06 | 1.88e-06 | 3.70e-06 | 1.84e-06
gpp124-1 | 9.15e-07 | 3.71e-07 | 9.20e-07 | 3.71e-07
gpp124-2 | 1.71e-06 | 1.11e-06 | 2.25e-06 | 1.09e-06
gpp124-3 | 4.48e-06 | 2.94e-06 | 5.65e-06 | 2.88e-06
gpp124-4 | 2.72e-05 | 1.20e-05 | 2.78e-05 | 1.20e-05
hinf1 | 1.03e-04 | 5.19e-05 | 1.03e-04 | 5.19e-05
hinf3 | 2.71e-02 | 1.36e-02 | 2.71e-02 | 1.36e-02
hinf4 | 1.92e-03 | 9.68e-04 | 1.92e-03 | 9.68e-04
hinf5 | 4.81e+00 | 4.47e+00 | 4.81e+00 | 4.47e+00
hinf6 | 1.60e-02 | 1.73e-02 | 2.89e-02 | 1.46e-02
hinf7 | 1.42e-02 | 7.48e-03 | 1.42e-02 | 7.48e-03
hinf8 | 4.06e-02 | 2.03e-02 | 4.06e-02 | 2.03e-02
hinf10 | 6.14e-02 | 5.28e-02 | 1.35e-01 | 6.73e-02
hinf11 | 7.37e-02 | 3.68e-02 | 7.37e-02 | 3.68e-02
hinf12 | 2.97e-05 | 1.45e-05 | -8.09e-06 | -1.29e-05
hinf13 | 1.03e-02 | 4.33e-03 | 6.81e-03 | 2.03e-03
hinf14 | 6.37e-05 | 8.15e-05 | 5.90e-05 | 3.54e-05
hinf15 | 1.80e-02 | 2.76e-02 | 1.80e-02 | 2.76e-02
qap5 | -3.08e-07 | 8.53e-08 | -4.55e-09 | 1.58e-09
qap6 | 4.52e-02 | 2.22e-02 | 4.52e-02 | 2.22e-02
qap7 | 3.15e-02 | 1.55e-02 | 3.15e-02 | 1.55e-02
qap8 | 1.13e-01 | 5.54e-02 | 1.13e-01 | 5.54e-02
qap9 | 2.19e-02 | 1.06e-02 | 2.19e-02 | 1.06e-02
qap10 | 5.11e-02 | 1.69e-02 | 5.11e-02 | 1.69e-02

From the results of the solvers for the ill-posed group, we can observe the
same thing as for the well-posed group. When using SDPA and Mosek with the
tight setting, they sometimes returned inaccurate primal solutions and
inaccurate dual solutions, respectively. In addition, SDPT3 with the tight
setting returned more accurate optimal solutions $(x_{\rm tight},y_{\rm
tight},z_{\rm tight})$ than $(x_{\rm def},y_{\rm def},z_{\rm def})$ for almost
all instances. We note that SDPT3 obtained the same results when
$\epsilon_{\rm primal}=\epsilon_{\rm dual}=\epsilon_{\rm gap}=$ 1e-13 as when
$\epsilon_{\rm primal}=\epsilon_{\rm dual}=\epsilon_{\rm gap}=$ 1e-12.

Let us compare the results for the well-posed and ill-posed groups. Comparing
Figures 2-4 and 5-7, we can observe the following:

  * •

Algorithm 6 did not consistently obtain an accurate optimal solution for the
ill-posed group compared to the results for the well-posed group.

  * •

Algorithm 6 returned $(x_{\rm pro},y_{\rm pro},z_{\rm pro})$ such that $|{\rm
err_{5}}(x_{\rm pro},y_{\rm pro},z_{\rm pro})|\simeq|{\rm err_{6}}(x_{\rm
pro},y_{\rm pro},z_{\rm pro})|$ for almost all instances of the well-posed
group, but such a relation did not hold for the ill-posed group.

The first observation is evident in Figures 2-7, (c). To clarify the second
observation, see Table 7. Table 7 summarizes the average value of $||{\rm
err_{5}}(x_{\rm pro},y_{\rm pro},z_{\rm pro})|-|{\rm err_{6}}(x_{\rm
pro},y_{\rm pro},z_{\rm pro})||$ for each group. Table 7 shows that for all
combinations of Algorithm 6 and the solvers, the average value of $||{\rm
err_{5}}(x_{\rm pro},y_{\rm pro},z_{\rm pro})|-|{\rm err_{6}}(x_{\rm
pro},y_{\rm pro},z_{\rm pro})||$ for the ill-posed group is greater than for
the well-posed group. From this table, we can see that the relation $|{\rm
err_{5}}(x_{\rm pro},y_{\rm pro},z_{\rm pro})|\simeq|{\rm err_{6}}(x_{\rm
pro},y_{\rm pro},z_{\rm pro})|$ holds for the well-posed group but not for the
ill-posed group.

Table 7: Comparison of the average values of $\left||{\rm err_{5}}(x_{\rm pro},y_{\rm pro},z_{\rm pro})|-|{\rm err_{6}}(x_{\rm pro},y_{\rm pro},z_{\rm pro})|\right|$ Method | Well-posed group | Ill-posed group
---|---|---
Mosek + Algorithm 6 | 5.34e-11 | 3.45e-03
SDPA + Algorithm 6 | 1.93e-11 | 3.75e-03
SDPT3 + Algorithm 6 | 2.51e-12 | 3.51e-04

There are two reasons for this. The first one is Algorithm 6 did not obtain
accurate optimal primal solutions for some instances of the ill-posed group.
As mentioned above, the projection and rescaling methods might not work stably
due to the feasibility status of the primal problem for the ill-posed group,
which results in the instability of Algorithm 6 in obtaining accurate optimal
primal solutions. Therefore, the relation $|{\rm err_{5}}(x_{\rm pro},y_{\rm
pro},z_{\rm pro})|\simeq|{\rm err_{6}}(x_{\rm pro},y_{\rm pro},z_{\rm pro})|$
did not hold for the ill-posed group compared to the results for the well-
posed group. The second one is the existence of a reducing direction for (P).
Recall that $(f,-\mathcal{A}^{*}f)\in\mathbb{R}^{m}\times\mathbb{E}$ is called
a reducing direction for (P) if $(f,-\mathcal{A}^{*}f)$ satisfies
$b^{\top}f=0$ and $-\mathcal{A}^{*}f\in\mathcal{K}\setminus\\{0\\}$. Thus, if
(D) has an optimal solution and a reducing direction for (P) exists, then the
optimal solution set of (D) is unbounded. Suppose that optimal dual solutions
and reducing directions for (P) exist for all instances of the ill-posed
group. Then, for any optimal dual solution $(y,z)$ and any positive value
$k>0$, there exists a dual optima solution $(y_{opt},z_{opt})$ and a reducing
direction $(f,-\mathcal{A}^{*}f)$ such that $y=y_{opt}+f$ and
$z_{opt}+k\mathcal{A}^{*}f\notin\mathcal{K}$, and we have

$\displaystyle{\rm err_{5}}(x,y,z)$ $\displaystyle=\cfrac{\langle
c,x\rangle-b^{\top}(y_{opt}+f)}{1+|\langle c,x\rangle|+|b^{\top}y|}$
$\displaystyle=\cfrac{\langle c,x\rangle-b^{\top}y_{opt}}{1+|\langle
c,x\rangle|+|b^{\top}y|},\ \ \ \mbox{(since $b^{\top}f=0$)}$

and

$\displaystyle{\rm err_{6}}(x,y,z)$ $\displaystyle=\cfrac{\langle
x,c-\mathcal{A}^{*}y_{opt}-\mathcal{A}^{*}f\rangle}{1+|\langle
c,x\rangle|+|b^{\top}y|}$ $\displaystyle=\cfrac{\langle
x,c\rangle-(\mathcal{A}x)^{\top}y_{opt}+\langle
x,-\mathcal{A}^{*}f\rangle}{1+|\langle c,x\rangle|+|b^{\top}y|}.$

Therefore, even if $\|\mathcal{A}x_{\rm pro}-b\|_{2}$ is sufficiently small,
the value of $\langle x_{pro},-\mathcal{A}^{*}f\rangle$ might not be
negligibly small as long as $y_{\rm pro}\simeq y_{opt}+f$ holds for some
reducing direction $(f,-\mathcal{A}^{*}f)$ such that the value of
$\|-\mathcal{A}^{*}f\|$ is very large.

Table 8 summarizes the average values of
$f(z)=\frac{\|z\|}{1+\displaystyle\max_{i}|c_{i}|+\displaystyle\max_{i,j}|A_{ij}|}$
for $z_{\rm def}$ and $z_{\rm pro}$ by the well-posed and the ill-posed
groups, where $A\in\mathbb{R}^{m\times d}$ is a matrix representation of the
linear operator $\mathcal{A}$. Since Mosek and SDPA with the tight settings
did not work stably in this experiment, the average values of $f(z_{\rm
tight})$ for each solver are omitted. In Table 8, the first, third, and fifth
rows show the average values of $f(z)$ for $z_{\rm def}$ obtained from Mosek,
SDPA, and SDPT3, respectively. The second, fourth, and sixth rows show the
average values of $f(z)$ for $z_{\rm pro}$ obtained from Algorithm 6 using
$(x_{\rm def},y_{\rm def},z_{\rm def})$ returned from Mosek, SDPA, and SDPT3,
respectively. We note that Algorithm 6 returned $z_{\rm pro}$ such that
$\|z_{\rm pro}\|>>\|z_{\rm def}\|$ for hinf2. Since hinf2 is included in the
well-posed group, in Table 8, the values in rows 2, 4, and 6 of the “Well-
posed group” column are significantly different from the values in rows 1, 3
and 5 of the same column, respectively. The average values of $f(z_{\rm def})$
and $f(z_{\rm pro})$ for the well-posed group, excluding hinf2, are summarized
in the second column. From Table 8, we can see that the average values of
$f(z_{\rm def})$ and $f(z_{\rm pro})$ in the third column are greater than in
the second column. Thus, the dual solutions $z_{\rm def}$ and $z_{\rm pro}$
obtained by the solvers and Algorithm 6 for the ill-posed group are likely to
include a reducing direction $(f,-\mathcal{A}^{*}f)$ such that the value of
$\|-\mathcal{A}^{*}f\|$ is very large. Furthermore, we can expect that the
dual optimal solution set of hinf2 includes solutions $(y_{min},z_{min})$ and
$(y_{max},z_{max})$ such that $\|z_{min}\|<<\|z_{max}\|$.

Table 8: Comparison of the average values of $f(z)=\frac{\|z\|}{1+\displaystyle\max_{i}|c_{i}|+\displaystyle\max_{i,j}|A_{ij}|}$ for $z_{\rm def}$ and $z_{\rm pro}$ Method | Well-posed group | Well-posed group excluding hinf2 | Ill-posed group
---|---|---|---
Mosek | 1.04e+03 | 6.85e+02 | 6.75e+05
Mosek + Algorithm 6 | 3.32e+03 | 6.85e+02 | 2.63e+10
SDPA | 1.04e+03 | 6.85e+02 | 8.81e+04
SDPA + Algorithm 6 | 3.32e+03 | 6.85e+02 | 2.74e+11
SDPT3 | 6.70e+02 | 6.85e+02 | 4.66e+08
SDPT3 + Algorithm 6 | 3.33e+03 | 6.85e+02 | 7.87e+09
(a) The values of $\log_{10}{\rm err_{1}}(x,y,z)$ |  (b) The values of $\log_{10}{\rm err_{3}}(x,y,z)$
---|---
(c) The values of $\log_{10}|{\rm err_{5}}(x,y,z)|$ |  (d) The values of $\log_{10}|{\rm err_{6}}(x,y,z)|$
Figure 2: Numreical results for the well-posed group with Mosek (a) The values of $\log_{10}{\rm err_{1}}(x,y,z)$ |  (b) The values of $\log_{10}{\rm err_{3}}(x,y,z)$
---|---
(c) The values of $\log_{10}|{\rm err_{5}}(x,y,z)|$ |  (d) The values of $\log_{10}|{\rm err_{6}}(x,y,z)|$
Figure 3: Numreical results for the well-posed group with SDPA (a) The values of $\log_{10}{\rm err_{1}}(x,y,z)$ |  (b) The values of $\log_{10}{\rm err_{3}}(x,y,z)$
---|---
(c) The values of $\log_{10}|{\rm err_{5}}(x,y,z)|$ |  (d) The values of $\log_{10}|{\rm err_{6}}(x,y,z)|$
Figure 4: Numreical results for the well-posed group with SDPT3 (a) The values of $\log_{10}{\rm err_{1}}(x,y,z)$ |  (b) The values of $\log_{10}{\rm err_{3}}(x,y,z)$
---|---
(c) The values of $\log_{10}|{\rm err_{5}}(x,y,z)|$ |  (d) The values of $\log_{10}|{\rm err_{6}}(x,y,z)|$
Figure 5: Numreical results for the ill-posed group with Mosek (a) The values of $\log_{10}{\rm err_{1}}(x,y,z)$ |  (b) The values of $\log_{10}{\rm err_{3}}(x,y,z)$
---|---
(c) The values of $\log_{10}|{\rm err_{5}}(x,y,z)|$ |  (d) The values of $\log_{10}|{\rm err_{6}}(x,y,z)|$
Figure 6: Numreical results for the ill-posed group with SDPA (a) The values of $\log_{10}{\rm err_{1}}(x,y,z)$ |  (b) The values of $\log_{10}{\rm err_{3}}(x,y,z)$
---|---
(c) The values of $\log_{10}|{\rm err_{5}}(x,y,z)|$ |  (d) The values of $\log_{10}|{\rm err_{6}}(x,y,z)|$
Figure 7: Numreical results for the ill-posed group with SDPT3

Finally, let us compare the execution time of Algorithm 6 for each instance.
Figures 8, 9 and 10 show the execution time of Algorithm 6 using $(x_{\rm
def},y_{\rm def},z_{\rm def})$ returned from Mosek, SDPA and SDPT3,
respectively. In these figures, the horizontal axis shows instances sorted in
ascending order of problem size. In this study, we defined the problem size as
$m\times d$, where $m$ is the number of constraints and $d$ is the dimension
of the Euclidean space $\mathbb{E}$ corresponding to $\mathcal{K}$. The line
graph shows the execution time of Algorithm 6, and the bar graph shows the
problem size. From Figures 8-10, we can observe the following:

  * •

As the problem size increases, the execution time of Algorithm 6 tends to
increase.

  * •

There was a variation in the execution time of Algorithm 6, even for instances
of similar problem size.

The above results were observed due to the structure of the projective
rescaling method. The computational cost of the projection and rescaling
algorithm proposed in [12] is

$\mathcal{O}\left(-\frac{r}{\log\xi}\log\left(\frac{1}{\varepsilon}\right)\left(m^{3}+m^{2}d+\frac{1}{\xi^{2}}p^{2}r_{\max}^{2}\left(\max\left(C^{\rm
sd},md\right)\right)\right)\right)$ (7)

where $C^{\rm sd}$ is the computational cost of the spectral decomposition. In
(7), the maximum number of iterations of the main algorithm, the computational
cost of the projection matrix, the maximum number of iterations of the basic
procedure and the computational cost per iteration of the basic procedure
corresponds to $-\frac{r}{\log\xi}\log\left(\frac{1}{\varepsilon}\right)$,
$m^{3}+m^{2}d$, $\frac{1}{\xi^{2}}p^{2}r_{\max}^{2}$ and $\max\left(C^{\rm
sd},md\right)$, respectively. Since the increase in problem size will increase
the time required to compute the projection matrix, the execution time of
Algorithm 6 in Figures 8-10 tends to increase monotonically. Here, we note
that the main algorithm and the basic procedure can be terminated in less than
their maximum number of iterations. The main algorithm (and the basic
procedure) might terminate in a significantly different number of iterations
for problems with similar problem sizes. Therefore, there is expected to be
some variation in the execution time of Algorithm 6 for problems with similar
problem sizes.

Figure 8: Execution time of Algorithm 6 using $(x_{\rm def},y_{\rm def},z_{\rm
def})$ returned from Mosek Figure 9: Execution time of Algorithm 6 using
$(x_{\rm def},y_{\rm def},z_{\rm def})$ returned from SDPA Figure 10:
Execution time of Algorithm 6 using $(x_{\rm def},y_{\rm def},z_{\rm def})$
returned from SDPT3

###  4.3 Experiments to determine feasibility status

The numerical results in the previous section showed that Algorithm 6 obtained
more accurate approximate optimal solutions for the well-posed group than the
solvers. On the other hand, the numerical results for the ill-posed group
showed that Algorithm 6 failed to detect that the primal problem is in weak
status. Algorithm 6, theoretically, can obtain a reducing direction for (P)
(or (D)) if (P) (or (D)) is in weak status, but Algorithm 6 returned not a
reducing direction for (P) but an approximate interior feasible solution to
(P) for the instances of the ill-posed group. Detecting the feasibility status
of (P) and (D) and finding the reducing directions for (P) and (D) can be
replaced by solving a certain SDP whose primal and dual problems are strongly
feasible, [16]. Thus, using the results of [16], we tested whether Algorithm 6
can detect the feasibility status of SDP more accurately than the solvers. We
conducted this experiment with the ill-posed instances of SDPLIB and the
weakly infeasible cases [14]. These instances have the primal problems in weak
status.

####  4.3.1 Preparation to explain the flow of the experiment

Before going to the outline of our experiments, we derive the
$(P^{D}_{\mathcal{K}})$ and $(D^{D}_{\mathcal{K}})$, which is closely related
to the feasibility status of (P), by using the following lemma.

######  Lemma 4.1 (Lemma 3.4 in [16]).

For (P) and (D), consider the following pair of primal and dual problems.

$\begin{array}[]{lll}(P^{D}_{\mathcal{K}})&{\rm inf}_{x,t,w}&t\\\ &{\rm
sub.to}&-\langle c,x-te\rangle+t-w=0\\\ &&\langle ex\rangle+w=1\\\
&&\mathcal{A}x-t\mathcal{A}e=0\\\
&&(x,t,w)\in\mathcal{K}\times\mathbb{R}_{+}\times\mathbb{R}_{+}\end{array}$

$\begin{array}[]{lll}(D^{D}_{\mathcal{K}})&{\rm
sup}_{y_{1},y_{2},y_{3}}&y_{2}\\\ &{\rm
sub.to}&cy_{1}-ey_{2}-\mathcal{A}^{*}y_{3}\in\mathcal{K}\\\
&&1-y_{1}(1+\langle c,e\rangle)+\langle e,\mathcal{A}^{*}y_{3}\rangle\geq 0\\\
&&y_{1}-y_{2}\geq 0\end{array}$

The following properties hold.

  1. 1.

$(P^{D}_{\mathcal{K}})$ and $(D^{D}_{\mathcal{K}})$ are strongly feasible.

Let $(x^{*},t^{*},w^{*})$ be an optimal solution to $(P^{D}_{\mathcal{K}})$
and $(y_{1}^{*},y_{2}^{*},y_{3}^{*})$ be an optimal solution to
$(D^{D}_{\mathcal{K}})$.

  2. 2.

The optimal value is zero if and only if (D) is not strongly feasible. In this
case, one of the two alternatives below must hold:

    1. (a)

$x^{*}$ is an improving ray of (P), or

    2. (b)

$x^{*}$ is a reducing direction for (D).

  3. 3.

The optimal value is positive if and only if (D) is strongly feasible.

###### Proof.

See Lemma 3.4 in [16]. ∎

Note that [16] showed that Lemma 4.1 holds not only for the identity element
$e$ of $\mathcal{K}$ but also for any interior point of $\mathcal{K}$. Similar
to Lemma 4.1, Proposition 4.2, which is not in [16], can be easily proved.

######  Proposition 4.2.

For (P) and (D), consider the following pair of primal and dual problems.

$\begin{array}[]{lll}(P^{P}_{\mathcal{K}})&{\rm inf}_{y,z,t,w}&t\\\ &{\rm
sub.to}&b^{\top}y+t-w=0\\\ &&\langle e,z\rangle+w=1\\\
&&-\mathcal{A}^{*}y-z+te=0\\\
&&(y,s,t,w)\in\mathbb{R}^{m}\times\mathcal{K}\times\mathbb{R}_{+}\times\mathbb{R}_{+}\\\
\\\ (D^{P}_{\mathcal{K}})&{\rm sup}_{x_{1},x_{2},x_{3}}&x_{2}\\\ &{\rm
sub.to}&\mathcal{A}x_{3}-x_{1}b=0\\\ &&1-x_{1}-\langle e,x_{3}\rangle\geq 0\\\
&&x_{1}-x_{2}\geq 0\\\ &&x_{3}-x_{2}e\in\mathcal{K}\end{array}$

The following properties hold.

  1. 1.

Both $(P^{P}_{\mathcal{K}})$ and $(D^{P}_{\mathcal{K}})$ are strongly
feasible.

  2. 2.

Let $(y^{*},z^{*},t^{*},w^{*})$ be a primal optimal solution. The optimal
value is zero if and only if (P) is not strongly feasible. Moreover, if the
optimal value is zero, $(y^{*},z^{*})$ is an improving ray of (D) or a
reducing direction for (P).

  3. 3.

Let $(x_{1}^{*},x_{2}^{*},x_{3}^{*})$ be a dual optimal solution. If the
optimal value is positive, $\frac{1}{x_{1}^{*}}x_{3}^{*}$ is an interior
feasible solution to (P).

###### Proof.

(1): Let $(y^{0},z^{0},t^{0},w^{0})=\left(0,\frac{1}{\langle
e,e\rangle+1}e,\frac{1}{\langle e,e\rangle+1},\frac{1}{\langle
e,e\rangle+1}\right)$ and $(x_{1}^{0},x_{2}^{0},x_{3}^{0})=(0,-1,0)$. Then
$(y^{0},z^{0},t^{0},w^{0})$ and $(x_{1}^{0},x_{2}^{0},x_{3}^{0})$ are interior
feasible solutions to $(P^{P}_{\mathcal{K}})$ and $(D^{P}_{\mathcal{K}})$,
respectively.
(2): If $t^{*}=0$, then $b^{\top}y^{*}=w^{*}\geq 0$ and
$-\mathcal{A}^{*}y^{*}=z^{*}\in\mathcal{K}$ hold, which implies that (P) is
not strongly feasible by Proposition 2.4. Conversely, if (P) is not strongly
feasible, there exists $(y,z)\in\mathbb{R}^{m}\times\mathcal{K}$ such that
$-\mathcal{A}^{*}y=z$, $\mathcal{A}^{*}y\neq 0$ and $b^{\top}y\geq 0$ by
Proposition 2.4. Let $k=\langle e,z\rangle+b^{\top}y$. Since $b^{\top}y\geq
0$, $z\in\mathcal{K}$ and $z\neq 0$ hold, $k$ is positive. Letting
$(\bar{y},\bar{z},\bar{t},\bar{w})=(\frac{1}{k}y,\frac{1}{k}z,0,\frac{1}{k}b^{\top}y)$,
we can easily see that $(\bar{y},\bar{z},\bar{t},\bar{w})$ is an optimal
solution for $(P^{P}_{\mathcal{K}})$, which implies that the optimal value is
zero.

Suppose that $t^{*}=0$. If $b^{\top}y^{*}=0$, then $z^{*}\neq 0$ holds since
$w^{*}=0$ and hence, $(y^{*},z^{*})$ is a reducing direction for (P). If
$b^{\top}y^{*}>0$, we can easily see that $(y^{*},z^{*})$ is an improving ray
of (D).
(3): If $x_{2}^{*}>0$, then $x_{1}^{*}>0$ and
$x_{3}^{*}\in\mbox{int}\mathcal{K}$ hold. Since $x_{1}^{*}$ and $x_{3}^{*}$
satisfy $\mathcal{A}x_{3}^{*}-x_{1}^{*}b=0$, $\frac{1}{x_{1}^{*}}x_{3}^{*}$ is
an interior feasible solution to (P). ∎

With some modifications to $(P^{P}_{\mathcal{K}})$ and
$(D^{P}_{\mathcal{K}})$, we have optimization problems that Algorithm 6 can
handle well.

######  Proposition 4.3.