Kernel-function-based primal-dual interior-point methods for convex quadratic optimization over symmetric cone

Abstract

In this paper, we give a unified computational scheme for the complexity analysis of kernel-function-based primal-dual interior-point methods for convex quadratic optimization over symmetric cone. By using Euclidean Jordan algebras, the currently best-known iteration bounds for large- and small-update methods are derived, namely, $O(\sqrt{r}logrlog\frac{r}{\epsilon })$ and $O(\sqrt{r}log\frac{r}{\epsilon })$, respectively. Furthermore, this unifies the analysis for a wide class of conic optimization problems.

MSC:90C25, 90C51.

1 Introduction

Since the groundbreaking paper of Karmarkar, many researchers have proposed and analyzed various interior-point methods (IPMs) for linear optimization (LO) and a large amount of results have been reported [1–4]. However, there is a gap between the practical behavior of the IPMs and the theoretical performance results. The so-called small-update IPMs enjoy the best-known worst-case iteration bound $O(\sqrt{n}log\frac{n}{\epsilon })$ but their performance in computational practice is poor. In practice, however, the so-called large-update IPMs are much more efficient than small-update IPMs but with relatively weak theoretical result $O(nlog\frac{n}{\epsilon })$. Recently, Peng et al. [5] introduced so-called self-regular barrier functions for primal-dual IPMs for LO, the iteration bound for large-update methods for LO was improved from $O(nlog\frac{n}{\epsilon })$ to $O(\sqrt{n}lognlog\frac{n}{\epsilon })$, which almost closes the gap between the iteration bounds for large- and small-update methods. Bai et al. [6] presented a large class of eligible kernel functions, which is fairly general and includes the classical logarithmic function and the self-regular functions, as well as many non-self-regular functions as special cases. The best-known iteration bounds for LO obtained are as good as the ones in [5] for appropriate choices of the eligible kernel functions. For some other related kernel-based IPMs we refer to [7–27].

In this paper, we present a unified kernel-function approach to primal-dual IPMs for convex quadratic optimization over symmetric cone (CQSCO), which is a generalization of symmetric cone optimization (SCO) (when $\mathcal{Q}=0$), which contains LO, second-order cone optimization (SOCO) and semidefinite optimization (SDO) as special case. CQSCO also includes convex quadratic optimization (CQO) and convex quadratic semidefinite optimization (CQSDO). Let $(\mathcal{V},\circ )$ be an n-dimensional Euclidean Jordan algebra (EJA) with rank r equipped with the standard inner product $〈x,s〉=tr(x\circ s)$, and $\mathcal{K}$ be the corresponding symmetric cone. The primal problem of CQSCO is given by

$$\begin{array}{r}minf(x)=\frac{1}{2}〈x,\mathcal{Q}(x)〉+〈c,x〉\\ \text{s.t. }\mathcal{A}(x)=b,x\in \mathcal{K},\end{array}$$

(P)

where $c\in \mathcal{V}$ and $b\in {\mathbf{R}}^m$ are given data, $\mathcal{A}:\mathcal{V}\to {\mathbf{R}}^m$ is a given linear map, and $\mathcal{Q}$ is a given self-adjoint positive semidefinite (with respect to $〈\cdot ,\cdot 〉$) linear operator on $\mathcal{V}$, i.e., for any $x,s\in \mathcal{V}$, then $〈\mathcal{Q}(x),s〉=〈x,\mathcal{Q}(s)〉$ and $〈\mathcal{Q}(x),x〉\ge 0$. The dual problem of (P) is given by

$$\begin{array}{r}max-\frac{1}{2}〈x,\mathcal{Q}(x)〉+b^{T}y\\ \text{s.t. }{\mathcal{A}}^T(y)+s=\mathrm{\nabla }f(x)=\mathcal{Q}(x)+c,s\in \mathcal{K},\end{array}$$

(D)

where ${\mathcal{A}}^T$ is the adjoint of $\mathcal{A}$. Many researchers have studied CQSCO and achieved plentiful and beautiful results. For an overview of these results we refer to [28–34].

Without loss of generality, we assume that the linear map $\mathcal{A}$ is surjective, which implies that $\mathcal{A}{\mathcal{A}}^T$ is nonsingular. Furthermore, we also assume that both (P) and (D) satisfy the interior-point condition (IPC), i.e., there exists $(x^0,y^0,s^0)$ such that

$$\mathcal{A}\left(x^0\right)=b,x^0\in int\mathcal{K},{\mathcal{A}}^T\left(y^0\right)+s^0-\mathcal{Q}\left(x^0\right)=c,s^0\in int\mathcal{K}.$$

The perturbed Karush-Kuhn-Tucker optimality conditions for the problems (P) and (D) are given as follows:

$$\begin{array}{c}\mathcal{A}(x)=b,x\in \mathcal{K},\hfill \\ {\mathcal{A}}^T(y)+s-\mathcal{Q}(x)=c,s\in \mathcal{K},\hfill \\ x\circ s=\mu e,\hfill \end{array}$$

(1)

where μ is a positive parameter that is to be driven to zero explicitly. Since the IPC holds and $\mathcal{A}$ is surjective, the parameterized system (1) has a unique solution $(x(\mu ),y(\mu ),s(\mu ))$ for each $\mu >0$, and we call $x(\mu )$ the μ-center of $(P)$ and $(y(\mu ),s(\mu ))$ the μ-center of (D). The set of μ-centers gives a homotopy path (with μ running through all the positive real numbers), which is called the central path. If $\mu \to 0$, then the limit of the central path exists and since the limit points satisfy the complementarity condition, i.e., $x\circ s=0$, it naturally yields an optimal solution for (P) and (D) (see, e.g., [29, 35]).

IPMs follow the central path approximately and find an approximate solution of the underlying problems (P) and (D) as μ go to zero. Just like the case of a linear SDO, linearizing the third equation in (1) may not lead to an element in $\mathcal{V}$. Thus it is necessary to symmetrize that equation before linearizing it. For this purpose, we can apply the following scaling scheme (cf. Lemma 28 in [36]): Let $u\in int\mathcal{K}$. Then

$$x\circ s=\mu e\iff P(u)x\circ P\left(u^{-1}\right)s=\mu e.$$

Thus, we replace the third equation of the system (1) by

$$P(u)x\circ P\left(u^{-1}\right)s=\mu e.$$

Applying Newton’s method, and neglecting the term $P(u)\mathrm{\Delta }x\circ P(u^{-1})\mathrm{\Delta }s$, we have

$$\begin{array}{c}\mathcal{A}(\mathrm{\Delta }x)=0,\hfill \\ {\mathcal{A}}^T(\mathrm{\Delta }y)+\mathrm{\Delta }s-\mathcal{Q}(\mathrm{\Delta }x)=0,\hfill \\ P(u)x\circ P\left(u^{-1}\right)\mathrm{\Delta }s+P\left(u^{-1}\right)s\circ P(u)\mathrm{\Delta }x=\mu e-P(u)x\circ P\left(u^{-1}\right)s.\hfill \end{array}$$

(2)

The appropriate choices of u that lead to obtaining the unique search directions from the above system are called commutative class of search directions (see, e.g., [36]). In this paper, we consider the so-called NT-scaling scheme, the resulting direction is called NT search direction. This scaling scheme was first proposed by Nesterov and Todd [37, 38] for self-scaled cones and then adapted by Faybusovich [35, 39] for symmetric cones.

Lemma 1.1 (Lemma 3.2 in [39])

Let $x,s\in int\mathcal{K}$. Then there exists a unique $w\in int\mathcal{K}$ such that

$$x=P(w)s.$$

Moreover,

$$w=P{(x)}^{\frac{1}{2}}{\left(P\left(x^{\frac{1}{2}}\right)s\right)}^{-\frac{1}{2}}[=P\left(s^{-\frac{1}{2}}\right){\left(P\left(s^{\frac{1}{2}}\right)x\right)}^{\frac{1}{2}}].$$

The point w is called the scaling point of x and s (in this order). As a consequence there exists $\tilde{v}\in int\mathcal{K}$ such that

$$\tilde{v}=P{(w)}^{-\frac{1}{2}}x=P{(w)}^{\frac{1}{2}}s.$$

(3)

Let $u=w^{-\frac{1}{2}}$, where w is the NT-scaling point of x and s. We define

$$v:=\frac{P{(w)}^{-\frac{1}{2}}x}{\sqrt{\mu }}[=\frac{P{(w)}^{\frac{1}{2}}s}{\sqrt{\mu }}],$$

(4)

and the scaled search directions as follows:

$$d_x:=\frac{P{(w)}^{-\frac{1}{2}}\mathrm{\Delta }x}{\sqrt{\mu }}\text{and}{d}_s:=\frac{P{(w)}^{\frac{1}{2}}\mathrm{\Delta }s}{\sqrt{\mu }}.$$

(5)

It follows from (4) and (5) that

$$\begin{array}{c}\overline{\mathcal{A}}(d_x)=0,\hfill \\ {\overline{\mathcal{A}}}^T(\mathrm{\Delta }y)+d_s-\overline{\mathcal{Q}}(d_x)=0,\hfill \\ d_x+d_s=v^{-1}-v,\hfill \end{array}$$

(6)

where $\overline{\mathcal{A}}=\frac{\mathcal{A}P{(w)}^{\frac{1}{2}}}{\sqrt{\mu }}$ and $\overline{\mathcal{Q}}=P{(w)}^{\frac{1}{2}}\mathcal{Q}P{(w)}^{\frac{1}{2}}$. We can easily verify that the system (6) has a unique solution (see, e.g., [29, 35]).

In this paper, we replace the right-hand side of the third equation in (6) by $-{\psi }^{\prime }(v)$, i.e., $-\mathrm{\nabla }\mathrm{\Psi }(v)$, as defined by (26) (see Section 3), where $\psi (t)$ is any eligible kernel function. This yields the following system:

$$\begin{array}{c}\overline{\mathcal{A}}(d_x)=0,\hfill \\ {\overline{\mathcal{A}}}^T(\mathrm{\Delta }y)+d_s-\overline{\mathcal{Q}}(d_x)=0,\hfill \\ d_x+d_s=-{\psi }^{\prime }(v).\hfill \end{array}$$

(7)

Since (7) has the same matrix of coefficients as (6), also (7) has a unique solution.^a It follows that the eligible kernel function ${\psi }^{\prime }(t)$ determines in a natural way search directions for an interior-point algorithm.

The new search directions $d_x$ and $d_s$ are computed by solving (7), thus Δx and Δs are obtained from (5). If $(x,y,s)\ne (x(\mu ),y(\mu ),s(\mu ))$, then $(\mathrm{\Delta }x,\mathrm{\Delta }y,\mathrm{\Delta }s)$ is nonzero. The new iteration point is obtained according to

$$x_{+}:=x+\alpha \mathrm{\Delta }x,y_{+}:=y+\alpha \mathrm{\Delta }y\text{and}{s}_{+}:=s+\alpha \mathrm{\Delta }s.$$

(8)

Similarly to the LO case, we require that the step size α should be taken so that the proximity measure function $\mathrm{\Psi }(v)$ decreases sufficiently. A default bound for such a step size α will be given later by (38).

Furthermore, we can conclude that

$$x\circ s=\mu e\iff v=e\iff \mathrm{\nabla }\mathrm{\Psi }(v)=0\iff \mathrm{\Psi }(v)=0.$$

(9)

Hence, the value of $\mathrm{\Psi }(v)$ can be considered as a measure for the distance between the given iterate $(x,y,s)$ and the corresponding μ-center $(x(\mu ),y(\mu ),s(\mu ))$.

The algorithm considered in this paper is described in Figure 1.

Figure 1

Algorithm.

Given any eligible kernel function $\psi (t)$, the parameters τ, θ and the step size α should be chosen in such a way that the algorithm is ‘optimized’ in the sense that the number of iterations required by the algorithm is as small as possible. We will prove that the resulting iteration bounds depend on the eligible kernel functions in Section 5.

The purpose of the paper is to propose a unified analysis of kernel-function-based primal-dual IPMs for CQSCO and give a general scheme on how to calculate the iteration bounds for the entire class of eligible kernel functions. The obtained complexity results match the best-known iteration bounds known for large-update methods, $O(\sqrt{r}logrlog\frac{r}{\epsilon })$ and small-update methods, $O(\sqrt{r}log\frac{r}{\epsilon })$. The order of the iteration bounds are derived as good as the ones for the LO case except that n is replaced by r, the rank of EJA. Although expected, these results were not obvious and, at certain steps of the analysis, they were not trivial and/or straightforward generalization of the LO case. Furthermore, this unifies the analysis for a wide class of conic optimization problems, which includes LO, CQO, SOCO, SDO, CQSDO, SCO and so on.

The outline of the paper is as follows. In Section 2, we provide some basic concepts and useful results on EJAs and symmetric cones. In Section 3, we recall and develop some useful properties of the eligible kernel functions and the corresponding barrier functions. In Section 4, we uniformly analyze the primal-dual IPMs for CQSCO. In Section 5, we derive the complexity bounds for large- and small-update methods. In Section 6, we report some preliminary numerical experiments. Finally, some conclusions and remarks are made in Section 7.

The following notations are used throughout the paper. ${\mathbf{R}}^n$, ${\mathbf{R}}_{+}^n$, and ${\mathbf{R}}_{++}^n$ denote the set of all vectors (with n components), the set of non-negative vectors and the set of positive vectors, respectively. ${\mathbf{R}}^{m\times n}$ is the space of all $m\times n$ matrices. ${\mathbf{S}}^n$, ${\mathbf{S}}_{+}^n$ and ${\mathbf{S}}_{++}^n$ denote the cones of symmetric, symmetric positive semidefinite and symmetric positive definite $n\times n$ matrices, respectively. We use the matrix inner product $A•B=tr(A^{T}B)$, i.e., the trace of the matrix $A^{T}B$. The largest eigenvalue and the smallest eigenvalue of x are defined by ${\lambda }_{\max }(x)$ and ${\lambda }_{\min }(x)$, respectively. The Löwner partial ordering ‘${⪰}_{\mathcal{K}}$’ of $\mathcal{V}$ defined by a symmetric cone $\mathcal{K}$ is defined by $x{⪰}_{\mathcal{K}}s$ if $x-s\in \mathcal{K}$. The interior of $\mathcal{K}$ is denoted as $int\mathcal{K}$ and we write $x{\succ }_{\mathcal{K}}s$ if $x-s\in int\mathcal{K}$. Finally, if $g(x)\ge 0$ is a real-valued function of a real non-negative variable, the notation $g(x)=O(x)$ means that $g(x)\le \overline{c}x$ for some positive constant $\overline{c}$, and $g(x)=\mathrm{\Theta }(x)$ that $c_{1}x\le g(x)\le c_{2}x$ for the two positive constants $c_1$ and $c_2$.

2 Preliminaries

For any $x,y\in \mathcal{V}$, the Lyapunov transformation $L(x)$ and the quadratic representation $P(x)$ are given by

$$L(x)y:=x\circ y$$

(10)

and

$$P(x):=2L{(x)}^2-L\left(x^2\right),$$

(11)

where $L{(x)}^2=L(x)L(x)$, respectively.

For any EJA $\mathcal{V}$, the corresponding cone of squares

$$\mathcal{K}(\mathcal{V}):=\{x^2:x\in \mathcal{V}\}$$

(12)

is indeed a symmetric cone (cf. Theorem III.2.1 in [40]). In the sequel, $\mathcal{K}$ will always denote a symmetric cone, and $\mathcal{V}$ an EJA with $rank(\mathcal{V})=r$ for which $\mathcal{K}$ is its cone of squares.

The following theorem gives an important decomposition, the spectral decomposition, on the space $\mathcal{V}$.

Theorem 2.1 (Theorem III.1.2 in [40])

Let $x\in \mathcal{V}$. Then there exists a Jordan frame $\{c_1,\dots ,c_r\}$ and real numbers ${\lambda }_1(x),\dots ,{\lambda }_r(x)$ such that

$$x=\sum _{i=1}^r{\lambda }_i(x)c_i.$$

(13)

The numbers ${\lambda }_i(x)$ (with their multiplicities) are called the eigenvalues of x. Furthermore, the trace and the determinant of x are given by

$$tr(x)=\sum _{i=1}^r{\lambda }_i(x)\mathit{\text{and}}det(x)=\prod _{i=1}^r{\lambda }_i(x),$$

respectively.

Let $x\in \mathcal{K}$ with the spectral decomposition given by (13), the vector-valued function $\psi (x)$ is defined by

$$\psi (x):=\psi ({\lambda }_1(x))c_1+\cdots +\psi ({\lambda }_r(x))c_r.$$

(14)

Furthermore, if $\psi (t)$ is differentiable, the derivative ${\psi }^{\prime }(t)$ exists, and we also have the vector-valued function ${\psi }^{\prime }(x)$, namely

$${\psi }^{\prime }(x)={\psi }^{\prime }({\lambda }_1(x))c_1+\cdots +{\psi }^{\prime }({\lambda }_r(x))c_r.$$

(15)

It should be noted that ${\psi }^{\prime }(x)$ is just a vector-valued function induced by the derivative ${\psi }^{\prime }(t)$ of the function $\psi (t)$ rather than the derivative of the vector-valued function $\psi (x)$ defined by (14).

The following theorem provides another important decomposition, the Peirce decomposition, on the space $\mathcal{V}$.

Theorem 2.2 (Theorem IV.2.1 in [40])

Let $x\in \mathcal{V}$ with the spectral decomposition given by (13). Then we have

$$\mathcal{V}=\underset{i\le j}{⨁}{\mathcal{V}}_{ij},$$

where

$${\mathcal{V}}_{ii}:=\{x|x\circ c_i=x\}\mathit{\text{and}}{\mathcal{V}}_{ij}:=\left\{x\right|x\circ c_i=\frac{1}{2}x=x\circ c_j\},1\le i<j\le r,$$

are Peirce spaces of $\mathcal{V}$. Then, for any $x\in \mathcal{V}$, there exist $x_i\in \mathbf{R}$, $c_i\in {\mathcal{V}}_{ii}$, and $x_{ij}\in {\mathcal{V}}_{ij}$ ($i<j$) such that

$$x=\sum _{i=1}^r{x}_i{c}_i+\sum _{i<j}{x}_{ij}.$$

For any $x,s\in \mathcal{V}$, we define

$$〈x,s〉:=tr(x\circ s),$$

(16)

and we refer to it as the trace inner product. The Frobenius norm induced by this trace inner product, namely ${\parallel \cdot \parallel }_F$, is defined by

$${\parallel x\parallel }_F:=\sqrt{〈x,x〉}.$$

(17)

Thus, we have

$${\parallel x\parallel }_F=\sqrt{tr\left(x^2\right)}=\sqrt{\sum _{i=1}^r{\lambda }_i^2(x)}.$$

(18)

Furthermore, we can easily verify that

$$|{\lambda }_{\min }(x)|\le {\parallel x\parallel }_F\text{and}|{\lambda }_{\max }(x)|\le {\parallel x\parallel }_F.$$

(19)

Lemma 2.3 (Lemma 14 in [36])

Let $x,s\in \mathcal{V}$. Then

$${\lambda }_{\min }(x+s)\ge {\lambda }_{\min }(x)+{\lambda }_{\min }(s)\ge {\lambda }_{\min }(x)-{\parallel s\parallel }_F$$

and

$${\lambda }_{\max }(x+s)\le {\lambda }_{\max }(x)+{\lambda }_{\max }(s)\le {\lambda }_{\max }(x)+{\parallel s\parallel }_F.$$

Let $f:D\to \mathbf{R}$ be a univariate function on the open set $D\subseteq \mathbf{R}$ that is differentiable or even continuously differentiable if necessary, and $x={\sum }_{i=1}^r{\lambda }_i(x)c_i$ be the spectral decomposition of $x\in \mathcal{V}$ with respect to the Jordan frame $\{c_1,\dots ,c_r\}$. The real-valued separable spectral function $F:\mathcal{V}\to \mathbf{R}$ and the vector-valued separable spectral function $G:\mathcal{V}\to \mathcal{V}$ are defined by

$$F(x):=\sum _{i=1}^{r}f({\lambda }_i(x))$$

(20)

and

$$G(x):=\sum _{i=1}^{r}f({\lambda }_i(x))c_i,$$

(21)

respectively.

The following two theorems give explicitly the first derivatives of $F(x)$ and $G(x)$, respectively.

Theorem 2.4 (Theorem 38 in [41])

Let f is continuously differentiable in D. Then $F(x)$ is continuously differentiable at x and

$$D_{x}F(x)=\sum _{i=1}^r{f}^{\prime }({\lambda }_i(x))c_i.$$

Theorem 2.5 (Lemma 1 in [42])

Let f is a continuously differentiable in D. Then $G(x)$ is continuously differentiable at x and

$$D_{x}G(x)=\sum _{i=1}^r{f}^{\prime }({\lambda }_i(x))x_i{c}_i+\underset{{\lambda }_i(x)={\lambda }_j(x)}{\sum _{i<j}}{f}^{\prime }({\lambda }_i(x))x_{ij}+\underset{{\lambda }_i(x)\ne {\lambda }_j(x)}{\sum _{i<j}}\frac{f({\lambda }_i(x))-f({\lambda }_j(x))}{{\lambda }_i(x)-{\lambda }_j(x)}{x}_{ij},$$

where $1\le i<j\le r$.

3 Properties of the eligible kernel (barrier) functions

We call a univariate $\psi :(0,\mathrm{\infty })\to [0,\mathrm{\infty })$ a kernel function [5] if it satisfies the following three conditions:

$${\psi }^{\prime }(1)=\psi (1)=0,$$

(22a)

$${\psi }^{″}(t)>0,$$

(22b)

$$\underset{t\downarrow 0}{lim}\psi (t)=\underset{t\to \mathrm{\infty }}{lim}\psi (t)=\mathrm{\infty }.$$

(22c)

This means that $\psi (t)$ is strictly convex and minimal at $t=1$, with $\psi (1)=0$. Moreover, (22c) implies that $\psi (t)$ has the barrier property.

In this paper, we consider the so-called eligible kernel function [6], i.e., the kernel function satisfies four of the following five conditions, namely the first and the last three conditions:

$$t{\psi }^{″}(t)+{\psi }^{\prime }(t)>0,t<1,$$

(23a)

$$t{\psi }^{″}(t)-{\psi }^{\prime }(t)>0,t>1,$$

(23b)

$${\psi }^{‴}(t)<0,t>0,$$

(23c)

$$2{\psi }^{″}{(t)}^2-{\psi }^{\prime }(t){\psi }^{‴}(t)>0,t<1,$$

(23d)

$${\psi }^{″}(t){\psi }^{\prime }(\beta t)-\beta {\psi }^{\prime }(t){\psi }^{″}(\beta t)>0,t>1,\beta >1.$$

(23e)

Note that the first four conditions are logically independent, and the fifth condition is a consequence of (23b) and (23c). Since (23b) is much simpler to check than (23e), in many cases it is easy to know that $\psi (t)$ is eligible if it satisfies the first four conditions [6].

The following lemma is cited from [6] to state the exponential convexity, which plays an important role in the analysis of kernel-function-based primal-dual IPMs [5, 6].

Lemma 3.1 (Lemma 2.1 in [6])

Let $t_1>0$ and $t_2>0$. Then

$$\psi (\sqrt{t_1{t}_2})\le \frac{1}{2}(\psi (t_1)+\psi (t_2)).$$

Now, we define the barrier function $\mathrm{\Psi }(v):int\mathcal{K}\to {\mathbf{R}}_{+}$ as

$$\mathrm{\Psi }(x,s,\mu ):=\mathrm{\Psi }(v):=tr(\psi (v)).$$

(24)

It follows from Theorem 2.1 and (14) that

$$\mathrm{\Psi }(v)=\sum _{i=1}^r\psi ({\lambda }_i(v)).$$

(25)

Furthermore, we have, by Theorem 2.4,

$$\mathrm{\nabla }\mathrm{\Psi }(v)={\psi }^{\prime }(v):={\psi }^{\prime }({\lambda }_1(v))c_1+\cdots +{\psi }^{\prime }({\lambda }_r(v))c_r,$$

(26)

where $\mathrm{\nabla }\mathrm{\Psi }(v)$ denotes the derivative of the barrier function $\mathrm{\Psi }(v)$.

As a consequence of Lemma 3.1, we have the following important result.

Theorem 3.2 (Theorem 4.3.2 in [23])

Let $x,s\in int\mathcal{K}$. Then

$$\mathrm{\Psi }\left({\left(P{(x)}^{1/2}s\right)}^{1/2}\right)\le \frac{1}{2}(\mathrm{\Psi }(x)+\mathrm{\Psi }(s)).$$

Note that during the course of the algorithm the largest values of $\mathrm{\Psi }(v)$ occur just after the update of μ. So next we derive an estimate for the effect of a μ-update on the value of $\mathrm{\Psi }(v)$. It follows from (24) that

$$\mathrm{\Psi }(\beta v)=\sum _{i=1}^r\psi (\beta {\lambda }_i(v)).$$

By applying Theorem 3.2 in [6], with x being the vector in ${\mathbf{R}}^r$ consisting of all the eigenvalues of the symmetric cone v, the theorem below immediately follows.

Theorem 3.3 Let $v\in int\mathcal{K}$ and $\beta \ge 1$. Then

$$\mathrm{\Psi }(\beta v)\le r\psi \left(\beta \varrho \left(\frac{\mathrm{\Psi }(v)}{r}\right)\right).$$

Corollary 3.4 Let $0\le \theta <1$ and $v_{+}=\frac{v}{\sqrt{1-\theta }}$. If $\mathrm{\Psi }(v)\le \tau$, then

$$\mathrm{\Psi }(v_{+})\le r\psi \left(\frac{\varrho (\frac{\tau }{r})}{\sqrt{1-\theta }}\right).$$

Proof With $\beta =\frac{1}{\sqrt{1-\theta }}\ge 1$ and $\mathrm{\Psi }(v)\le \tau$, the corollary follows immediately from Theorem 3.3. □

The norm-based proximity measure $\delta (v):int\mathcal{K}\to {\mathbf{R}}_{+}$ is defined by

$$\delta (v):=\frac{1}{2}{\parallel \mathrm{\nabla }\mathrm{\Psi }(v)\parallel }_F.$$

(27)

It follows from (17) and (26) that

$$\delta (v)=\frac{1}{2}{\parallel {\psi }^{\prime }(v)\parallel }_F=\frac{1}{2}\sqrt{\sum _{i=1}^r{\psi }^{\prime }{({\lambda }_i(v))}^2}.$$

(28)

Hence, we can conclude that $\delta (v)\ge 0$, and $\delta (v)=0$ if and only if $\mathrm{\Psi }(v)=0$.

It follows from (25) and (28) that $\delta (v)$ and $\mathrm{\Psi }(v)$ depend only on the eigenvalues ${\lambda }_i(v)$ of the symmetric cone $\mathcal{V}$. This observation makes it possible to apply Theorem 4.8 in [6], with x being the vector in ${\mathbf{R}}^r$ consisting of all the eigenvalues of the symmetric cone v. This gives the following theorem, which yields a lower bound on $\delta (v)$ in terms of $\mathrm{\Psi }(v)$.

Theorem 3.5 Let $v\in int\mathcal{K}$. Then

$$\delta (v)\ge \frac{1}{2}{\psi }^{\prime }\left(\varrho (\mathrm{\Psi }(v))\right).$$

In what follows, we consider the derivatives of the function $\mathrm{\Psi }(x(t))$ with respect to t, where $x(t)=x_0+tu\in int\mathcal{K}$ with $t\in \mathbf{R}$ and $u\in \mathcal{V}$. It follows from Theorem 2.1 and Theorem 2.2 that the spectral decomposition of $x(t)$ with respect to the Jordan frame $\{c_1,\dots ,c_r\}$ can be defined by

$$x(t)=\sum _{i=1}^r{\lambda }_i(x(t))c_i,$$

(29)

and the Peirce decomposition of u can be defined by

$$u=\sum _{i=1}^r{u}_i{c}_i+\sum _{i<j}{u}_{ij}.$$

(30)

From Theorem 2.4 and Theorem 2.5, after some elementary reductions, we can derive the first two derivatives of the general function $\mathrm{\Psi }(x(t))$ with respect to t as follows:

$$D_t\mathrm{\Psi }(x(t))=tr(D_x\mathrm{\Psi }(x(t))\circ x^{\prime }(t))=tr(\sum _{i=1}^r{\psi }^{\prime }\left({\lambda }_i(x(t))\right)c_i\circ u)$$

(31)

and

$$\begin{array}{rl}{D}_t^2\mathrm{\Psi }(x(t))=& \sum _{i=1}^r{\psi }^{″}\left({\lambda }_i(x(t))\right){(u_i)}^2+\underset{{\lambda }_i(x(t))={\lambda }_j(x(t))}{\sum _{i<j}}{\psi }^{″}\left({\lambda }_i(x(t))\right)tr\left({(u_{ij})}^2\right)\\ +\underset{{\lambda }_i(x(t))\ne {\lambda }_j(x(t))}{\sum _{i<j}}\frac{{\psi }^{\prime }({\lambda }_i(x(t)))-{\psi }^{\prime }({\lambda }_j(x(t)))}{{\lambda }_i(x(t))-{\lambda }_j(x(t))}tr\left({(u_{ij})}^2\right).\end{array}$$

(32)

The condition (23c) implies that ${\psi }^{″}(t)$ is monotonically decreasing in $t\in (0,+\mathrm{\infty })$. Under the assumption that $i<j$ implies ${\lambda }_i(x(t))\ge {\lambda }_j(x(t))$, we can conclude that

$$D_t^2\mathrm{\Psi }(x(t))\le \sum _{i=1}^r{\psi }^{″}\left({\lambda }_i(x(t))\right){(u_i)}^2+\sum _{i<j}{\psi }^{″}\left({\lambda }_j(x(t))\right)tr\left({(u_{ij})}^2\right),$$

(33)

which bounds the second-order derivative of $\mathrm{\Psi }(x(t))$ with respect to t (see, e.g., [23]).

4 Analysis of the algorithms

From (8) and (5), after some elementary reductions, we have

$$x_{+}=\sqrt{\mu }P{\left(w^{(j)}\right)}^{1/2}(v+\alpha d_x)\text{and}{s}_{+}=\sqrt{\mu }P{(w)}^{-1/2}(v+\alpha d_s).$$

Note that during an inner iteration the parameter μ is fixed. Hence, after the default step the new scaled vector $v_{+}$ is given by

$$v_{+}=P{(w_{+})}^{-1/2}P{(w)}^{1/2}(v+\alpha d_x)=P{(w_{+})}^{1/2}P{(w)}^{-1/2}(v+\alpha d_s),$$

where, according to Lemma 1.1,

$$w_{+}=P{(x_{+})}^{1/2}\left({\left(P{(x_{+})}^{1/2}{s}_{+}\right)}^{-1/2}\right).$$

To calculate a decrease of the barrier function $\mathrm{\Psi }(v)$ during an inner iteration it is standard to consider a decrease as a function of α defined by

$$f(\alpha ):=\mathrm{\Psi }(v_{+})-\mathrm{\Psi }(v).$$

However, working with $f(\alpha )$ may not be easy because in general $f(\alpha )$ is not convex. Thus, we are searching for the convex function $f_1(\alpha )$ that is an upper bound of $f(\alpha )$ and whose derivatives are easier to calculate than those of $f(\alpha )$. The key element in this process is replacing $v_{+}$ with a similar element that will allow the use of exponential-convexity of the barrier function. By Proposition 5.9.3 in [23], we have

$$v_{+}\sim {(P{(v+\alpha d_x)}^{\frac{1}{2}}(v+\alpha d_s))}^{\frac{1}{2}}$$

and therefore

$$\mathrm{\Psi }(v_{+})=\mathrm{\Psi }\left({(P{(v+\alpha d_x)}^{\frac{1}{2}}(v+\alpha d_s))}^{\frac{1}{2}}\right).$$

Theorem 3.2 implies that

$$\mathrm{\Psi }(v_{+})\le \frac{1}{2}(\mathrm{\Psi }(v+\alpha d_x)+\mathrm{\Psi }(v+\alpha d_s)).$$

Hence, we have

$$f(\alpha )\le f_1(\alpha ):=\frac{1}{2}(\mathrm{\Psi }(v+\alpha d_x)+\mathrm{\Psi }(v+\alpha d_s))-\mathrm{\Psi }(v),$$

which means that $f_1(\alpha )$ gives an upper bound for the decrease of the barrier function $\mathrm{\Psi }(v)$. Furthermore, we can easily verify that $f(0)=f_1(0)=0$.

It follows from (31) that

$$f_1^{\prime }(\alpha )=\frac{1}{2}(tr({\psi }^{\prime }(v+\alpha d_x)\circ d_x)+tr({\psi }^{\prime }(v+\alpha d_s)\circ d_s)).$$

This gives, by (7),

$$f_1^{\prime }(0)=\frac{1}{2}tr({\psi }^{\prime }(v)\circ (d_x+d_s))=-\frac{1}{2}tr({\psi }^{\prime }(v)\circ {\psi }^{\prime }(v))=-\frac{1}{2}{\parallel {\psi }^{\prime }(v)\parallel }_F^2=-2\delta {(v)}^2<0.$$

Let $\eta =v+\alpha d_x$ and $\gamma =v+\alpha d_s$. To simplify the notations we used (and will use below), the following conventions:

$$d_{xi}:={(d_x)}_i,d_{si}:={(d_s)}_i,d_{xij}:={(d_x)}_{ij}\text{and}{d}_{sij}:={(d_s)}_{ij}.$$

(34)

It follows directly from (32) and (33) that

$$\begin{array}{rcl}{f}_1^{″}(\alpha )& \le & \frac{1}{2}\sum _{i=1}^r{\psi }^{″}({\lambda }_i(\eta )){(d_{xi})}^2+\sum _{i<j}{\psi }^{″}({\lambda }_j(\eta ))tr\left({(d_{xij})}^2\right)\\ +\frac{1}{2}\sum _{i=1}^r{\psi }^{″}({\lambda }_i(\gamma )){(d_{si})}^2+\sum _{i<j}{\psi }^{″}({\lambda }_j(\gamma ))tr\left({(d_{sij})}^2\right).\end{array}$$

(35)

In the sequel, we use the short notation $\delta :=\delta (v)$.

Lemma 4.1 One has

$${\parallel d_x\parallel }_F^2+{\parallel d_s\parallel }_F^2\le 4{\delta }^2.$$

Proof Since $\mathcal{Q}$ is a given self-adjoint positive semidefinite linear operator, we have

$$\begin{array}{rcl}〈d_x,d_s〉& =& 〈d_x,\overline{\mathcal{Q}}(d_x)-{\overline{\mathcal{A}}}^T(\mathrm{\Delta }y)〉=〈d_x,\left(P{(w)}^{\frac{1}{2}}\mathcal{Q}P{(w)}^{\frac{1}{2}}\right)(d_x)〉-〈d_x,{\overline{\mathcal{A}}}^T(\mathrm{\Delta }y)〉\\ =& 〈P{(w)}^{\frac{1}{2}}{d}_x,\mathcal{Q}\left(P{(w)}^{\frac{1}{2}}{d}_x\right)〉\ge 0.\end{array}$$

(36)

Hence, we have

$$4{\delta }^2={\parallel d_x+d_s\parallel }_F^2={\parallel d_x\parallel }_F^2+{\parallel d_s\parallel }_F^2+2〈d_x,d_s〉\ge {\parallel d_x\parallel }_F^2+{\parallel d_s\parallel }_F^2.$$

(37)

This completes the proof of the lemma. □

Similar to the proof of Lemma 4.1 in [6], we have the following lemma, which gives an upper bound of $f_1^{″}(\alpha )$ in terms of δ and ${\psi }^{″}(t)$.

Lemma 4.2 One has

$$f_1^{″}(\alpha )\le 2{\delta }^2{\psi }^{″}({\lambda }_{\min }(v)-2\alpha \delta ).$$

Let $\rho (s):[0,\mathrm{\infty })\to (0,1]$ be the inverse function of $-\frac{1}{2}{\psi }^{\prime }(t)$ for $t\le 1$, where $\psi (t)$ is the eligible kernel function. Similar to the LO case, the default step size is chosen. In this paper, we use

$$\tilde{\alpha }:=\frac{1}{{\psi }^{″}(\rho (2\delta ))}$$

(38)

as the default step size.

In what follows, we will show that the barrier function $\mathrm{\Psi }(v)$ in each inner iteration with the default step size $\tilde{\alpha }$, as defined by (38), is decreasing. For this, we need the following technical result.

Lemma 4.3 (Lemma 3.12 in [5])

Let $h(t)$ be a twice differentiable convex function with $h(0)=0$, $h^{\prime }(0)<0$ and let $h(t)$ attain its (global) minimum at $t^{\ast }>0$. If $h^{″}(t)$ is increasing for $t\in [0,t^{\ast }]$, then

$$h(t)\le \frac{t{h}^{\prime }(0)}{2},0\le t\le t^{\ast }.$$

As a consequence of Lemma 4.3 and the fact that $f(\alpha )\le f_1(\alpha )$, which is a twice differentiable convex function with $f_1(0)=0$, and $f_1^{\prime }(0)=-2{\delta }^2<0$, we can easily prove the following lemma.

Lemma 4.4 Let the step size α satisfies the condition $\alpha \le \tilde{\alpha }$. Then

$$f(\alpha )\le -\alpha {\delta }^2.$$

Combining the results of Lemma 4.4 and (38), we have the following theorem, which shows that the default step size (38) yields a sufficient decrease of the barrier function value during each inner iteration.

Theorem 4.5 Let $\tilde{\alpha }$ be the default step size, as given by (38). Then

$$f(\tilde{\alpha })\le -\frac{{\delta }^2}{{\psi }^{″}(\rho (2\delta ))}.$$

(39)

By using the condition (23d), we can conclude that the right-hand side expression in (39) is monotonically decreasing in δ (cf. Lemma 4.7 in [6]). Thus, combining the results of Theorems 4.5 and 3.5, we have

$$f(\tilde{\alpha })\le -\frac{{({\psi }^{\prime }(\varrho (\mathrm{\Psi }(v))))}^2}{4{\psi }^{″}(\rho ({\psi }^{\prime }(\rho (\mathrm{\Psi }(v)))))}.$$

(40)

This expresses the decrease in $\mathrm{\Psi }(v)$ during an inner iteration completely in ψ, its first and second derivatives, and the inverse functions ρ and ϱ.

5 Complexity of the algorithms

In this section, we first derive an upper bound for the number of the iteration bounds by the algorithm depicted in Figure 1. Then we conclude this section by applying the iteration bound to a wide variety of kernel functions.

5.1 Iteration bounds for the algorithms

We need to count how many inner iterations are required to return to the situation where $\mathrm{\Psi }(v)\le \tau$. We use the value of $\mathrm{\Psi }(v)$ after the μ-update by ${\mathrm{\Psi }}_0$, and the subsequent values in the same outer iteration are denoted as ${\mathrm{\Psi }}_k$, $k=1,2,\dots ,K$, where K denotes the total number of inner iterations in the outer iteration.

Let the constants $\beta >0$ and $\gamma \in (0,1]$ be such that for $\mathrm{\Psi }(v)\ge \tau$. We have

$$\frac{{({\psi }^{\prime }(\varrho (\mathrm{\Psi }(v))))}^2}{4{\psi }^{″}(\rho ({\psi }^{\prime }(\varrho (\mathrm{\Psi }(v)))))}\ge \beta \mathrm{\Psi }{(v)}^{1-\gamma }.$$

Note that the left-hand side expression is increasing in $\mathrm{\Psi }(v)$. Therefore, such numbers β and γ certainly exist (take, e.g., $\gamma =1$ and β equals the value of the left-hand side expression for $\mathrm{\Psi }(v)=\tau$). In addition, the appropriate values of β and γ will vary for each eligible kernel function and finding them may not always be straightforward.

The following lemma provides an estimate for the number of inner iterations between two successive barrier parameter updates, in terms of ${\mathrm{\Psi }}_0$ and the parameters β and γ.

Lemma 5.1 One has

$$K\le \frac{{\mathrm{\Psi }}_0^{\gamma }}{\beta \gamma }.$$

Proof The definition of K implies ${\mathrm{\Psi }}_{K-1}>\tau$ and ${\mathrm{\Psi }}_K\le \tau$ and

$${\mathrm{\Psi }}_{k+1}\le {\mathrm{\Psi }}_k-\beta {({\mathrm{\Psi }}_k)}^{1-\gamma },k=0,1,\dots ,K-1.$$

Thus, the conclusion of the lemma follows immediately from Lemma 14 in [5] with $t_k={\mathrm{\Psi }}_k$. This completes the proof of the lemma. □

The number of outer iterations coincides with the number of barrier parameter θ updates until we obtain $r\mu <\epsilon$. It is well known (cf. Lemma Π.17 in [1]) that the number of outer iterations is bounded above by $\frac{1}{\theta }log\frac{r}{\epsilon }$. Thus, an upper bound on the total number of iterations is obtained by multiplying the number of outer iterations and the number of inner iterations.

Theorem 5.2 The total number of iterations is bounded above by

$$\frac{{\mathrm{\Psi }}_0^{\gamma }}{\theta \beta \gamma }log\frac{r}{\epsilon }.$$

From Theorem 5.2 and Corollary 3.4, we obtain the iteration bound for the algorithm depicted in Figure 1 follows:

$$\frac{1}{\theta \beta \gamma }{\left(r\psi \left(\frac{\varrho (\frac{\tau }{r})}{\sqrt{1-\theta }}\right)\right)}^{\gamma }log\frac{r}{\epsilon },$$

(41)

which means that the total number of iterations is completely determined by the parameters θ, β, γ, τ, and the eligible kernel function $\psi (t)$.

5.2 Application to the eligible kernel functions

It follows from Theorem 5.2 that the iteration bound of the algorithms depends on the parameters β and γ and the upper bound on ${\mathrm{\Psi }}_0$. Since these are different for different eligible kernel functions, the iteration bounds will also vary. Similar to the analysis considered in [6], Section 6.1, for the LO case, the iteration bounds for large- and small-update methods based on the eligible kernel functions can be performed in a systematic way by using the following scheme.

• Step 0: Input an eligible kernel function $\psi (t)$; an update parameter θ, $0<\theta <1$; a threshold parameter τ; and an accuracy parameter ε.

• Step 1: Solve the equation $-\frac{1}{2}{\psi }^{\prime }(t)=s$ to get $\rho (s)$, the inverse function of $-\frac{1}{2}{\psi }^{\prime }(t)$, $t\in (0,1]$. If the equation is hard to solve, derive a lower bound for $\rho (s)$.

• Step 2: Calculate the decrease of $\mathrm{\Psi }(v)$ in terms of δ for the default step size $\tilde{\alpha }$ from

  $$f(\tilde{\alpha })\le -\frac{{\delta }^2}{{\psi }^{″}(\rho (2\delta ))}.$$

• Step 3: Solve the equation $\psi (t)=s$ to get $\varrho (s)$, the inverse function of $\psi (t)$, $t\ge 1$. If the equation is hard to solve, derive the lower and upper bounds for $\varrho (s)$.

• Step 4: Derive a lower bound for $\delta (v)$ in terms of $\mathrm{\Psi }(v)$ by using

  $$\delta (v)\ge \frac{1}{2}{\psi }^{\prime }\left(\varrho (\mathrm{\Psi }(v))\right).$$

• Step 5: Using the results of Step 3 and Step 4 find positive constants β and γ, with $\gamma \in (0,1]$, such that

  $$f(\tilde{\alpha })\le -\beta \mathrm{\Psi }{(v)}^{1-\gamma }.$$

• Step 6: Calculate the uniform upper bound ${\mathrm{\Psi }}_0$ for $\mathrm{\Psi }(v)$ from

  $${\mathrm{\Psi }}_0\le L_{\psi }(r,\theta ,\tau ):=r\psi \left(\frac{\varrho (\frac{\tau }{r})}{\sqrt{1-\theta }}\right).$$

• Step 7: Derive an upper bound for the total number of iterations from

  $$\frac{{\mathrm{\Psi }}_0^{\gamma }}{\theta \beta \gamma }log\frac{r}{\epsilon }.$$

• Step 8: Set $\tau =O(r)$ and $\theta =\mathrm{\Theta }(1)$ so as to calculate an iteration bound for large-update methods, or set $\tau =O(1)$ and $\theta =\mathrm{\Theta }(\frac{1}{\sqrt{r}})$ to get an iteration bound for small-update methods.

The resulting iteration bounds for a wide class of eligible kernel functions have been outlined in a series of papers [5–10, 14, 15, 17–19, 21] starting with [6] for LO, we immediately get the iteration bounds for large- and small-update methods for CQSCO. The resulting iteration bounds are summarized in the 3rd and 4th columns of Table 1. For the detailed analysis of the algorithms can be refereed to the given references.

Table 1 Complexity results for the eligible kernel functions

Remark 5.3 For large-update methods, the currently best-known iteration bound is

$$O(\sqrt{r}logrlog\frac{r}{\epsilon }).$$

In particular, for ${\psi }_3(t)$ and ${\psi }_4(t)$ this bound is obtained if we choose $q=\frac{1}{2}logr$, and for ${\psi }_7(t)$ and ${\psi }_8(t)$ this bound is obtained if we choose $q=logr$. The same bound is achieved for ${\psi }_{18}(t)$, also by taking $p=1$ and $q=\frac{1}{2}logr$.

Remark 5.4 For small-update methods, the currently best-known iteration bound is

$$O(\sqrt{r}log\frac{r}{\epsilon }).$$

In particular, for ${\psi }_3(t)$, ${\psi }_4(t)$, ${\psi }_7(t)$, ${\psi }_8(t)$, ${\psi }_{16}(t)$ and ${\psi }_{18}(t)$, this bound is derived is we take $q=O(1)$.

Both for large- and small-update methods, the order of the iteration bounds are obtained as good as the bounds for the LO case except that n is replaced by r, the rank of the EJA. Thus, the iteration bounds are as good as they can be in the current state-of-the-art.

6 Numerical results

In this section, we report the computational performance of the algorithm depicted in Figure 1 for CQSDO, which is an important cases of CQSCO. The numerical experiments are carried out on a PC with Intel (R) Core (TM) i5-2500 Duo CPU at 3.30 GHz and 8 GB of physical memory. The PC runs MATLAB Version: 7.11.0.584 (R2010b) on a Windows 7 Enterprise 64-bit operating system.

We consider the primal problem of CQSDO in the standard form

$$min\{\frac{1}{2}X•\mathcal{Q}(X)+C•X:A_i•X=b_i,i=1,2,\dots ,m,X⪰0\},$$

and its dual problem

$$max\{-\frac{1}{2}X•\mathcal{Q}(X)+b^{T}y:\sum _{i=1}^m{y}_i{A}_i-\mathcal{Q}(X)+S=C,S⪰0\}.$$

Here, $\mathcal{Q}:{\mathbf{S}}^n\to {\mathbf{S}}^n$ is a given self-adjoint positive semidefinite linear operator on ${\mathbf{S}}^n$, i.e., for any $A,B\in {\mathbf{S}}^n$, $\mathcal{Q}(A)•B=A•\mathcal{Q}(B)$ and $\mathcal{Q}(A)•A\ge 0$. $b\in {\mathbf{R}}^m$ is a given vector, $C\in {\mathbf{R}}^{n\times n}$ is a given matrix. Without loss of generality we assume that the matrices $A_i$, $i=1,2,\dots ,m$ are linearly independent and CQSDO satisfy the IPC. The detailed discussion and analysis of primal-dual IPMs for CQSDO can be found in [24, 28, 34].

Let us define

$$P:=X^{\frac{1}{2}}{\left(X^{\frac{1}{2}}S{X}^{\frac{1}{2}}\right)}^{-\frac{1}{2}}{X}^{\frac{1}{2}}[=S^{-\frac{1}{2}}{\left(S^{\frac{1}{2}}X{S}^{\frac{1}{2}}\right)}^{\frac{1}{2}}{S}^{-\frac{1}{2}}],$$

(42)

and also we define $D:=P^{\frac{1}{2}}$. This leads to the definition of the following variance matrix:

$$V:=\frac{1}{\sqrt{\mu }}{D}^{-1}X{D}^{-1}[=\frac{1}{\sqrt{\mu }}DSD].$$

(43)

Furthermore, we define the scaled search directions as follows:

$$D_X:=\frac{1}{\sqrt{\mu }}{D}^{-1}\mathrm{\Delta }X{D}^{-1}\text{and}{D}_S:=\frac{1}{\sqrt{\mu }}D\mathrm{\Delta }SD.$$

(44)

The scaled search direction $(D_X,\mathrm{\Delta }y,D_S)$ is computed through solving the following linear system:

$$\begin{array}{c}{\overline{A}}_i•D_X=0,i=1,2,\dots ,m,\hfill \\ \sum _{i=1}^m\mathrm{\Delta }{y}_i{\overline{A}}_i-\overline{\mathcal{Q}}(D_X)+D_S=0,\hfill \\ D_X+D_S=V^{-1}-V,\hfill \end{array}$$

(45)

where

$${\overline{A}}_i:=\frac{1}{\sqrt{\mu }}D{A}_{i}D,i=1,2,\dots ,m\text{and}\overline{\mathcal{Q}}(D_X):=D\mathcal{Q}(D{D}_{X}D)D.$$

Then the new search direction $(\mathrm{\Delta }X,\mathrm{\Delta }y,\mathrm{\Delta }S)$ is obtained from (44). If $(X,y,S)\ne (X(\mu ),y(\mu ),S(\mu ))$, then $(\mathrm{\Delta }X,\mathrm{\Delta }y,\mathrm{\Delta }S)$ is nonzero. The new iterate is obtained by taking a default step size α along the search directions as follows:

$$X_{+}:=X+\alpha \mathrm{\Delta }X,y_{+}:=y+\alpha \mathrm{\Delta }y\text{and}{S}_{+}:=S+\alpha \mathrm{\Delta }S.$$

(46)

It should be noted that the default step size (38) selected during each inner iteration is small enough for analyzing the algorithm, while in practice it should be chosen large enough for the efficiency of the algorithm. In the following test problem, we choose the maximum allowed step size such that the next iterate satisfying the positive semidefiniteness condition, i.e., $X+\alpha D_X⪰0$ and $S+\alpha D_S⪰0$.

We consider the following special CQSDO example with $\mathcal{Q}(X)=E$:

$$\begin{array}{c}{A}_1=\left(\begin{array}{cccccccc}2& 1& 0& 1& 1& 1& 2& 1\\ 1& 0& 0& 2& 1& 2& 1& 0\\ 0& 0& 2& 1& 0& 2& 1& 0\\ 1& 2& 1& 0& 1& 1& 2& 2\\ 1& 1& 0& 1& 2& 0& 1& 1\\ 1& 2& 2& 1& 0& 0& 2& 1\\ 2& 1& 1& 2& 1& 2& 2& 2\\ 1& 0& 0& 2& 1& 1& 2& 0\end{array}\right),A_2=\left(\begin{array}{cccccccc}2& 1& 0& 1& 1& 1& 2& 1\\ 1& 0& 0& 2& 1& 2& 1& 0\\ 0& 0& 2& 1& 0& 2& 1& 0\\ 1& 2& 1& 0& 1& 1& 2& 2\\ 1& 1& 0& 1& 2& 0& 1& 1\\ 1& 2& 2& 1& 0& 0& 2& 1\\ 2& 1& 1& 2& 1& 2& 2& 2\\ 1& 0& 0& 2& 1& 1& 2& 0\end{array}\right),\hfill \\ A_3=\left(\begin{array}{cccccccc}2& 1& 0& 1& 1& 1& 2& 1\\ 1& 0& 0& 2& 1& 2& 1& 0\\ 0& 0& 2& 1& 0& 2& 1& 0\\ 1& 2& 1& 0& 1& 1& 2& 2\\ 1& 1& 0& 1& 2& 0& 1& 1\\ 1& 2& 2& 1& 0& 0& 2& 1\\ 2& 1& 1& 2& 1& 2& 2& 2\\ 1& 0& 0& 2& 1& 1& 2& 0\end{array}\right),A_4=\left(\begin{array}{cccccccc}2& 1& 0& 1& 1& 1& 2& 1\\ 1& 0& 0& 2& 1& 2& 1& 0\\ 0& 0& 2& 1& 0& 2& 1& 0\\ 1& 2& 1& 0& 1& 1& 2& 2\\ 1& 1& 0& 1& 2& 0& 1& 1\\ 1& 2& 2& 1& 0& 0& 2& 1\\ 2& 1& 1& 2& 1& 2& 2& 2\\ 1& 0& 0& 2& 1& 1& 2& 0\end{array}\right),\hfill \\ b=\left(\begin{array}{c}8\\ 4\\ 12\\ 8\end{array}\right),\mathcal{Q}(X)=\left(\begin{array}{cccccccc}1& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right),\hfill \\ C=\left(\begin{array}{cccccccc}6& 4& 4& 6& 5& 1& 5& 3\\ 4& 4& 6& 5& 6& 4& 3& 4\\ 4& 6& 6& 3& 4& 6& 5& 5\\ 6& 5& 3& 0& 5& 5& 3& 4\\ 5& 6& 4& 5& 4& 2& 3& 2\\ 1& 4& 6& 5& 2& 2& 6& 4\\ 5& 3& 5& 3& 3& 6& 8& 6\\ 3& 4& 5& 4& 2& 4& 6& 2\end{array}\right).\hfill \end{array}$$

In the test problems, we use the threshold parameter $\tau =3$, the accuracy parameter $\epsilon ={10}^{-7}$, and the update parameter $\theta =\frac{1}{2\sqrt{n}}$ with $n=8$ in the implementation. In this case, the algorithm depicted in Figure 1 is indeed a small-update method. We choose $X=E$, $y=e$ and $S=E$ as the starting point for our algorithm. Here E and e denote the identity matrix of dimension 8 and the identity vector of dimension 4, respectively. Note that the point is strictly feasible. The initial value of the barrier parameter μ is $X•S/n$ with $n=8$, i.e., 1. We can easily check that $\mathrm{\Psi }(X,S;\mu )=0<\tau =3$. So these data can, indeed, be used to initialize our algorithm.

An optimal solution of the primal problem is given by

$$X^{\ast }=\left(\begin{array}{cccccccc}0.1788& 0.1173& 0.0605& 0.0769& -0.0034& -0.0088& 0.1511& -0.1572\\ 0.1173& 0.3874& 0.0684& 0.1492& 0.1188& 0.3490& 0.2450& 0.1194\\ 0.0605& 0.0684& 0.5162& 0.0046& 0.1338& -0.0294& 0.1971& -0.0010\\ 0.0769& 0.1492& 0.0046& 0.0789& 0.0746& 0.1092& 0.0864& -0.0252\\ -0.0034& 0.1188& 0.1338& 0.0746& 0.2652& 0.1015& 0.0095& -0.0208\\ -0.0088& 0.3490& -0.0294& 0.1092& 0.1015& 0.4175& 0.1547& 0.2746\\ 0.1511& 0.2450& 0.1971& 0.0864& 0.0095& 0.1547& 0.2667& 0.0321\\ -0.1572& 0.1194& -0.0010& -0.0252& -0.0208& 0.2746& 0.0321& 0.3780\end{array}\right)$$

and for the dual problem an optimal solution is given by

$$\begin{array}{c}{y}^{\ast }=(1.0564;0.8583;1.1417;1.0416),\hfill \\ S^{\ast }=\left(\begin{array}{cccccccc}0.0660& 0.0192& -0.0228& -0.0212& 0.0403& -0.0652& -0.0033& 0.0698\\ 0.0192& 0.0206& -0.0149& -0.0053& 0.0207& -0.0473& 0.0052& 0.0361\\ -0.0228& -0.0149& 0.0366& 0.0483& -0.0496& 0.0579& -0.0427& -0.0427\\ -0.0212& -0.0053& 0.0483& 0.0789& -0.0652& 0.0527& -0.0681& -0.0379\\ 0.0403& 0.0207& -0.0496& -0.0652& 0.0691& -0.0819& 0.0532& 0.0645\\ -0.0652& -0.0473& 0.0579& 0.0527& -0.0819& 0.1340& -0.0414& -0.1069\\ -0.0033& 0.0052& -0.0427& -0.0681& 0.0532& -0.0414& 0.0705& 0.0194\\ 0.0698& 0.0361& -0.0427& -0.0379& 0.0645& -0.1069& 0.0194& 0.0945\end{array}\right).\hfill \end{array}$$

The respective objective values are $\frac{1}{2}X•\mathcal{Q}(X)+C•X=32.959138158$ and $-\frac{1}{2}X•\mathcal{Q}(X)+b^{T}y=32.959138116$, and the duality gap $X•S$ is $4.2163\times {10}^{-8}$, which is less than 10⁻⁷.

The numerical results of IPM for the sample problem of CQSDO based on ${\psi }_1(t)$ with $\theta =\frac{1}{2\sqrt{n}}$ are summarized in Table 2. For our small-update method, we need 19 main iterations to reach our accuracy. To save space, we show the primal and dual objective value at the moments when the duality gap is reduced again with a factor 10, until the desired accuracy is achieved. The numerical results are summarized in Table 2.

Table 2 Output of IPM for the sample problem of CQSDO based on ${\mathit{\psi }}_{\mathbf{1}}\mathbf{(}\mathit{t}\mathbf{)}$ with $\mathit{\theta }\mathbf{=}\frac{\mathbf{1}}{\mathbf{2}\sqrt{\mathit{n}}}$

It is clear from Table 2 that the small-update method presented in this paper is not efficient from a practical point of view, just as the feasible IPMs with the best theoretical performance are far from practical. In fact, our algorithm suffers from the usual drawback of primal-dual IPMs that the number of iterations needed for convergence leads to be close to the upper bound, namely, $O(\sqrt{n}log\frac{n}{\epsilon })$. This is due to the small, fixed μ-updates (i.e., ${\mu }_{+}=(1-\theta )\mu$ with $\theta =\frac{1}{2\sqrt{n}}$ for CQSDO). It is desirable to make the largest possible update θ at each iteration, albeit at the cost of extra computation.

In order to reveal the impact of the update parameter θ on the performance of the algorithm, we take the larger possible update parameter $\theta =0.9$ in the implementation. In this case, the algorithm depicted in Figure 1 is indeed a large-update method. We only need 14 main iterations to reach our accuracy. The outputs of IPMs for the sample problem of CQSDO based on ${\psi }_1(t)$ with $\theta =0.9$ are shown in Table 3.

Table 3 Output of IPM for the sample problem of CQSDO based on ${\mathit{\psi }}_{\mathbf{1}}\mathbf{(}\mathit{t}\mathbf{)}$ with $\mathit{\theta }\mathbf{=}\mathbf{0.9}$

It is clear from Table 3 that the iteration number of the algorithm depend on the update parameter θ. A larger value of the update parameter θ gives rise to better results. However, it should be pointed out that the update parameter θ would be too large to solve the problem in the computational procedure. In the solution procedure, we might use the dynamic updates of the barrier parameter, as described in [1]. This may significantly enhance the practical performance of the proposed algorithm.

7 Conclusions and remarks

In this paper, we presented a unified approach and comprehensive treatment of primal-dual IPMs for CQSCO based on the entire class of the eligible kernel functions. For large-update methods the best iteration bound is $O(\sqrt{r}logrlog\frac{r}{\epsilon })$ and for small-update methods all iteration bounds have the same order of magnitude, namely, $O(\sqrt{r}log\frac{r}{\epsilon })$, which almost closes the gap between the iteration bounds for large- and small-update methods. Some preliminary numerical results are provided to demonstrate the computational performance of the algorithm depicted in Figure 1.

The paper generalizes results obtained in the following papers, [6] where Bai et al. consider kernel-based primal-dual IPMs for LO, and [11, 16, 30] and [23] where Bai et al., El Ghami et al., Wang et al. and Vieira consider the same type of IPMs for SOCO, SDO, CQSDO and SCO, respectively. It turns out that the iterations bounds are the same as for the non-negative orthant except that n is replaced by r, the rank of the EJA. However, the analysis of the proposed algorithm is far more complicated in [6, 11, 16, 23]. This is due to the following fact that we lose the orthogonality of search directions that exist in LO, SOCO, SDO, and SCO cases does not hold for CQSCO.

Some interesting topics for further research remain. First, the search direction used in this paper is based on the NT-symmetrization scheme and it is natural to ask if other symmetrization schemes can be used. Second, although we present a simple examples to show the computational performance of the proposed algorithm, more numerical experiments are desired to compare the behavior of our algorithm with other existing IPMs. Finally, the extension to general nonlinear optimization over symmetric cone deserves to be investigated.

Endnote

^a It may be worth mentioning that if we use the kernel function of the classical logarithmic barrier function, i.e., $\psi (t)=\frac{1}{2}(t^2-1)-logt$, then ${\psi }^{\prime }(t)=t-t^{-1}$, whence $-{\psi }^{\prime }(v)=v^{-1}-v$, and hence system (7) then coincides with the classical system (6).