Primal-dual interior-point algorithms for convex quadratic circular cone optimization

In this paper we focus on a class of special nonsymmetric cone optimization problem called circular cone optimization problem, which has a convex quadratic function as the objective function and an intersection of a non-self-dual circular cone and linear equations as the constraint condition. Firstly we establish the algebraic relationships between the circular cone and the second-order cone and translate the original problem from the circular cone optimization problem to the second-order cone optimization problem. 
Then we present kernel-function based primal-dual interior-point algorithms for solving this special circular cone optimization 
and derive the iteration bounds for large- and small-update methods. 
Finally, some preliminary numerical results are provided to demonstrate the computational performance of the proposed algorithms.


1.
Introduction. The symmetric conic optimization (SCO) is a wide class of problems that contains linear optimization (LO), second-order cone optimization (SOCO) and semidefinite optimization (SDO) as special cases. There are many solution approaches for SCO. Among them, the interior-point methods (IPMs) gain much more attention. Faybusovich [8] made the first attempt to generalize IPMs to SCO by using Euclidean Jordan algebras. Since then, many researchers have studied IPMs for SCO and achieved plentiful and beautiful results [3,8,9,14]. For an overview of these and related results, we refer to the recent monograph [2] and the references within.
Nonsymmetric conic optimization (NSCO) is a kind of optimization model which has a linear objective function and constraint conditions with a non-self-dual cone and an affine space. It is applied to many real-life problems [19] and theoretical problems [20]. NSCO has the characteristic of SCO in form, but because of the non-self-dual of the constraint cone in essence, it is difficult to study the theory and design the algorithm. Recently, much attention has been paid to NSCO. For example, Nesterov [11] addressed a new IPM, which is based on an extension of the idea of self-scaled optimization to the NSCO case. Motivated by the work of Nesterov, Skajaa et al. [15] presented a homogeneous interior-point algorithm for nonsymmetric convex conic optimization. Matsukawa and Yoshise [10] proposed a primal barrier function Phase I algorithm for NSCO. Chen et al. [7,22,23,24] considered a special nonsymmetric cone, the circular cone. It is a pointed closed convex cone having hyperspherical sections orthogonal to its axis of revolution about which the cone is invariant to rotation. In addition, the second-order cone is a special case of the circular cone when the rotation angle is 45 degrees. They explored some important features about circular cone, such as characterizing its tangent cone, normal cone, and second-order regularity, etc.
In this paper, we consider convex quadratic optimization over a circular cone (CQO-CC) with linear equality constraints as follows where Q ∈ S n ++ , A ∈ R m×n with rank(A) = m, c ∈ R n , b ∈ R m and C n θ is an n-dimensional circular cone. The CQO-CC is actually LO over a circular cone (LO-CC) when Q is a null matrix, i.e., CQO-CC is a generalization of LO-CC. And it is also a generalization of convex quadratic optimization over a second-order cone (CQO-SOC). Both of the perspectives give us great ideas about the following studies.
The main reason why we choose to study this problem is that the CQO-CC problem is applied to many real-life engineering problems [19,25]. And some nonconvex optimization, combination optimization and NP-hard problems can be translated to this problem. Moreover, it plays a crucial role in the force optimization problem (FOP) [6]. The FOP comes up in several applications and setting, such as grasp optimization, or real-time grasp control and force optimization for the legs of a quadruped robot [6].
Because of the self-duality of the second-order cone, the CQO-SOC is a special case of convex quadratic SDO (CQSDO) with many kinds of polynomial interiorpoint algorithms [17,18]. The circular cone optimization (CCO), however, is NSCO and unfortunately we still can not find a good polynomial interior-point algorithm that applies to it. So the basic idea in this paper is to establish the algebraic relationships between the circular cone and the second-order cone and then convert the CCO problems into the SOCO problems by the algebraic relationships.
The purpose of the paper is to present primal-dual IPMs for CQO-CC based on the kernel function of the classic barrier function, i.e., which is a special case of the eligible kernel function [4]. We first translate the CQO-CC problem as the CQO-SOC problem by using the algebraic relationships between the circular cone and the second-order cone. Then we propose kernelfunction based primal-dual IPMs for CQO-CC and establish the iteration bounds for large-and small-update methods. Some preliminary numerical results are provided to demonstrate the computational performance of the proposed algorithms. The paper is organized as follows. In Section 2, the algebraic relationships between the circular cone and the second-order cone are developed. In Section 3, we present the framework of kernel-function based primal-dual IPMs for CQO-CC. The analysis and complexity of the algorithms for large-and small-update methods are presented in Section 4. In Section 5, we report some preliminary numerical results. Finally, some conclusions and remarks are made in Section 6.
Some of the notations used throughout the paper are as follows: || · || denotes the Euclidean norm, i.e., 2-norm. R n , R n + and R n ++ denote the set of vectors with n components, the set of nonnegative vectors and the set of positive vectors, respectively. R m×n denotes the m × n real matrices. S n ++ denotes symmetric positive definite n × n matrices. We use rank(A) to denote the rank of the matrix A. Finally, if g(x) ≥ 0 is a real valued function of a real nonnegative variable, the notation g(x) =O(x) means that g(x) ≤cx for some positive constantc and g(x) = Θ(x) that c 1 x ≤ g(x) ≤ c 2 x for two positive constants c 1 and c 2 .
2. The relationships between the circular cone and the second-order cone.
2.1. The second-order cone. The second-order (also called Lorentz or ice−cream) cone in R n is given by where n (n ≥ 2) is some natural number and x 2:n denotes (x 2 , x 3 , . . . , x n ) T . For any two vectors x, s ∈ R n , the bilinear operator • is defined by This operator is commutative and, moreover, (R n , •) is a Jordan algebra. Note that the map s → x • s is linear. The matrix of this linear map, with respect to the standard basis, is denoted as mat(x), and one may easily verify that it is an arrow-shaped (symmetric) matrix where I n−1 denotes the n − 1-dimensional identity matrix. Some properties of the arrow-shaped matrix mat(x) are as follows [3] .
, the interior of K n ), then mat(x −1 )e = x −1 . The maximal and minimal eigenvalues of mat(x) are denoted as λ max (x) and λ min (x), respectively. These are given by One may easily verify that the vectors are the eigenvectors for the eigenvalues λ max (x) and λ min (x), respectively. The so-called spectral decomposition of a vector x ∈ R n is given by The importance of the decomposition (7) is that it enables us to extend the definition of any function ψ(t) to a function that maps K into K. Let x ∈ R n with the spectral decomposition as defined by (7). The vector-valued function ψ(x) is given by Let ψ(t) be differentiable. Replacing ψ(λ max (x)) and ψ(λ min (x)) in (8) by ψ (λ max (x)) and ψ (λ min (x)), respectively, we can conclude that the vector-valued function ψ (x) = ψ (λ max (x)) z 1 + ψ (λ min (x)) z 2 (9) is defined as well.
In addition, the definition of the trace function for x ∈ R n is given by Then for any x, s ∈ R n , one has Tr(x • s) = 2x T s.

2.2.
The circular cone. The circular cone is a pointed closed convex cone having hyperspherical sections orthogonal to its axis of revolution about which the cone is invariant to rotation. Let its half-aperture angle be θ where θ ∈ (0, π 2 ), see Figure  1. Then, we denote the n-dimensional circular cone by C n θ which is expressed as It should be pointed out that the second-order cone is actually a special case of the circular cone when the rotation angle is π 4 . In fact, we have  i.e., C n π 4 = K n . This is illustrated in figure 2.
Definition 2.1. Let K be a cone. The set is called the dual cone of K.
Definition 2.2. If K is a cone satisfied with i.e., K * = K, where K * is given by Definition 2.1, then it is called self-dual.
It is well known the second-order cone K n is self-dual, i.e., (K n ) * = K n . However, the circular cone C n θ is not. In fact, the dual cone of the circular cone C n θ is given by . , x n ) T ∈ R n | sin θ ||x|| ≤ x 1 }. This is exactly the hard point in solving the CQO-CC problems.
2.3. The algebraic relationships. From (11), we have which yields

YANQIN BAI, XUERUI GAO AND GUOQIANG WANG
For simplicity, let us denote It is clear that the matrix B is positive definite and its inverse matrix is Then, the expression (13) is equivalent to The following theorem implies the relationships among C n θ , (C n θ ) * and K n . Theorem 2.3 (Throrem 2.1 in [23]). Let K n and C n θ be defined as in (2) and (11), respectively. Thus, C n π 2 −θ is expressed as: According to Theorem 2.3 we have the following simplified but important relational expressions: So far, we have deduced the relationships which will help us to translate the problems we study to the SOCO problems with better properties.
3.1. The CQO-SOC problem. The primal problem of CQO-SOC in standard form is where Q ∈ S n ++ , A ∈ R m×n with rank(A) = m, c ∈ R n and b ∈ R m . We notice that x ∈ K n is a geometric constraint, and it is difficult to operate. So we rewrite it in the form which is an algebraic constraint. Therefore the problem (18) can be rewritten as follows The Lagrangian function of the problem (19) is given by where λ ∈ R + , ν ∈ R m . Let Then s ∈ K n . Hence Thus the above Lagrangian function (20) can be simplified as It is clear that the Lagrangian function is convex with respect to x, therefore if its partial derivative with respect to x is equal to 0, i.e., then L(x, s, ν) will achieve a minimum And the dual function of (23) is given by Hence, the Lagrangian dual problem of CQO-SOC is expressed as where y = −ν ∈ R m , x ∈ R n . Solving the pair of problems (18) and (27) is equivalent to solving the optimality conditions, i.e., KKT conditions (22), and then we have an equivalence transformation of the KKT conditions s T x = 0.
Since the CQO-SCO problem is a special case of the CQSDO problem, we can solve it by using kernel-function based primal-dual IPMs for CQSDO (e.g., [17,18]).

3.2.
The CQO-CC problem. The geometrical constraint x ∈ C n θ can be rewritten as Then the equivalent form of (P) is given by Hence, its Lagrangian function is Lemma 3.1. For x, λ given by (P ) and (30), respectively. Let Proof. For given x ∈ C n θ , first we have Multiplying both sides of this inequality by λ 2 sin 2 θ cos 2 θ, one has Moving cos 2 θ λ 2 sin 4 θ x 2 1 to the right hand side of the inequality, we obtain Plugging λ 2 cos 4 θ into n i=2 x 2 i , the equivalence of the inequality is as follows This implies that s ∈ (C n θ ) * . Furthermore, we can easily verify that This completes the proof.
From Lemma 3.1, we can conclude that (30) is equivalent to The dual function corresponding to (31) can be written as Then we have the dual problem of CQO-CC as follows Similar to the CQO-SOC problem (e.g., [1,5,12,16]), finding optimal solution of the pair of (P) and (D) is equivalent to solving the following KKT conditions s T x = 0.
Because C n θ and (C n θ ) * in the above conditions are not the same cone, it is difficult to apply the primal-dual interior-point algorithm in (32). Whereupon according to the relationships (17) we have established in Section 2, the KKT conditions (32) can be rewritten as s T x = 0.
Further more (P) and (D) are transformed as follows Thus we can apply kernel-function based primal-dual IPMs to solving CQO-CC.
3.3. The central path. Throughout the paper, we assume that both (P 0 ) and (D 0 ) satisfy the interior-point condition (IPC), i.e., there exists (x 0 , y 0 , s 0 ) such that It is well known that the IPC implies that (P 0 ) and (D 0 ) have optimal solutions with duality gap zero [1]. Under the IPC assumption, the system (33) is equivalent to the following system (see [1,12]) mat(x)s = 0.
The basic idea of primal-dual IPMs is to replace the third equation in (34) by the parameterized equation mat(x)s = µe with µ > 0. Thus we consider the following system mat(x)s = µe.
For each µ > 0, the parameterized system (35) has a unique solution (x(µ), y(µ), s(µ)) and we call x(µ) the µ-center of (P 0 ) and (y(µ), s(µ)) the µ-center of (D 0 ). The set of µ-centers (with µ running through all the positive real numbers) gives a homotopy path, which is called the central path of (P 0 ) and (D 0 ). Note that at the µ-center we have Then we can derive the duality gap as follows If µ → 0, then the limit of the central path exists and since the limit points satisfy the complementarity condition mat(x)s = 0, the limit yields an ε-approximate solution for (P 0 ) and (D 0 ).

3.4.
The new search directions for CQO-CC. IPMs follow the central path approximately and approach the optimal set of CQO-CC by letting µ go to zero. The search direction is usually derived from a certain Newton-type system. First, we should point out that a straightforward approach to obtain such a system fails to define unique search directions. By linearizing the system (35), we have mat(s)∆x + mat(x)∆s = µe − mat(x)s.
The system (38) has a unique solution if and only if the coefficient matrix of (38) (or that of the system equivalent to (38)) is an invertible matrix. It requires that the matrix AB −1 mat(x)A T should be a positive definite matrix. However, considering that mat(x) may not be positive definite, the matrix AB −1 mat(x)A T may not either. Thus, we can not guarantee that the system (38) has a unique solution.
Fortunately, it is well known that this difficulty can be solved by applying a scaling scheme. This goes as follows.
In this paper, we consider the NT-symmetrization scheme (e.g., [12,21]). Let u := − 1 2 , where is the NT-scaling point of x and s . Let us define andĀ Plugging (39) and (40) Obviously, the matrix D is invertible. Thus the corresponding system has a unique solution and so is (41). It may be worth mentioning that if we use the kernel function of the classical logarithmic barrier function given by (1) Thus the system (41) can be rewritten as Since the system (43) has the same coefficient matrix as the system (41), the system (43) has a unique solution.

YANQIN BAI, XUERUI GAO AND GUOQIANG WANG
Before ending this section, we define the barrier function as follows which is in terms of the original variables x and s, and the barrier parameter µ. It is easy to verify that Ψ(v) ≥ 0 and Hence, the value of Ψ(v) can be considered as a measure for the distance between the given iterate (x, y, s) and the µ-center (x(µ), y(µ), s(µ)).

The generic kernel-function based primal-dual IPMs for CQO-CC.
As we mentioned in Section 3.4., the closeness between (x, y, s) and (x(µ), y(µ), s(µ)) is measured by the value of Ψ(v), with τ ≥ 1 as a threshold value. If Ψ(v) ≤ τ , then we start a new outer iteration by performing a µ-update, otherwise we enter an inner iteration by computing the search directions at the current iterates with respect to the current value of µ and use (44) to get new iterates. If necessary, we will repeat the procedure until we find the iterate that is in the neighborhood of (x(µ), y(µ), s(µ)). Then µ is reduced by the factor 1 − ω with 0 < ω < 1 and we apply Newton's method targeting at the new µ-centers, and so on. Note that at the µ-center the duality gap equals µ, according to (37). This process is repeated until µ is small enough, i.e., 2µ < ε. At this stage we have found an ε-approximate solution of (P 0 ) and (D 0 ). The above discussion can be summarized in the form of the generic kernelfunction based primal-dual IPMs for CQO-CC presented in Figure 3.
The following lemma implies that any eligible kernel function ψ(t) is exponential convexity, which has been proven to be very useful in the analysis of interior-point algorithms based on the kernel functions [4,12,17].
The kernel function ψ(t) we selected is actually a special case of that in [17] when p = 1 and q = 1. It is therefore stated without proof.
As a consequence of Lemma 4.1, we have the following lemma, which is crucial for the analysis of the interior-point algorithms for CQO-CC.
By Lemma 4.2, we have Denote the decrease in Ψ(v) during an inner iteration as Then we have f (α) ≤ f 1 (α), where One can easily verify that f (0) = f 1 (0) = 0. The first and second derivatives of f 1 (α) are as follows and In the analysis of the algorithms, we need a measure for the distance between the given pair (x, y, s) and the µ-center (x(µ), y(µ), s(µ)). For this we define a norm-based proximity measure δ(x, s; µ) as follows Using (49), the third equation in the system (43), and (51), we have Below we use the short notation δ := δ(v), and we have the main important result (see, Lemma 3.3 in [5]).

4.2.
Choice of a suitable step size α. Now we choose the suitable step size for the algorithms, which makes x + and s + be feasible and Ψ(v + ) − Ψ(v) decrease sufficiently. The following way for selecting the step size is an extension of that in the LO case. Thus for the proof of some following lemmas we can refer to [3].
In what follows, we use the notatioñ and we will useα as the default step size. Obviously, we haveᾱ ≥α.
Lemma 4.7. Let the step size α be such that α ≤ᾱ. Then Combining Lemma 4.6 and Lemma 4.7, we can derive the following theorem.
Theorem 4.8. Withα being the default step size, as given by (56), one has .
It should be noted that the right-hand side expression in (58) is monotonically decreasing in δ [4].
In the sequel, we want to express the decrease of f (α) as a function of Ψ(v). To this end, we need a lower bound on δ(v) in terms of Ψ(v). Such a bound can be obtained from the following theorem and lemmas.
It is easy to find the function t − 1 t is monotonically increasing for t ≥ 1. Replacing (Ψ(v)) in the above expression by (1 + 2Ψ(v)) 1 2 , we have This completes the proof.
Consider the kernel function we choose is: The inverse function t = ρ(u) follows by solving t from the equation − 1 2 ψ (t) = u : This enable us to compute ρ(u) exactly: Combining Theorem 4.8, (62) and Lemma 4.11, we have 4.3. Iteration bounds. It should be mentioned that during the course of the algorithms the largest value of Ψ(v) occur just after the update of µ. So next we derive an estimate for the effect of a µ-update on the value of Ψ(v). The following lemma is a reformulation of Theorem 3.2 in [4].
Let v ∈ K n + and β ≥ 1. Then One may easily verify that ψ(t) = t 2 −1 2 − log t ≤ t 2 2 for t ≥ 1. Along with the right-hand inequality of (60), we have

YANQIN BAI, XUERUI GAO AND GUOQIANG WANG
The expression (64) provides a uniform upper bound for Ψ(v) during the execution of the algorithms. In the sequel, we use the symbol L to denote this uniform upper bound of Ψ(v).
We need to count how many inner iterations are required after a µ-update to return to the situation where Ψ(v) ≤ τ . We denote the value of Ψ(v) after the µ-update as Ψ 0 and the subsequent values in the same outer iteration as Ψ k , k = 1, 2, . . . , K, where K denotes the number of inner iterations between two successive µ-updates. Obviously, Ψ 0 ≤ L.
This lemma provides an estimate for the number of inner iterations between two successive barrier parameter updates, in terms of L and constants η and γ.
Applying the values of η and γ of (65) to (66), we have The number of outer iterations is bounded above by (see, Lemma II.17 in [13]) An upper bound for the total number of iterations is given by multiplying the number of outer iterations and the number of inner iterations, i.e., After some elementary reductions, we have the following theorems which yield the iteration bounds for the algorithms with large-and small-update methods.  iterations. The output gives an ε-approximate solution of (P) and (D).

5.
Numerical results. In this section, we report the computational performance of the algorithm in Figure 3 for For θ 1 = π 3 , θ 2 = π 4 , and θ 3 = π 6 , let the corresponding b be respectively.
The optimal solutions of the primal problems for θ i (i = 1, 2, 3) are given by Table  1.   Notice that the third duality gap 1.4205 × 10 −5 does not reach the ideal accuracy 10 −6 . In fact, for θ 3 , when µ reach the ideal accuracy the errors and biases between QB −1 x(µ) and A T y(µ) + Bs(µ) − c are not enough small, which leads to the third duality gap in Table 3 not reaching the ideal accuracy.
In the test problems, we use the same threshold parameter, the accuracy parameter and the update parameter as the Case I in the implementation. We choose .7876) T , respectively, as the starting points. Note that these points are also strictly feasible for θ i (i = 1, 2, 3), respectively. The initial value of the barrier parameter µ 0 = 1. And Ψ(x 0 , s 0 ; µ 0 ) = 0.5377 < τ.
The optimal solutions of the primal problems for θ i (i = 1, 2, 3) are given by Table  4.  Table 4. The optimal solutions for 6-dimensional CQO-CC And for the dual problems the optimal solutions are given by Table 5.   Table 5. The optimal solutions for 6-dimensional dual problems Furthermore, the following Table 6 shows the respective objective values and the corresponding duality gaps.   Table 6. The objective values and duality gaps for 6-dimensional CQO-CC Similar to Case I the third duality gap 3.5773 × 10 −5 does not reach the ideal accuracy 10 −6 for θ 3 < π 4 . Although the accuracy of the duality gap can not reach 10 −6 when θ < π 4 , it has already reach 10 −5 , which is also acceptable. The above numerical results show that the primal-dual interior-point algorithms can also apply to the CQO-CC problems as long as we transform the CQO-CC problems to the CQO-SOC problems. 6. Conclusions and remarks. In this paper, we have presented primal-dual IPMs for CQO-CC based on the kernel function of the classic barrier function. The obtained iteration bounds for large-and small-update methods are derived, namely, O(2 log 2 ε ) and O( √ 2 log 2 ε ), respectively, which are as good as the ones for the SOCO case. Although expected, these results were not obvious and, at certain steps of the analysis, they were not trivial or straightforward generalization of the SOCO case.
The research of CQO-CC problems is necessary and practical. Consider the friction cone constraint ||(f x , f y )|| = f 2 x + f 2 y ≤μf z in the grasp FOP (see [6]), which has been mentioned briefly, where f = (f x , f y , f z ) ∈ R 3 is the force applied at a contact point p andμ is the friction coefficient at the contact point p. Obviously, there exists a unique θ so that tan θ =μ. Thus, the friction cone constraint is a circular cone constraint and we can solve the grasp FOP by the algorithms we have presented. Some interesting topics for further research remain. Firstly, the search directions used in this paper are based on the NT-scaling scheme. It may be possible to design similar algorithms using other scaling schemes and still obtain polynomialtime iteration bounds. Secondly, the extension to the other NSCO problems deserve to be investigated.