An Inexact Interior-Point Lagrangian Decomposition Algorithm with Inexact Oracles

TL;DR

This paper introduces an inexact interior-point Lagrangian decomposition algorithm for convex optimization, combining dual decomposition, barrier smoothing, and proximal-Newton methods, with proven polynomial-time complexity and practical numerical validation.

Contribution

It presents a novel inexact interior-point method that handles inexact primal and dual solutions, extending the applicability and efficiency of convex optimization algorithms.

Findings

Algorithm achieves polynomial-time complexity.

Effective primal recovery strategy provided.

Numerical experiments demonstrate competitive performance.

Abstract

We develop a new inexact interior-point Lagrangian decomposition method to solve a wide range class of constrained composite convex optimization problems. Our method relies on four techniques: Lagrangian dual decomposition, self-concordant barrier smoothing, path-following, and proximal-Newton technique. It also allows one to approximately compute the solution of the primal subproblems (called the slave problems), which leads to inexact oracles (i.e., inexact gradients and Hessians) of the smoothed dual problem (called the master problem). The smoothed dual problem is nonsmooth, we propose to use an inexact proximal-Newton method to solve it. By appropriately controlling the inexact computation at both levels: the slave and master problems, we still estimate a polynomial-time iteration-complexity of our algorithm as in standard short-step interior-point methods. We also provide a strategy to recover primal solutions and establish complexity to achieve an approximate primal solution. We illustrate our method through two numerical examples on well-known models with both synthetic and real data and compare it with some existing state-of-the-art methods.

Figures

Tables

Table 1. Table 1: Numerical results of three solvers on 10 problem instances of ( 54 ).

Problem size		CPU time [s]			Feasibility violation			Objective value $f^{\star }$
​​$n$	​​$p$	​​IPLD	​​CP	​​SDPT3	​​IPLD	​​CP	​​SDPT3	​​IPLD	​​CP	​​SDPT3
​​96	​​17,686	​​0.70	​​3.80	​​4.24	​​4.747e-09	​​9.986e-08	​​0.000e+00	​​-60.1375	​​-60.1375	​​-60.1375
​​188	​​29,502	​​1.31	​​4.44	​​9.59	​​5.126e-09	​​9.974e-08	​​0.000e+00	​​-116.8216	​​-116.8216	​​-116.8216
​​239	​​38,050	​​1.85	​​6.64	​​11.75	​​7.227e-10	​​9.983e-08	​​0.000e+00	​​-171.4066	​​-171.4066	​​-171.4066
​​306	​​53,048	​​2.43	​​9.45	​​227.32	​​2.266e-10	​​9.995e-08	​​0.000e+00	​​-228.6001	​​-228.6001	​​-228.6001
​​242	​​72,016	​​2.78	​​10.25	​​809.57	​​9.055e-09	​​9.982e-08	​​0.000e+00	​​-288.6970	​​-288.6970	​​-272.1732
​​324	​​125,848	​​5.61	​​26.98	*	​​1.107e-08	​​9.987e-08	*	​​-569.2405	​​-569.2405	​​*
​​658	​​243,936	​​19.37	​​78.58	*	​​5.478e-09	​​9.987e-08	*	​​-1133.8747	​​-1133.8747	​​*
​​833	​​432,218	​​47.86	​​203.95	*	​​8.007e-08	​​9.999e-08	*	​​-2124.3265	​​-2124.3265	​​*
​​1,383	​​1,194,500	​​206.88	​​571.17	*	​​9.473e-08	​​9.983e-08	*	​​-3236.1724	​​-3236.1724	​​*
​​1,619	​​2,389,000	​​556.34	​​1297.86	*	​​4.422e-09	​​9.994e-08	*	​​-6474.3812	​​-6474.3812	​​*

Equations

$$\Vert u{\Vert }_x:={\left(u^{\top }{\nabla }^{2}f(x)u\right)}^{1/2},\text{and}\Vert v{\Vert }_x^{\ast }:={\left(v^{\top }{\nabla }^{2}f(x{)}^{-1}v\right)}^{1/2}.$$

$$P^{\star}:=\min_{x\in\mathcal{K}}\Big{\{}P(x):=g(x)+\phi(Ax)\Big{\}},$$

$$g^{\star}:=\min_{x\in\mathcal{K}}\Big{\{}g(x)~{}~{}~{}\mathrm{s.t.}~{}~{}Ax\in\mathcal{C}\Big{\}}.$$

$$g(x):=\sum _{i=1}^N{g}_i(x_i)\text{and}\mathcal{K}:={\mathcal{K}}_1\times \cdots \times {\mathcal{K}}_N,\text{vspace}-0.75ex$$

$$0\in \mathrm{ri}\left(\mathrm{dom}(\varphi )-A(\mathrm{dom}(g)\cap \mathcal{K})\right),\text{vspace}-0.75ex$$

$$D^{\star}:=\min_{y\in\mathbb{R}^{n}}\bigg{\{}D(y):=\underbrace{\max_{x\in\mathcal{K}}\bigg{\{}\langle Ax,y\rangle-g(x)\bigg{\}}}_{d(y)}~{}+~{}\phi^{\ast}(-y)\bigg{\}},\vspace{-0.75ex}$$

$$d(y):=\sum_{i=1}^{N}d_{i}(y)~{}~{}~{}~{}\text{with}~{}~{}~{}d_{i}(y):=\max_{x_{i}\in\mathcal{K}_{i}}\Big{\{}\langle A_{i}x_{i},y\rangle-g_{i}(x_{i})\Big{\}}.\vspace{-0.75ex}$$

$$\left\{\begin{array}[]{lll}0&\in\nabla{g}(x^{\star})-A^{\top}y^{\star}+\mathcal{N}_{\mathcal{K}}(x^{\star}),{}{}{}&\text{(primal optimality)}\vspace{1ex}\\ 0&\in y^{\star}+\partial{\phi}(Ax^{\star})&\text{(dual optimality)}\vspace{1ex}\\ x^{\star}&\in\mathcal{K},&\text{(primal feasibility)}.\end{array}\right.\vspace{-0.75ex}$$

$$0\in A{x}^{\star }-\partial {\varphi }^{\ast }(-y^{\star })\equiv \nabla d(y^{\star })-\partial {\varphi }^{\ast }(-y^{\star }).\text{vspace}-0.75ex$$

$$|⟨{\nabla }^{3}f(x)[v]u,u⟩|\le M_f\Vert u{\Vert }_x^2\Vert v{\Vert }_x^{\theta -2}\Vert v{\Vert }_2^{3-\theta },x\in \mathrm{dom}(f),u,v\in {\mathbb{R}}^p,\text{vspace}-0.75ex$$

$$f^{\ast }(y):=\underset{x\in \mathrm{dom}(f)}{\sup }\left\{-y^{\top }x-f(x)\right\},$$

$${\mu }_g:=\min \left\{{\lambda }_{\min }({\nabla }^{2}g(x))\mid x\in \mathrm{dom}(f)\cap \mathrm{dom}(g)\right\}\in (0,+\infty ),$$

$$\underset{u\in {\mathbb{R}}^p}{\sup }\left\{\nabla f(x{)}^{\top }u-\Vert u{\Vert }_x^2\right\}\le {\nu }_f,\forall x\in \mathrm{dom}(f).$$

$${\!\!\!\!\!}\overline{D}_{t}^{\star}:={\!\!}\min_{y\in\mathbb{R}^{n}}\Big{\{}\overline{D}_{t}(y):={\!\!\!}\underbrace{\max_{x\in\mathrm{int}\left(\mathcal{K}\right)}{\!\!}\Big{\{}y^{\top}Ax-g(x)-tf(x)\Big{\}}}_{\bar{d}_{t}(y)}+\underbrace{\phi^{\ast}(-y)}_{\bar{h}(y)}\Big{\}}.{\!\!\!\!\!}\vspace{-0.75ex}$$

$$D_{t}^{\star}:=\min_{y\in\mathbb{R}^{n}}\Bigg{\{}D_{t}(y):=\underbrace{\tfrac{M_{t}^{2}}{4}\bar{d}_{t}(y)}_{d_{t}(y)}{~{}}+{~{}}\underbrace{\tfrac{M_{t}^{2}}{4}\bar{h}(y)}_{h_{t}(y)}\Bigg{\}}.\vspace{-0.75ex}$$

$$\psi_{t}(x;y):=\frac{M_{t}^{2}}{4}\big{[}g(x)+tf(x)-y^{\top}Ax\big{]}.\vspace{-0.75ex}$$

$$|u{|}_{x,t}:={\left(u^{\top }{\nabla }^2{\psi }_t(x)u\right)}^{1/2},\text{and}|v{|}_{x,t}^{\ast }:={\left(v^{\top }{\nabla }^2{\psi }_t(x{)}^{-1}v\right)}^{1/2},$$

$$x^{\ast}_{t}(y):=\mathrm{arg}\!\!\!\!\!\min_{x\in\mathrm{int}\left(\mathcal{K}\right)}\Big{\{}\psi_{t}(x;y):=\tfrac{M_{t}^{2}}{4}\big{[}g(x)+tf(x)-y^{\top}Ax\big{]}\Big{\}},$$

$$\nabla{\psi}_{t}(x^{\ast}_{t}(y);y)\equiv\tfrac{M_{t}^{2}}{4}\big{[}\nabla{g}(x^{\ast}_{t}(y))+t\nabla{f}(x^{\ast}_{t}(y))-A^{\top}y\big{]}=0,$$

$$\textbf{$($Exact oracles$)$:}~{}~{}\left\{\begin{array}[]{ll}d_{t}(y)&=-\psi_{t}(x^{\ast}_{t}(y);y),\vspace{1ex}\\ \nabla{d_{t}}(y)&=\tfrac{M_{t}^{2}}{4}Ax^{\ast}_{t}(y),\vspace{1ex}\\ \nabla^{2}{d_{t}}(y)&=\frac{M_{t}^{4}}{16}A\nabla^{2}{\psi_{t}}(x^{\ast}_{t}(y))^{-1}A^{\top}.\end{array}\right.$$

$$\|u\|_{y,t}:=\big{(}u^{\top}{\nabla}^{2}{d_{t}}(y)u\big{)}^{1/2}~{}~{}~{}\text{and}~{}~{}~{}\|v\|_{y,t}^{\ast}:=\big{(}v^{\top}{\nabla}^{2}{d_{t}}(y)^{-1}v\big{)}^{1/2}.$$

$$\textbf{$($Inexact oracles$)$:}~{}~{}\left\{\begin{array}[]{ll}\widetilde{d}_{t}(y)&=-\psi_{t}(\widetilde{x}^{\ast}_{t}(y);y),\vspace{1ex}\\ \widetilde{\nabla}{d_{t}}(y)&=\tfrac{M_{t}^{2}}{4}A\widetilde{x}^{\ast}_{t}(y),\vspace{1ex}\\ \widetilde{\nabla}^{2}{d_{t}}(y)&=\frac{M_{t}^{4}}{16}A\nabla^{2}{\psi_{t}}(\widetilde{x}^{\ast}_{t}(y))^{-1}A^{\top}.\end{array}\right.$$

$$|\!\|u|\!\|_{y,t}:=\big{(}u^{\top}\widetilde{\nabla}^{2}{d_{t}}(y)u\big{)}^{1/2}~{}~{}~{}~{}\text{and}~{}~{}~{}~{}|\!\|v|\!\|_{y,t}^{\ast}:=\big{(}v^{\top}\widetilde{\nabla}^{2}{d_{t}}(y)^{-1}v\big{)}^{1/2}.$$

$$\text{if}|\nabla {\psi }_t({\tilde{x}}_t^{\ast }(y);y){|}_{{\tilde{x}}_t^{\ast }(y),t}^{\ast }\le \frac{\delta }{1+\delta }\text{then}{\delta }_t(y):=|{\tilde{x}}_t^{\ast }(y)-x_t^{\ast }(y){|}_{{\tilde{x}}_t^{\ast }(y),t}\le \delta .$$

$$\left\{\begin{array}[]{ll}&0\leq\omega\left(\frac{\delta_{t}(y)}{1+\delta_{t}(y)}\right)\leq d_{t}(y)-\widetilde{d}_{t}(y)\leq\omega_{\ast}\left(\frac{\delta_{t}(y)}{1-\delta_{t}(y)}\right),\vspace{1.25ex}\\ &\big{(}1-\delta_{t}(y)\big{)}^{2}\widetilde{\nabla}^{2}{d_{t}}(y)\preceq\nabla^{2}{d_{t}}(y)\preceq\big{(}1-\delta_{t}(y)\big{)}^{-2}\widetilde{\nabla}^{2}{d_{t}}(y),\vspace{1.25ex}\\ &|\!\|\widetilde{\nabla}{d_{t}}(y)-\nabla{d_{t}}(y)|\!\|_{y,t}^{\ast}\leq\delta_{t}(y),\end{array}\right.$$

$$|\nabla g({\tilde{x}}_t^{\ast }(y))+t\nabla f({\tilde{x}}_t^{\ast }(y))-A^{\top }y{|}_{{\tilde{x}}_t^{\ast }(y),t}^{\ast }\le \frac{4\delta }{M_t^2(1+\delta )}.$$

$$\left\{\begin{array}[]{lll}&|A^{\top}\tilde{y}^{\star}-\nabla{g}(\tilde{x}^{\star})|_{\tilde{x}^{\star},t}^{\ast}\leq\varepsilon{}{}{}{}{}{}{}&\text{($\varepsilon$-primal optimality)},\vspace{1ex}\\ &r\in\tilde{y}^{\star}+\partial{\phi}(A\tilde{x}^{\star}+e){}{}{}{}{}&\text{($\varepsilon$-dual optimality)},\vspace{1ex}\\ &\tilde{x}^{\star}\in\mathrm{int}\left(\mathcal{K}\right){}{}{}{}{}{}&\text{(primal feasibility)},\vspace{1ex}\\ &~{}~{}|\!\|e|\!\|_{\tilde{y}^{\star},t}^{\ast}\leq\varepsilon~{}~{}~{}\text{and}~{}~{}~{}|A^{\top}r|_{\tilde{x}^{\star},t}^{\ast}\leq\varepsilon.\end{array}\right.$$

$$0\in \nabla d_t(y)+\partial h_t(y)=\frac{M_t^2}{4}A{x}_t^{\ast }(y)+\partial h_t(y).$$

$$0\in {\tilde{\nabla }}^2{d}_{t_{k+1}}(y^k)({\bar{y}}^{k+1}-y^k)+\tilde{\nabla }{d}_{t_{k+1}}(y^k)+\partial h_{t_{k+1}}({\bar{y}}^{k+1}).$$

$$Q_{t_{k+1}}(y):=⟨\tilde{\nabla }{d}_{t_{k+1}}(y^k),y-y^k⟩+\frac{1}{2}⟨{\tilde{\nabla }}^2{d}_{t_{k+1}}(y^k)(y-y^k),y-y^k⟩+h_{t_{k+1}}(y),$$

Full text

1 Introduction

The Lagrangian dual decomposition framework is a classical technique to handle constrained convex optimization problems with separable structures such as conic, multi-stage stochastic, network, and distributed optimization problems Bertsekas1989b ; Birge1985 ; Connejo2006 ; Dantzig1963 ; Palomar2006 . This approach has been incorporated with interior-point methods to obtain a dual decomposition interior-point framework in early 1990s Kojima1993 . Since then, many researchers have regularly applied this approach to different problems. For example, Fukuda2002 exploited this idea to develop a dual decomposition algorithm for semidefinite programming, and Zhao2005 considered this method for general convex and multi-stage stochastic programming. The authors in Necoara2009 further investigated the method from Zhao2005 to solve a more general class of problems and obtained more intensive and rigorous complexity guarantees. The work TranDinh2012c studied this framework under the effect of inexact oracle computed by inexactly solving the primal subproblems up to a given accuracy. Other related theoretical results include bitlislioglu2017interior ; fukuda2000interior ; gros2014newton ; halldorsson2003interior ; Pakazad2017 ; Shida2008 ; Yamashita2011 . In particular, Pakazad2017 solved loosely coupled problems using message passing, and bitlislioglu2017interior applied it to multi-agent optimization problems. However, none of these works has investigated general constrained composite convex optimization settings involving linear operators and allows both inexactness in the slave problems and master problem altogether. In addition, existing methods do not handle directly nonsmooth objectives but often introduce auxiliary variables to reformulate the underlying problem into a smooth problem which may significantly increase problem size and loose their theoretical guarantee.

Motivation and goals:

Although the Lagrangian decomposition method is classical, it is very useful to handle large-scale constrained convex problems with separable structure by means of parallel and distributed computational architectures. In this paper, we conduct an intensive study on the interior-point Lagrangian decomposition (IPLD) framework considered in many existing works, especially Necoara2009 ; TranDinh2012c ; Zhao2005 , from the following aspects.

• (a)

Firstly, we consider a more general problem class than Kojima1993 ; Necoara2009 ; TranDinh2012c ; Zhao2005 by handling directly a nonsmooth composite convex function with a linear operator (see (2)) instead of couple linear equality constraints as in existing methods by means of proximal Newton-type methods (see Subsection 5.1).

• (b)

Secondly, our method works with inexact oracles of the dual problem arising from inexactly solving the primal subproblems (the slave problems). We explicitly describe the range of accuracies to flexibly control the tolerance of the subproblems (see Subsection 4.2).

• (c)

Thirdly, we also exploit inexact proximal-Newton method to handle general nonsmooth terms of the dual problems.

• (d)

Fourthly, we provide a thorough analysis for both the primal and dual problems and derive concrete iteration-complexity bounds for our method.

• (e)

Finally, we incorporate our approach with a recent concept called “generalized self-concordance” developed in SunTran2017gsc to handle new applications.

We are interested in the class of constrained composite convex problems where $g$ is smooth and satisfies some additional properties so that existing methods often do not have a theoretical convergence guarantee. For instance, the objective function does not have Lipschitz gradient or is not “tractably proximal”. We also consider a generic convex set where the projection onto it may not be tractable to compute such as general polyhedra. Under such assumptions, our problem setting covers a wide range class of applications ranging from optimal control, operations research, and networks to machine learning, statistics, and signal processing BenTal2001 ; Boyd2004 . It also covers standard conic programming such as linear programming, second-order cone programming, and semidefinite programming.

Our contribution:

We exploit the approach from Kojima1993 ; Necoara2009 ; TranDinh2012c ; Zhao2005 to develop a new algorithm for solving a class of constrained convex optimization problems. The main idea is to smooth the dual problem using a self-concordant barrier function Nesterov1994 associated with the constraint set, and apply a path-following scheme to solve the smoothed dual problem. While Kojima1993 ; Necoara2009 ; Zhao2005 exactly follow this main stream, TranDinh2012c proposed another path-following scheme and analyzed its convergence under inexact computation. It also provides a strategy to recover an approximate primal solution from its approximate dual solution. Compared to TranDinh2012c , this work studies a much more general problem class than TranDinh2012c . In addition, it is different from existing works, including TranDinh2012c , in several aspects as previously mentioned. To this end, we can summarize our contribution as follows:

• (a)

We exploit the approach in Kojima1993 ; Necoara2009 ; TranDinh2012c ; Zhao2005 and combine it with recent new mathematical tools in Nesterov2004c ; SunTran2017gsc to develop a new algorithm. The new mathematical tools allow us to cover much broader class of models than Kojima1993 ; Necoara2009 ; TranDinh2012c ; Zhao2005 , and to analyze polynomial-time iteration-complexity. In addition, we handle a more general class of problems than Kojima1993 ; Necoara2009 ; TranDinh2012c ; Zhao2005 by allowing general composite convex objectives involving linear operators (see (2)).

• (b)

We propose a new inexact interior-point Lagrangian decomposition algorithm to solve this class of problems. Our algorithm can deal with inexact oracles of the dual problems arising from approximating the primal subproblem solutions. It also uses an inexact proximal-Newton scheme to approximate the search direction in the dual problem. We characterize explicitly the choice of all related parameters and accuracies based on our analysis.

• (c)

We establish a polynomial-time iteration-complexity estimate of our method to find an approximate optimal solution. Our algorithm can be viewed as a short-step interior-point methods for general convex problems involving Nesterov and Nemirovskii’s self-concordance structures. Our complexity bound is the same as in standard interior-point methods (up to a constant factor), while it is able to directly handle nonsmooth objective by means of proximal operator.

In addition to the above main contribution, let us highlight some technical contribution of our methods. Firstly, unlike other methods involving inexact oracles in the literature Devolder2010 , our inexact oracle is rendered from inexact solution of the subproblem. The accuracy level can be adaptively chosen instead fixing as in existing methods to flexibly trade-off the computation cost by choosing rough accuracy at the early iterations and decrease it in the last iterations. Secondly, solving the primal subproblem (slave problem) is reduced to solve a nonlinear equation instead of a general convex problem as in some existing decomposition methods. As a result, we can characterize an implementable criterion to control the inexactness of the primal subproblems by using Newton-type schemes. Thirdly, instead of using unspecified parameters such as the radius of quadratic convergence region and contraction factor, we compute these parameters explicitly using the theory of self-concordant barriers as often seen in interior-point methods Nesterov2004 ; Nesterov1994 . Finally, combining inexact oracle and inexact methods make our algorithm practical since this computation is unavoidable in iterative methods, especially, in decomposition approaches when handling complex models.

Paper organization:

The rest of this paper is organized as follows. Section 2 states the problem of interest, basic assumptions, and its dual form. Section 3 recalls some preliminary results on (generalized) self-concordance and self-concordant barriers Nesterov1994 . Section 4 focuses on barrier smoothing techniques and inexact oracles. Section 5 presents our main algorithm and its complexity analysis as well as convergence guarantees. Section 6 provides two numerical examples to verify the theoretical results. For the sake of presentation, we move all the technical proofs to the appendix.

2 Problem statement, basic assumptions, and dual formulation

Notation and terminologies:

We work with finite dimensional vector space $\mathbb{R}^{p}$ or $\mathbb{R}^{n}$ endowed with standard inner product $x^{\top}y$ or $\langle x,y\rangle$ and Euclidean norm $\left\|x\right\|_{2}:=\sqrt{x^{\top}x}$. We denote by $\mathbb{S}^{p}_{+}$ (resp., $\mathbb{S}^{p}_{++})$ the set of symmetric positive semidefinite matrices (resp., symmetric positive definite matrices). Given $H\in\mathbb{S}^{p}_{++}$, we define a weighted norm $\|u\|_{H}:=\left(u^{\top}Hu\right)^{1/2}$ and its dual norm $\|v\|_{H}^{\ast}:=\left(v^{\top}H^{-1}v\right)^{1/2}$ for any vectors $u$ and $v$ in $\mathbb{R}^{p}$. For $X,Y\in\mathbb{S}^{p}_{+}$, $X\preceq Y$ means that $Y-X\in\mathbb{S}^{p}_{+}$ and $X\succeq Y$ stands for $X-Y\in\mathbb{S}^{p}_{+}$.

Given a three-time differentiable and strictly convex function $f$, we define the following local norms for any $u$ and $v$ in $\mathbb{R}^{p}$:

[TABLE]

They also satisfy the Cauchy-Schwarz inequality, i.e. $u^{\top}v\leq\|u\|_{x}\|v\|_{x}^{\ast}$. We say that $f$ is $\mu_{f}$-strongly convex if $f(\cdot)-(\mu_{f}/2)\|\cdot\|^{2}$ remains convex. We also often use the following two convex functions: $\omega(\tau):=\tau-\ln(1+\tau)$ for $\tau\geq 0$, and $\omega_{\ast}(\tau):=-\tau-\ln(1-\tau)$ for $\tau\in[0,1)$. These functions are smooth and strictly convex. We also use $\mathcal{O}\left(\cdot\right)$ to denote big-O complexity notion.

2.1 The primal problem and basic assumptions

Consider the following constrained composite convex optimization problem:

[TABLE]

where $g:\mathbb{R}^{p}\to\mathbb{R}$ is a smooth and convex function, $\phi:\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$ is a proper, closed, and convex function, $A\in\mathbb{R}^{n\times p}$, and $\mathcal{K}$ is a nonempty, closed, and convex set in $\mathbb{R}^{p}$. As a special case of (2), if we choose $\phi:=\delta_{\mathcal{C}}$, the indicator of a nonempty, closed, and convex set $\mathcal{C}$ in $\mathbb{R}^{n}$, then (2) reduces to the following general constrained convex problem:

[TABLE]

Without loss of generality, we can also assume that $g$ and $\mathcal{K}$ possess a separable structure as follows:

[TABLE]

for $N\geq 1$, where $x_{i}\in\mathbb{R}^{p_{i}}$, $\mathcal{K}_{i}\subseteq\mathbb{R}^{p_{i}}$, and $\sum_{i=1}^{N}p_{i}=p$ for $i=1,\cdots,N$.

Note that the separable structure (4) frequently appears in graph and network optimization. It is also a natural structure in conic programming such as linear programming and monotropic programming Rockafellar1985 . Another example is convex empirical minimization models in statistical learning, which can also be reformulated into (2) by duplicating variables.

Basic assumptions:

Our approach relies on the following assumptions:

Assumption 2.1

The optimal solution set $\mathcal{X}^{\star}$ of (2) is nonempty, and hence the optimal value $P^{\star}$ is finite. The following Slater condition holds:

[TABLE]

where $\mathrm{ri}\left(\mathcal{Z}\right)$ is the relative interior of $\mathcal{Z}$, and $\mathrm{dom}(\cdot)$ is the domain of $(\cdot)$.

Assumption 2.2

The function $g$ is standard self-concordant as in Definition 1. $\mathcal{K}$ is endowed with a $\nu_{f}$-self-concordant barrier $f$ as in Definition 2 and $A$ is full-row rank.

Note that Assumption 2.1 is standard and required in any primal-dual optimization method to guarantee strong duality. Assumption 2.2 is also not restrictive. First, the self-concordance of $g$ can be relaxed to a broader class called generalized self-concordant function as shown in Proposition 1 with additional structures. Next, the full-row rankness of $A$ can always be obtained by eliminating redundant rows. Finally, the self-concordant barrier of $\mathcal{K}$ is always guaranteed under mild condition as discussed in Nesterov1994 .

Throughout this paper, we assume that both Assumptions 2.1 and 2.2 hold without recalling them in the sequel.

2.2 Dual problem and optimality condition

The dual problem associated with (2) can be written as

[TABLE]

where $\phi^{\ast}(\cdot):=\sup_{u}\left\{\langle\cdot,u\rangle-\phi(u)\right\}$ is the Fenchel conjugate of $\phi$. Under the separable structure (4), we can decompose the dual function $d$ into $N$ functions $d_{i}$ on smaller spaces $\mathbb{R}^{p_{i}}$. That is

[TABLE]

This computation can be carried out in parallel. Moreover, under Assumption 2.1, the dual optimal solution set $\mathcal{Y}^{\star}$ of (6) is nonempty, and the strong duality holds, i.e. $P^{\star}+D^{\star}=0$. The optimality condition of the primal problem (2) can be written as

[TABLE]

Under Assumption 2.1, (7) is the necessary and sufficient condition for $x^{\star}\in\mathcal{X}^{\star}$ to be a primal optimal solution of (2), and $y^{\star}\in\mathcal{Y}^{\star}$ to be a dual optimal solution of (6). Note that $0\in y^{\star}+\partial{\phi}(Ax^{\star})$ can be written as

[TABLE]

This is exactly the optimality condition of the dual problem (6). Our goal is to approximate a primal-dual solution of (2) and (6) in the sense of Definition 4.

3 Generalized self-concordance and self-concordant barriers

Let us review the theory of generalized self-concordant functions SunTran2017gsc and self-concordant barriers Nesterov2004 ; Nesterov1994 , which will be used in the sequel.

Generalized self-concordance and standard self-concordance:

Assume that $f:\mathrm{dom}(f)\subseteq\mathbb{R}^{p}\to\mathbb{R}$ is a three-time continuously differentiable convex function, i.e. $f\in\mathbb{C}^{3}(\mathrm{dom}(f))$, we use $\nabla^{3}{f}(x)[u]$ to denote the third order derivative of $f$ at $x\in\mathrm{dom}(f)$ along a direction $u\in\mathbb{R}^{p}$. We recall the following definition SunTran2017gsc .

Definition 1 (SunTran2017gsc )

A $\mathbb{C}^{3}$-convex function $f$ is said to be $(M_{f},\theta)$-generalized self-concordant with the parameter $M_{f}\geq 0$, and order $\theta>0$ if

[TABLE]

where we use the convention $\frac{0}{0}=0$ for the case $\theta<2$ and $\theta>3$. If $\theta=3$, then $f$ reduces to the self-concordant function defined by Nesterov and Nemirovskii in Nesterov1994 . If $\theta=3$ and $M_{f}=2$, then $f$ is said to be standard self-concordant.

Basic properties:

Basic and fundamental properties as well as examples of generalized self-concordant functions can be found in SunTran2017gsc . We recall the following Legendre conjugate of a generalized self-concordant function. Let $f:\mathrm{dom}(f)\to\mathbb{R}$ be an $(M_{f},\theta)$-generalized self-concordant function, we define

[TABLE]

the Legendre conjugate of $f$ (i.e. $f^{\ast}(-y)$ is the Fenchel conjugate of $f$). For generalized self-concordant functions and their conjugates, we have the following result.

Proposition 1

• $\mathrm{(a)}$

If $f$ is $(M_{f},\theta)$-generalized self-concordant with $\theta\in(0,3)$ and $\mu_{f}$-strongly convex w.r.t. the Euclidean norm $\|\cdot\|_{2}$, then $f$ is $\hat{M}_{f}$-self-concordant with $\hat{M}_{f}:=\mu_{f}^{\frac{\theta-3}{2}}M_{f}$.

• $\mathrm{(b)}$

If $f$ is an $(M_{f},\theta)$-generalized self-concordant function with $\theta\in[3,6)$, then its Legendre conjugate $f^{\ast}(\cdot)$ is also $(M_{f},\theta_{\ast})$-generalized self-concordant with $\theta_{\ast}:=6-\theta$.

• $\mathrm{(c)}$

Assume that $f$ is $M_{f}$-self-concordant on $\mathrm{dom}(f)$ and $g$ is nonlinear and $(M_{g},\theta)$-generalized self-concordant on $\mathrm{dom}(g)$ with $\theta\in(0,3]$. If $\mathrm{dom}(f)\cap\mathrm{dom}(g)$ is nonempty, closed, and bounded, then $h:=f+g$ is $M_{h}$-self-concordant with $M_{h}:=\max\left\{M_{f},\hat{M}_{g}\right\}$, where $\hat{M}_{g}:=\mu_{g}^{\frac{\theta-3}{2}}M_{g}$ and

[TABLE]

if $\theta<3$, and $\hat{M}_{g}:=M_{g}$ if $\theta=3$.

Proof

The proof of statements (a) and (b) can be found in (SunTran2017gsc, , Propositions 4 and 6). If $\mathrm{dom}(f)\cap\mathrm{dom}(g)$ is nonempty, closed, and bounded, then $g$ is also $\mu_{g}$-strongly convex on $\mathrm{dom}(f)\cap\mathrm{dom}(g)$ with $\mu_{g}$ defined by (11). Applying statement (a) to the strongly convex function $g$, we obtain statement (c). $\square$

Discussion:

Proposition 1 shows that the class of self-concordant functions can be extended to cover at least three classes of smooth convex functions. The first one is the class of smooth and strongly convex functions that is also generalized self-concordant as studied in SunTran2017gsc . In the case it is not strongly convex, one can add a small quadratic regularizer to obtain this property. The second class is the conjugate of generalized self-concordant functions with Lipschitz continuous gradient. The third class of functions is generalized self-concordant functions on bounded domain. We believe that these three classes of functions cover a sufficiently large class of applications, see SunTran2017gsc for more detailed examples and additional properties.

Standard self-concordant barriers:

Next, we recall the class of standard self-concordant barriers, and its properties.

Definition 2

Given a nonempty, closed, and convex set $\mathcal{K}$ in $\mathbb{R}^{p}$, we say that $f$ is a $\nu_{f}$-self-concordant barrier of $\mathcal{K}$ if $f$ is standard self-concordant on $\mathrm{dom}(f)\equiv\mathrm{int}\left(\mathcal{K}\right)$, $f(x)\to+\infty$ as $x$ approaches the boundary $\partial{\mathcal{K}}$ of $\mathcal{K}$, and

[TABLE]

The self-concordant barrier $f$ is said to be a logarithmically homogeneous self-concordant barrier if $f(\tau x)=f(x)-\nu_{f}\ln(\tau)$ for any $\tau>0$ and $x\in\mathrm{dom}(f)$.

Given a self-concordant barrier of $\mathcal{K}$, we define $x_{f}^{\star}:=\mathrm{arg}\!\displaystyle\min_{x\in\mathcal{K}}f(x)$ the analytical center of $\mathcal{K}$ if $x_{f}^{\star}$ exists. Clearly, if $\mathcal{K}$ is bounded, then $x_{f}^{\star}$ exists. In addition to these properties, we also have $\|x-x^{\star}_{f}\|_{x^{\star}_{f}}\leq\rho_{f}$ for any $x\in\mathrm{dom}(f)$, where $\rho_{f}:=\nu_{f}+2\sqrt{\nu_{f}}$ for general self-concordant barrier $f$ and $\rho_{f}:=\nu_{f}$ if $f$ is logarithmically homogeneous.

4 Barrier smoothing technique and inexact oracles

In this section, we describe a barrier smoothing technique for (2) which has been used in Kojima1993 ; Necoara2009 ; Nesterov2011c ; TranDinh2012c ; Zhao2005 . Without loss of generality, we can assume that $M_{g}=2$, since any self-concordant function $g$ with the parameter $M_{g}>0$, $(M_{g}^{2}/4)g$ is standard self-concordant.

4.1 Smoothed dual problem

Under Assumption 2.2, we consider the following self-concordant barrier smoothed dual problem of (2) (shortly, smoothed dual problem):

[TABLE]

Note that $g(\cdot)+tf(\cdot)$ is self-concordant with the parameter $M_{t}:=\max\left\{2,\frac{2}{\sqrt{t}}\right\}$ on $\mathrm{dom}(f)\cap\mathrm{dom}(g)$. To make it standard self-concordant, we rescale (13) as follows:

[TABLE]

From Necoara2009 or Zhao2005 , $d_{t}$ is standard self-concordant. Clearly, if $t\in(0,1]$, then $M_{t}=\frac{2}{\sqrt{t}}$. In this case, we have $\bar{d}_{t}(y)=td_{t}(y)$ and $\bar{h}(y)=th_{t}(y)$.

To evaluate the (normalized) smoothed dual function $d_{t}$ and its derivative, we consider the following standard self-concordant function:

[TABLE]

Primal local norms:

Note that \nabla^{2}{\psi_{t}}(x;y)=\tfrac{M_{t}^{2}}{4}\big{[}\nabla^{2}{g}(x)+t\nabla^{2}{f}(x)\big{]}=\nabla^{2}{\psi_{t}}(x) is symmetric positive definite on $\mathrm{dom}(g)\cap\mathrm{dom}(f)$ and independent of $y$. Therefore, we define the following local norms on the primal space:

[TABLE]

for any $u,v\in\mathbb{R}^{p}$. If $t\in(0,1]$, then |u|_{x,t}=\left(u^{\top}\big{[}\nabla^{2}{f}(x)+\frac{1}{t}\nabla^{2}{g}(x)\big{]}u\right)^{1/2}.

Exact oracles of the dual function $d_{t}$:

We can summarize the properties of $d_{t}$ defined in (14) into the following proposition which we omit the proof.

Proposition 2

Under Assumption 2.2, $\psi_{t}(\cdot;y)$ defined by (15) and $d_{t}(\cdot)$ defined by (14) are standard self-concordant. Moreover, if the following primal subproblem has optimal solution

[TABLE]

then its solution is unique. The optimality condition of this subproblem is

[TABLE]

which is necessary and sufficient for $x^{\ast}_{t}(y)$ to be optimal to (17). The function value and derivatives of $d_{t}$ in (14) can be evaluated as (see Nesterov1994 )

[TABLE]

Dual local norms:

Since $\nabla^{2}{d_{t}}(y)\succ 0$, we can define the following local norms in the dual space:

[TABLE]

4.2 Inexact oracles of the smoothed dual function

When $g$ and $\mathcal{K}$ are not trivial, solving the smoothed slave subproblem (17) exactly is impractical. We can only approximately solve (17) or (18) up to a given accuracy as defined in the following.

Definition 3

Let $x^{\ast}_{t}(y)$ be the exact solution of (17) at $y\in\mathbb{R}^{n}$. We call $\widetilde{x}^{\ast}_{t}(y)$ a $\delta$-(approximate) solution of (17) if $\delta_{t}(y):=|\widetilde{x}^{\ast}_{t}(y)-x^{\ast}_{t}(y)|_{\widetilde{x}^{\ast}_{t}(y),t}\leq\delta$, where $|\cdot|_{x,t}$ is defined by (16).

Given an inexact solution $\widetilde{x}^{\ast}_{t}(y)$ of (17) as defined in Definition 3, we define an inexact oracle of $d_{t}$ as follows:

[TABLE]

Since $\nabla^{2}\psi_{t}(\cdot)$ is positive definite and $A$ is full-row rank, $\widetilde{\nabla}^{2}{d_{t}}(y)$ is positive definite. Now we define the following local norms using inexact oracles:

[TABLE]

We first prove some properties of inexact solution $\widetilde{x}^{\ast}_{t}(y)$ and inexact oracles of $d_{t}$ defined by (21) in the following proposition, whose proof can be found in Appendix A.1.

Proposition 3

For any $\delta\in[0,1)$, we have:

[TABLE]

In addition, for $d_{t}$ and its derivatives defined by (19), and its inexact oracle defined by (21), the following properties hold

[TABLE]

where $\omega(\tau):=\tau-\ln(1+\tau)$ for $\tau\geq 0$ and $\omega_{\ast}(\tau):=-\tau-\ln(1-\tau)$ for $\tau\in[0,1)$.

Discussion: The first estimate (23) shows that to obtain an approximate solution $\widetilde{x}^{\ast}_{t}(y)$ such that $|\widetilde{x}^{\ast}_{t}(y)-x^{\ast}_{t}(y)|_{\widetilde{x}^{\ast}_{t}(y),t}\leq\delta$, we need to solve the slave problem (17) such that

[TABLE]

This condition is implementable, e.g., when we apply a Newton-type method to solve the nonlinear system (18). The estimates in (24) show us how the inexact oracles in (21) approximate the exact ones in (19).

Approximate primal-dual solutions:

Given an accuracy $\varepsilon>0$, our goal is to compute an $\varepsilon$-approximate primal-dual solution $(\tilde{x}^{\star},\tilde{y}^{\star})$ to $(x^{\star},y^{\star})$ of (7) in the following sense:

Definition 4

A pair $(\tilde{x}^{\star},\tilde{y}^{\star})$ is called an $\varepsilon$-approximate primal-dual solution to an exact primal-dual one $(x^{\star},y^{\star})$ of (7) if

[TABLE]

Here, the errors are measured through local norms in primal and dual spaces defined in (16) and (22). These norms are computable since they are defined through $\tilde{x}^{\star}$ and $\tilde{y}^{\star}$. In addition, since $A$ is full-row rank, $A^{\top}r=0$ if and only if $r=0$. Because $\tilde{x}^{\star}\in\mathrm{int}\left(\mathcal{K}\right)$, we have $\mathcal{N}_{\mathcal{K}}(\tilde{x}^{\star})=\left\{\boldsymbol{0}\right\}$. Therefore the first line of (26) can approximate the first line of (7). Similarly, the second line of (26) approximates the second line of (7), i.e. $0\in y^{\star}+\partial{\phi}(Ax^{\star})$. Therefore, Definition 4 is consistent with the optimality condition (7).

5 Inexact IPLD Method with Inexact Oracles

We develop an inexact interior-point Lagrangian decomposition method to solve (2) by using the inexact oracles (21).

5.1 Inexact proximal-Newton method for (14)

The optimality condition of (14): Recall the smoothed dual problem (14), its optimality condition is

[TABLE]

Any $y_{t}^{\ast}$ satisfies (27) is an optimal solution of (14). The sequence $\left\{(x^{\ast}_{t}(y^{\ast}_{t}),y^{\ast}_{t})\right\}_{t\geq 0}$ forms a central path, which converges to $(x^{\star},y^{\star})$ a primal-dual solution of (2).

Exact Proximal Newton scheme: Suppose that we are currently at $y^{k}$, since $d_{t}$ is twice differentiable, we will apply proximal-Newton method to compute $\bar{y}^{k+1}$, which leads to

[TABLE]

If we define

[TABLE]

then we can write $\bar{y}^{k+1}:=\mathrm{arg}\!\displaystyle\min_{y}Q_{t_{k+1}}(y)$. Introducing the notation $\mathrm{prox}_{h_{t}}^{\widetilde{\nabla}^{2}{d_{t}}}(\cdot)$, we can write (28) in the following form (see TranDinh2016c for a concrete definition)

[TABLE]

Inexact Proximal Newton scheme: Similarly, we can also approximately solve (28) up to a given accuracy as.

[TABLE]

Here, the approximation “$:\approx$” is defined in the following sense:

Definition 5

For a given $\epsilon\geq 0$ and $Q_{t_{k+1}}$ defined by (29), a vector $y^{k+1}$ given in (31) is said to be an $\epsilon$-approximate solution to $\bar{y}^{k+1}$ of (28) if

[TABLE]

Note that (32) implies $|\!\|y^{k+1}-\bar{y}^{k+1}|\!\|_{y^{k},t_{k+1}}\leq\epsilon$. There exists several convex optimization methods to compute $y^{k+1}$ in (31). For example, we can apply accelerated proximal gradient methods such as FISTA Beck2009 ; Nesterov2007 to compute this point. We can also apply semi-smooth Newton-CG augmented Lagrangian methods in li2018highly ; zhao2010newton to solve this problem. We will discuss the computation of $y^{k+1}$ in detail in Section 6.

Generalized gradient mapping: Now let us define the following inexact generalized gradient mapping

[TABLE]

Using $\widetilde{\nabla}^{2}{d_{t}}(\cdot)$ defined by (21), we further define the following quantity:

[TABLE]

We call $\lambda_{t}(y)$ the inexact proximal-Newton decrement. In Subsection 5.4 we can show that this quantity can be used to characterize the optimality condition (7).

5.2 The algorithm

From the above analysis, we can combine all the steps together and describe an algorithm to solve (2) as in Algorithm 1. In this algorithm, we explicitly show how to choose the accuracy of inexact oracles and inexact proximal-Newton direction, and how to update the penalty parameter $t$.

Note that we have not specified how to find a starting point $(x^{0},y^{0})$ to guarantee (35) and how to set $k_{\max}$ in Algorithm 1. In Subsection 5.5, we will show that such an $(x^{0},y^{0})$ can be found in finite steps. In Subsection 5.4, we show how to set $k_{\max}$ to get an $\varepsilon$-approximate primal-dual solution of (2).

5.3 Convergence analysis

Our analysis consists of several steps and is organized as follows:

• •

Lemma 1 provides an estimate between $\lambda_{t_{k+1}}(y^{k+1})$ and $\lambda_{t_{k+1}}(y^{k})$ in (34).

• •

Lemma 2 bounds $\lambda_{t_{k+1}}(y^{k})$ in terms of $\tilde{\Delta}_{t_{k}}$, $\tilde{\Delta}_{t_{k+1}}$ and $\lambda_{t_{k}}(y^{k})$, where $\tilde{\Delta}_{t_{k}}$ and $\tilde{\Delta}_{t_{k+1}}$ measure the distances between $\tilde{x}_{t_{k}}^{\ast}(y^{k})$ and $\tilde{x}_{t_{k+1}}^{\ast}(y^{k})$.

• •

Lemma 3 shows how to upper bound $\tilde{\Delta}_{t_{k}}$ and $\tilde{\Delta}_{t_{k+1}}$.

• •

The main result of this section is Theorem 5.1 which provides an update rule of $t$ to maintain the point $y^{k}$ in the neighborhood of the central path. The proof of this theorem is obtained by combining all the above lemmas.

Firstly, we state the main estimate of the inexact Newton-type step at Step 6 of Algorithm 1 in Lemma 1, whose proof is given in Appendix A.2.1.

Lemma 1

Let $\left\{y^{k}\right\}$ be generated by Algorithm 1, and

[TABLE]

Then

[TABLE]

In particular, if $\tilde{\delta}_{k+1}=0$, $\delta_{k}=0$, and $\epsilon_{k}=0$, then (37) reduces to

[TABLE]

Note that if we solve both the slave problem at Step 4 and the master problem at Step 6 exactly, then we could obtain the estimate (38), which is the same as in standard interior-point path-following methods Nesterov2004 . Next, we show a relation between $\lambda_{t_{k+1}}(y^{k})$ and $\lambda_{t_{k}}(y^{k})$, whose proof is in Appendix A.2.2.

Lemma 2

Let $t_{k}$ be updated as $t_{k+1}:=\sigma t_{k}$ for given $\sigma\in(0,1)$. Define

[TABLE]

Then, the following estimate holds

[TABLE]

The following lemma shows how to bound $\tilde{\Delta}_{t_{k}}$ and $\tilde{\Delta}_{t_{k+1}}$, the distances between $\widetilde{x}^{\ast}_{t_{k}}(y^{k})$ and $\widetilde{x}^{\ast}_{t_{k+1}}(y^{k})$, whose proof is given in Appendix A.2.3.

Lemma 3

Let $\tilde{\Delta}_{t_{k}}$ and $\tilde{\Delta}_{t_{k+1}}$ be defined by (39), and $t_{k+1}:=\sigma t_{k}$ for some $\sigma\in(0,1)$. We define the following quantities:

[TABLE]

Then, we have

[TABLE]

where $\nu_{f}$ is the barrier parameter of $f$. In particular, for fixed $\delta\in[0,1)$, if we choose $\hat{\delta}_{t_{k}}\leq\delta$ and $\hat{\delta}_{t_{k+1}}\leq\delta$, then

[TABLE]

where c_{\nu}(\sigma):=\frac{\delta}{\sigma}+\big{(}\frac{1-\sigma}{\sigma}\big{)}\sqrt{\nu_{f}} is a decreasing function of $\sigma$ on $(0,1]$. As a consequence, we also have

[TABLE]

Utilizing the results of Lemma 1, Lemma 2 and Lemma 3, we can prove the following main result on the iteration-complexity of Algorithm 1.

Theorem 5.1

Let us choose $\beta\in(0,\frac{1}{10}]$. Suppose that we choose $\tilde{\delta}_{0},\delta_{k},\tilde{\delta}_{k+1},\epsilon_{k}\in[0,\frac{\beta}{100}]$ and update $t_{k}$ in Algorithm 1 as $t_{k+1}:=\sigma t_{k}$ with

[TABLE]

In addition, if $y^{0}\in\mathbb{R}^{n}$ and $x^{0}\in\mathbb{R}^{p}$ satisfy (35), then for all $k\geq 0$, we have

[TABLE]

Consequently, the number of iterations to obtain $t_{k}\leq\hat{\varepsilon}$ for a given $\hat{\varepsilon}>0$ and $\lambda_{t_{k}}(y^{k})\leq\beta$ does not exceed:

[TABLE]

where $\nu_{f}$ is the barrier parameter of $f$ and $t_{0}\in(0,1]$.

Proof

Let us first assume that $\lambda_{t_{k+1}}(y^{k})\leq 2.1\beta$. Using $\delta_{k},\tilde{\delta}_{k+1},\epsilon_{k}\in[0,10^{-2}\beta]$, after a few elementary calculations, we can overestimate (37) in Lemma 1 as

[TABLE]

when $\beta\in(0,\frac{1}{10}]$. Now, we prove that $\lambda_{t_{k+1}}(y^{k})\leq 2.1\beta$ is always satisfied. Indeed, since

[TABLE]

we can choose $\delta$ in Lemma 3 to be $\frac{\beta}{100}$. In addition, from $\sigma:=1-\frac{0.29\beta}{0.3\beta+\sqrt{\nu_{f}}}$, we can show that

[TABLE]

Next, using (44), we get

[TABLE]

Finally, combining these estimates and (40) we can show that

[TABLE]

when $\beta\in(0,\frac{1}{10}]$. Since $t_{k}:=\sigma^{k}t_{0}=\left(1-\frac{0.29\beta}{0.3\beta+\sqrt{\nu_{f}}}\right)^{k}t_{0}$, to guarantee $t_{k}\leq\hat{\varepsilon}$, we impose $\sigma^{k}t_{0}=\left(1-\frac{0.29\beta}{0.3\beta+\sqrt{\nu_{f}}}\right)^{k}t_{0}\leq\hat{\varepsilon}$. Note that $-\ln\left(1-\frac{0.29\beta}{0.3\beta+\sqrt{\nu_{f}}}\right)\sim\frac{0.29\beta}{0.3\beta+\sqrt{\nu_{f}}}\sim~{}\frac{1}{\sqrt{\nu_{f}}}$, we have

[TABLE]

as stated in (46). Here, $\sim$ means that two quantities can be approximated by the same order. $\square$

The worst-case iteration complexity:

Theorem 5.1 shows that for any $\hat{\varepsilon}>0$, the number of iterations $k$ to obtain $y^{k}$ such that $\lambda_{t_{k}}(y^{k})\leq\beta$ and $t_{k}\leq\hat{\varepsilon}$ does not exceed

[TABLE]

which is the same as in standard interior-point methods Nesterov2004 ; Nesterov1994 up to a constant factor. It depends on $\sqrt{\nu_{f}}$, where $\nu_{f}$ is the barrier parameter of $f$. Note that the parameter $\beta$ in Algorithm 1 represents the radius of the central path neighborhood as in standard path-following methods. While the range of $\beta$ in standard exact path-following methods Nesterov2004 is $(0,\frac{3-\sqrt{5}}{2}]$, it is $[0,\frac{1}{10}]$ in our method. Clearly, the latter is much smaller than the former one. However, this range was roughly estimated in our analysis and it is affected by the inexactness in our algorithm.

As we will show in Subsection 5.4, the conditions $\lambda_{t_{k}}(y^{k})\leq\beta$ and $t_{k}\leq\hat{\varepsilon}$ imply an approximate solution of (2) and (6).

5.4 Optimality certification

Our goal is to compute an approximate solution of the primal problem (2). The following theorem shows how we can find this approximate solution for both the primal and dual problem.

Theorem 5.2

Let $\{(x^{k+1},y^{k})\}$ be the sequence generated by Algorithm 1. Then, for $t_{k+1}\in(0,1]$ we have the following guarantees:

[TABLE]

Consequently, the number of iterations to obtain an $\varepsilon$-primal-dual solution $(x^{k+1},y^{k})$ in the sense of Definition 4 does not exceed:

[TABLE]

where $t_{0}\in(0,1]$ and $\nu_{f}$ is the barrier parameter of $f$.

Proof

From Step 4 of Algorithm 1, we can see that $x^{k+1}:=\widetilde{x}^{\ast}_{t_{k+1}}(y^{k})\in\mathrm{int}\left(\mathcal{K}\right)$. Moreover, Step 4 also leads to

[TABLE]

Next, for $t\in(0,1]$, it is obvious that \nabla^{2}{\psi_{t}}(x;y)=\tfrac{M_{t}^{2}}{4}\big{[}\nabla^{2}{g}(x)+t\nabla^{2}{f}(x)\big{]}=\nabla^{2}{f}(x)+\frac{1}{t}\nabla^{2}{g}(x). Consequently, one has $\nabla^{2}{\psi_{t}}(x;y)\succeq\nabla^{2}{f}(x)$. Using this fact, we can easily show that

[TABLE]

Combining this inequality and (50), we obtain the second estimate of (48).

Now, from (28), we have

[TABLE]

Using (21) and the definition of $h_{t}$, the last estimate becomes

[TABLE]

If we define $r^{k}:=y^{k}-\bar{y}^{k+1}$ and $e^{k}:=t_{k+1}\widetilde{\nabla}^{2}{d_{t_{k+1}}}(y^{k})(\bar{y}^{k+1}-y^{k})$, then the last expression leads to

[TABLE]

It is obvious to show that

[TABLE]

which is the first statement in the last line of (48).

Now, from (36) and $t_{k+1}\in(0,1]$, we have \widetilde{\nabla}^{2}d_{t_{k+1}}(y^{k})=\frac{1}{t_{k+1}^{2}}A\big{(}\nabla^{2}f(x^{k+1})+\frac{1}{t_{k+1}}\nabla^{2}g(x^{k+1})\big{)}^{-1}A^{\top}. This implies that

[TABLE]

Therefore, we have $|A^{\top}r^{k}|_{x^{k+1},t_{k+1}}^{\ast}=t_{k+1}\lambda_{t_{k+1}}(y^{k})$, which proves the second statement in the last line of (48).

From (48), to obtain an $\varepsilon$-primal-dual solution $(x^{k+1},y^{k})$ in the sense of Definition 4, we need to set $\left(\sqrt{\nu_{f}}+\frac{\delta_{k}}{1+\delta_{k}}\right)t_{k+1}\leq\varepsilon$ and $t_{k+1}\lambda_{t_{k+1}}(y^{k})\leq\varepsilon$. Since $\lambda_{t_{k+1}}(y^{k})\leq 2.1\beta$ (see the proof in Theorem 5.1) and $\delta_{k}\leq\frac{\beta}{100}$, we can set $t_{k+1}\leq\hat{\varepsilon}$ such that

[TABLE]

i.e., $\hat{\varepsilon}\leq\frac{\varepsilon}{(1+\sqrt{\nu_{f}})}$. Combining this expression and (46), we can show that the number of iterations to obtain an $\epsilon$-primal-dual solution does not exceed $\mathcal{O}\left(\sqrt{\nu_{f}}\ln\left(\frac{\sqrt{\nu_{f}}t_{0}}{\varepsilon}\right)\right)$, which is exactly (49). $\square$

Discussion:

Theorem 5.2 estimates the maximum iterations $k_{\max}$ to obtain an $\varepsilon$-primal-dual solution $(x^{k+1},y^{k})$ of (2) and (6). It shows that such a number of iterations remains the same as in standard path-following methods Nesterov2004 up to a constant factor. Although the norms in (48) are local norms, but this is the standard metric used in general interior-point methods Nesterov2004 ; Nesterov1994 .

5.5 Finding an initial point in Algorithm 1

We need to find $(x^{0},y^{0})$ such that the condition (35) holds. As in standard interior-point methods, we need to perform a damped proximal-Newton method. Such a method can be found in, e.g. TranDinh2013e ; Tran-Dinh2013a , but since we use inexact oracles, we need to customize this method in our context. More specifically, we describe this routine in Algorithm 2.

We terminate Algorithm 2 if we find $x^{0}:=\hat{x}^{j_{\max}}$ and $y^{0}:=\hat{y}^{j_{\max}}$ such that (35) holds. Since the constraint of $x^{0}$ in (35) is always satisfied from Step 3 of Algorithm 2, we only need to guarantee that $\lambda_{t_{0}}(y^{0})\leq\beta$.

The following theorem estimates the number of iterations to obtain $(x^{0},y^{0})$ satisfying (35).

Theorem 5.3

Let us define $\hat{\lambda}_{j}:=|\!\|\hat{s}^{j}-\hat{y}^{j}|\!\|_{\hat{y}^{j},t_{0}}$ and $\lambda_{j}:=|\!\|s^{j}-\hat{y}^{j}|\!\|_{\hat{y}^{j},t_{0}}$. Let $\{(\hat{x}^{j},\hat{y}^{j})\}$ be the sequence generated by Algorithm 2, where we choose \delta_{j},\epsilon_{j}\in\big{(}0,\frac{\beta}{100}\big{]} and the step-size

[TABLE]

Then, after at most finite number of iterations $j_{\max}$ as

[TABLE]

we obtain $y^{0}:=\hat{y}^{j_{\max}}$ and $x^{0}:=\hat{x}^{j_{\max}}$ such that $\lambda_{t_{0}}(y^{0})\leq\beta$ and (35) holds, where $y^{\ast}_{t_{0}}$ is the optimal solution of (14) at $t:=t_{0}$.

Proof

Note that at each iteration $j$ of Algorithm 2, we always have $\lambda_{j}>\beta$. By the triangle inequality and the choice of $\epsilon_{j}$, we can easily show that

[TABLE]

In addition, from Lemma 5 in Appendix A.3, we have

[TABLE]

where $\omega(\tau):=\tau-\ln(1+\tau)\geq 0$. Using $\epsilon_{j}\leq 10^{-2}\beta$, $\delta_{j}\leq 10^{-2}\beta$, and $\hat{\lambda}_{j}\geq(1-10^{-2})\beta$ in the above inequality, we get

[TABLE]

Summing up this inequality from $j=0$ to $j=j_{\max}$, we obtain

[TABLE]

which implies $j_{\max}\leq\frac{D_{t_{0}}(\hat{y}^{0})-D_{t_{0}}(y^{\star}_{t_{0}})}{\omega(0.97\beta(1-10^{-2}\beta))}$. Consequently, we obtain (53). $\square$

Discussion:

Theorem 5.3 shows that the number of iterations to obtain a starting point $(x^{0},y^{0})$ is finite even with inexact oracles and inexact proximal-Newton methods. However, the convergence rate of Algorithm 2 is sublinear in $j$. If $t_{0}$ is large (i.e., close to $1$), Algorithm 2 often requires a small number of iterations. Another possibility is to apply a path-following procedure as in TranDinh2015f to obtain a new variant with linear convergence rate. Note that the per-iteration complexity of Algorithm 2 is essentially the same as in Algorithm 1 since the computation of $\hat{\lambda}_{j}$ is neglectable. In particular, if we choose $\epsilon_{j}=\delta_{j}=0$, the steps size $\alpha_{j}$ will become the standard damped Newton step-size $\frac{1}{1+\lambda_{j}}$ in the theory of self-concordant function Nesterov1994 .

6 Numerical Experiments

We provide two numerical examples to illustrate our algorithm and compare it with some existing methods. We choose SDPT3 Toh2010 as a common used conic solver, and Chambolle-Pock’s (CP) primal-dual method Chambolle2011 as one of the most powerful first-order methods that can handle our problem. The first example is the well-known network utility maximization (NUM) problem, and the second one is the spectrum management problem for multi-user DSL networks studied in Tsiaflakis2008 . Our method and the CP method are implemented in Matlab 2018b, running on a Linux server with 3.4GHz Intel Xeon E5 and 16Gb memory.

6.1 Implementation remarks

We discuss how we implement two main steps of Algorithm 1 as follows. First, we need to solve the slave problem at Step 4 up to a given accuracy $\delta_{k}$ such that $\delta_{k}\leq 10^{-2}\beta$. Solving this problem is equivalent to solving the nonlinear equation $\nabla\psi_{t_{k+1}}(x;y^{k})=0$ in $x$. Since $\psi_{t_{k+1}}(\cdot;y^{k})$ is standard self-concordant, we can apply a damped-step Newton method to solve it. Combining this method and a warm-start strategy, we can solve this equation efficiently. Second, if $\phi=\delta_{\left\{b\right\}}$ in (2) for a given $b\in\mathbb{R}^{n}$, then the master problem at Step 6 reduces to a positive definite linear system $\widetilde{\nabla}^{2}d_{t_{k+1}}(y^{k})(y-y^{k})=-\widetilde{\nabla}d_{t_{k+1}}(y^{k})+0.25M_{t_{k+1}}^{2}b$, which can be efficiently solved by, e.g., preconditioned conjugate gradient methods. However, since $\phi$ usually does not have such a simple form, we need to apply iterative methods such as accelerated proximal gradient method Beck2009 ; Nesterov2004 to solve this problem which has a linear convergence rate. Note that we can also apply a semi-smooth Newton-type methods as in yang2015sdpnalp to solve this problem efficiently. In our numerical test, we use FISTA which seems working well.

6.2 Network Utility Maximization

Consider a network consisting of a finite set $\mathcal{S}$ of $N$ nodes and a finite set $\mathcal{E}$ of undirected capacitated edges. Let $x_{ij}$ denote the rate of sending data from node $i$ to node $j$. We assume that such a flow $f_{ij}$ from node $i$ to node $j$ is fixed and unique (we usually choose $f_{ij}$ to be the shortest path from $i$ to $j$).

Assume that each node $i$ is associated with a utility function u_{i}(x_{i}):=\log\big{(}d_{i}^{\top}x_{i}+\mu_{i}\big{)}, where $x_{i}:=(x_{i1},\cdots,x_{iN})^{\top}$, $d_{i}:=(d_{i1},\cdots,d_{iN})^{\top}$ and $\mu_{i}$ is a scalar. Since we ignore self-links from node $i$ to itself, we set $d_{ii}=0$ and $f_{ii}=\emptyset$. We further assume that the rate $x_{ij}$ is constrained to lie in a given interval $[0,M]$, where the scalar $M$ denotes the maximum capacity of flows.

Under this setting, we formulate the problem of interest into the following constrained convex optimization problem called NUM:

[TABLE]

Here, $L_{e}$ and $U_{e}$ are the lower bound and upper bound capacity of each edge, respectively, $r_{ij}$ is the initial designed rate from node $i$ to node $j$ and we do not want to have the rate $x_{ij}$ to be far away from our target $r_{ij}$, and $\rho$ is the corresponding penalty parameter to control the distance from $x_{ij}$ to $r_{ij}$. By defining $g(x):=-\sum_{i\in\mathcal{S}}\ln(d_{i}^{\top}x_{i}+\mu_{i})+\frac{\rho}{2}\|x_{i}-r_{i}\|^{2}$, $Ax=\sum_{e\in f_{ij}}x_{ij}$, $\phi(\cdot):=\delta_{[L_{e},U_{e}]}(\cdot)$, and $\mathcal{K}:=[0,M]$, we can reformulate (54) into (2). Clearly, this problem satisfies Assumptions 2.1 and 2.2.

We implement Algorithm 1 using Algorithm 2 to find an initial point using $t_{0}:=0.25$. We also implement the Chambolle-Pock method in Chambolle2011 and use SDPT3 to solve (54) as our competitors. Note that SDPT3 can directly handle $\log$-terms in $g$ compared to other interior-point solvers such as SeDuMi, SDPA, or Mosek. To avoid solving subproblems in the Chambolle-Pock method, we reformulate (54) by introducing auxiliary variables $z_{i}:=d_{i}^{\top}x_{i}+\mu_{i}$ for $i\in\mathcal{S}$. Since the Chambolle-Pock method has two step-sizes $\tau$ and $\sigma$, we tune $\tau$ for each run and let $\sigma:=0.99/(\tau\|K\|^{2})$, where $K$ is the linear operator obtained from reformulating (54) into a composite form. The best values of $\tau$ we found are between $10^{-6}$ and $10^{-7}$ depending on problem.

All algorithms are terminated when both infeasibility and relative duality gap reach $10^{-7}$ accuracy or the maximum number of iterations $k_{\max}:=20,000$ is exceeded. In the first case, we certify that the problem is “solved”, while in the second case, we mark it by “”. If problem is too big to solve by our computer, we also mark it by “”.

We use the “tech-router-rf” dataset from http://networkrepository.com/tech-routers-rf.php from rocketfuel , where we have approximately $2000$ nodes and $6000$ edges. In this network, each node is either a router or a computer IP. Each computer IP has to go through one or multiple routers to send data to another computer IP. For larger networks, we use the “tech-pgp” dataset from http://networkrepository.com/tech-pgp.php from boguna2004models , which is a social network with approximately $11000$ nodes and $24000$ edges. Given a network structure, we generate the input data as follows. The initial designed rate $r_{i}$ are generated from a uniform distribution $\mathcal{U}(0,1)$ between [math] and $1$. The upper and lower bounds of capacity are generated as $L_{e}:=(1-\mathcal{U}(0,0.5))\bar{b}$ and $U_{e}:=(1+\mathcal{U}(0,0.5))\bar{b}$, where $\bar{b}:=\sum_{i\in\mathcal{S}}A_{i}r_{i}$. The maximum limit of rate $M$ is $1$ and the penalty paramter $\rho$ is chosen to be $0.01$. Both $d_{i}$ and $\mu_{i}$ are generated randomly using $\mathcal{U}(0,1)$. To have different problem instances, we use different sub-networks of the original one.

We run three algorithms on $10$ problems instances of different sizes. The results are reported in Table 1, where $n$ is the number of linear inequality constraints, $p$ is the number of variables in (54), IPLD is Algorithm 1, and CP is the Chambolle-Pock method in Chambolle2011 .

From Table 1, we observe the following facts:

• •

IPLD can solve large-scale problems with huge variables and moderate number of couple linear inequality constraints relatively fast and accurate. IPLD outperforms SDPT3 and CP in a majority of problems in terms of CPU time and achieves the same accuracy in the objective value and constraint violation.

• •

It is not surprising that CP can also achieve high accuracy but requires very large number of iterations. The CP algorithm requires from $6500$ to $15200$ iterations to achieve our specified accuracy depending on problem.

• •

SDPT3 is quickly prohibited to handle larger instances due to the increase of variables and constraints when transforming it into a conic and $\log$ form. Therefore, the problem cannot be fit into our computer memory.

In summary, we believe that our method, IPLD, can potentially solve large-scale convex problems of the form (2) as long as they satisfy Assumptions 2.1 and 2.2. It can often achieve high accuracy within reasonably computational effort and can be easily parallelized. While primal-dual first-order methods require to tune the step-size to obtain good performance, our method is relatively robust to inexact oracles and inexact Newton-type methods as well as the choice of parameter $t_{0}\in(0,1]$.

6.3 Spectrum management of multi-user DSL networks

We consider the spectrum management problem of multi-user DSL networks studied in Tsiaflakis2008 , which can be cast into the following constrained problem:

[TABLE]

where $x_{i}\in\mathbb{R}^{m}$,  $a_{i}\in\mathbb{R}^{m}$, $c_{i}\in\mathbb{R}^{m}_{+}$, $b\in\mathbb{R}^{m}$, $L\in\mathbb{R}^{m}_{++}$, $g_{i}\in\mathbb{R}^{m}$, and $H_{i}\in\mathbb{R}^{m\times m}$. Here, $m$ is the number of users, and $M$ is the number of channels. For the detail explanation of this model, we refer the reader to Tsiaflakis2008 . Clearly, (55) can be cast into (2), where $g$ is self-concordant, $Ax=\sum_{i=1}^{M}x_{i}$, $\phi:=\delta_{(-\infty,b]}$ the indicator of $(-\infty,b]$, and $\mathcal{K}:=[0,L]^{M}$.

Our goal in this example is to verify the performance of Algorithm 1 using different accuracy levels both for inexact oracles and inexact proximal-Newton method. For this purpose, we use two real datasets to test our algorithm. More precisely, we first fix the tolerance $\delta_{k}$ of the inexact oracles at $10^{-5}$ and change the tolerance $\epsilon_{k}$ of the inexact proximal-Newton method from $10^{-2}$ to $10^{-11}$. Then, we fix the tolerance $\epsilon_{k}$ at $10^{-5}$ in the inexact proximal-Newton scheme and vary $\delta_{k}$ in the inexact oracles between $10^{-2}$ and $10^{-8}$. In all these cases, we terminate our algorithm whenever the feasibility violation is below $10^{-5}$ and the relative gap is below $10^{-6}$.

In the first test, we use a $7$-user asymmetric ADSL downstream dataset, where $m=7$ and $M=224$. Figure 1 shows how the number of iterations and the normalized CPU time depend on the tolerances, where the normalized CPU time is computed by $(T-T_{\min})/(T_{\max}-T_{\min})$ with the time $T$.

We can see from the top row of Figure 1 that with $\delta_{k}=10^{-5}$ fixed and $\epsilon_{k}\leq 10^{-4}$, the number of iterations is almost stable and the computational time does not decrease significantly. This suggests that the accuracy $\epsilon_{k}=10^{-4}$ is sufficiently for computing proximal-Newton direction in the dual problem. If $\epsilon_{k}>10^{-4}$, then the number of iterations and CPU time increase significantly. Similarly, if we fix $\epsilon_{k}=10^{-5}$ and increase $\delta_{k}$ from $10^{-8}$ to $10^{-2}$, then we can observe from the bottom row of Figure 1 that $\delta_{k}\leq 10^{-4}$ is sufficient to accommodate the inexact oracles.

To confirm our above statement, we again test our algorithm with the second dataset, 12-user VDSL upstream dataset, where $n=12$ and $M=1147$. Figure 2 provides the number of iterations and normalized CPU time by rescaling it between $[0,1]$ as in Figure 1. We again observe very similar behavior in both situations, but since the problem is relatively larger than that of the first dataset, the computational time increases significantly when we decrease the accuracy $\delta_{k}$ of the inexact oracles.

Acknowledgements.

This work was partly supported by the National Science Foundation (NSF), awarded number: DMS-1619884.

Appendix A Appendix: The proof of technical results in the main text

We provide all the missing proofs in the main text.

A.1 The proof of Proposition 3: Properties of inexact oracles.

Since $x^{\ast}_{t}(y)$ is the exact solution of (18), we have

[TABLE]

Therefore, using the standard self-concordance of $\psi_{t}$, we can show that

[TABLE]

By the Cauchy-Schwarz inequality, we have

[TABLE]

Combining the last two inequalities, we eventually get

[TABLE]

This implies that $|\widetilde{x}^{\ast}_{t}(y)-x^{\ast}_{t}(y)|_{\widetilde{x}^{\ast}_{t}(y),t}\leq\delta$.

Next, using (19) and (21), we have $d_{t}(y)-\widetilde{d}_{t}(y)=\psi_{t}(\widetilde{x}^{\ast}_{t}(y);y)-\psi_{t}(x^{\ast}_{t}(y);y)$. Therefore, applying (Nesterov2004, , Theorems 4.1.7 and 4.1.8) respectively, we obtain the first estimate of (24).

Note that since $\nabla^{2}{d_{t}}(y)=\frac{M_{t}^{4}}{16}A\nabla^{2}{\psi_{t}}(x^{\ast}_{t}(y);y)^{-1}A^{\top}$ and $\widetilde{\nabla}^{2}{d_{t}}(y)=\frac{M_{t}^{4}}{16}A\nabla^{2}{\psi_{t}}(\widetilde{x}^{\ast}_{t}(y))^{-1}A^{\top}$, using (Nesterov2004, , Theoryem 4.1.6), we obtain the second estimate of (24).

Finally, since $\nabla{d}_{t}(y)-\widetilde{\nabla}{d_{t}}(y)=\tfrac{M_{t}^{2}}{4}A(x^{\ast}_{t}(y)-\widetilde{x}^{\ast}_{t}(y))$, we have

[TABLE]

In the last inequality, we use $A^{\top}(AQ^{-1}A^{\top})^{-1}A\preceq Q$ for any symmetric positive definite matrix $Q$ and any full-row rank matrix $A$. Hence, we obtain the third estimate of (24). $\square$

A.2 The technical proofs of Subsection 5.3: Convergence analysis

We provide the proof of technical results in Subsection 5.3.

A.2.1 The proof of Lemma 1: Key estimate of the inexact PN scheme (31).

For simplicity of presentation, we redefine $t:=t_{k}$, $t_{+}:=t_{k+1}$, $y:=y^{k}$, $y_{+}:=y^{k+1}$, and $\bar{y}_{+}:=\bar{y}^{k+1}$, where $\bar{y}^{k+1}$ is defined by (28) or (30). Using these new notations, we also denote $\delta:=\delta_{t_{+}}(y)$, $\delta_{+}:=\delta_{t_{+}}(y_{+})$, $\epsilon:=\|\!|y_{+}-\bar{y}_{+}\|\!|_{y,t_{+}}$, and $\hat{\lambda}:=|\!\|y_{+}-y|\!\|_{y,t_{+}}$ to make our analysis more clean.

If we define $r_{t_{+}}(y):=\widetilde{\nabla}^{2}{d_{t_{+}}}(y)(\bar{y}_{+}-y)+\widetilde{\nabla}{d_{t_{+}}}(y)$, then from (28), we have

[TABLE]

which is equivalent to

[TABLE]

Utilizing the scaled proximal operator defined by (30), we can write the last statement as

[TABLE]

Using the definition of $\widetilde{G}_{t}(y)$ in (33) and of $\lambda_{t}(y)$ in (34), we can derive

[TABLE]

By the non-expansiveness of $\mathrm{prox}_{h_{t_{+}}}^{\widetilde{\nabla}^{2}{d_{t_{+}}}}(\cdot)$, see Tran-Dinh2013a , we can further estimate this term as

[TABLE]

Next, we decompose the following term $R_{t_{+}}(y)$ as

[TABLE]

Before we estimate the five terms of $R_{t_{+}}(y)$, we recall the following inequalities, which will be repeatedly used in our proof.

[TABLE]

[TABLE]

Here, the second last inequality of (60) is from (Nesterov2004, , Theoryem 4.1.6). Note that (60) also holds for $|\!\|\cdot|\!\|_{y_{+},t_{+}}$ and $\|\cdot\|_{y,t_{+}}$.

Using (60), we have

[TABLE]

Note that (61) also holds for $|\!\|\cdot|\!\|_{y_{+},t_{+}}$ and $|\!\|\cdot|\!\|_{y,t_{+}}$.

Now, we estimate the first term in $R_{t_{+}}(y)$ of (58) as

[TABLE]

For the second term of (58), we have

[TABLE]

To estimate the third term of (58), let S(y):=\big{[}\widetilde{\nabla}^{2}{d_{t_{+}}}(y)-\nabla^{2}{d_{t_{+}}}(y)\big{]}(y_{+}-y). We have

[TABLE]

However, $\|S(y)\|_{y,t_{+}}^{\ast}$ can be estimated as

[TABLE]

Using this estimate into (64), we finally get

[TABLE]

For the fourth term $\widetilde{\nabla}^{2}{d_{t_{+}}}(y)(\bar{y}_{+}-y_{+})$ of (58), we have

[TABLE]

Finally, we estimate last terms $\mathcal{T}_{5}:=|\!\|\nabla{d}_{t_{+}}(y_{+})-{\nabla}{d_{t_{+}}}(y)-\nabla^{2}{d_{t_{+}}}(y)(y_{+}-y)|\!\|_{y_{+},t_{+}}^{\ast}$. Note that

[TABLE]

where the second inequality follows from (TranDinh2016c, , Theorem 1).

Plugging (62), (63), (65), (66), and (67) into (58), we can estimate

[TABLE]

Note that

[TABLE]

Substituting this estimate and (68) into (57), we finally obtain

[TABLE]

However, from the definition of $\lambda_{t_{+}}(y)$, we have

[TABLE]

This implies $\hat{\lambda}:=|\!\|y_{+}-y|\!\|_{y,t_{+}}\leq\lambda_{t_{+}}(y){\,}+{\,}\epsilon$. Substituting this estimate into (69), we obtain (37). In particular, if $\delta=\delta_{+}=\epsilon=0$, then we can simplify (37) to obtain (38). $\square$

A.2.2 The proof of Lemma 2: The relationship between $\lambda_{t_{+}}(y)$ and $\tilde{\Delta}$.

We again redefine $t:=t_{k}$, $t_{+}:=t_{k+1}$, $y:=y^{k}$, $y_{+}:=y^{k+1}$, and $\bar{y}_{+}:=\bar{y}^{k+1}$ as in Lemma 1. In addition, we also define \bar{u}:=\mathrm{prox}_{h_{t}}^{\widetilde{\nabla}^{2}{d_{t}}(y)}\big{(}y-\widetilde{\nabla}^{2}{d_{t}}(y)^{-1}\widetilde{\nabla}{d_{t}}(y)\big{)}.

First, we show that $\lambda_{t_{+}}(y)$ and $\lambda_{t}(y)$ can be respectively expressed as

[TABLE]

If we denote by $\bar{h}(y):=\phi^{\ast}(-y)$, then $h_{t}(y)=\frac{M_{t}^{2}}{4}\bar{h}(y)$. By the definition of $\bar{u}$ and $\bar{y}_{+}$, we can write

[TABLE]

Using the monotonicity of $\partial{\bar{h}}(\cdot)$, we can show that

[TABLE]

Rearranging this inequality, we obtain

[TABLE]

By the Cauchy-Schwarz inequality, we can derive that

[TABLE]

To estimate $\mathcal{T}_{1}$, we first show the relationship between $\nabla^{2}{\psi_{t}}(\widetilde{x}^{\ast}_{t}(y))$ and $\nabla^{2}{\psi_{t_{+}}}(\widetilde{x}^{\ast}_{t_{+}}(y))$. Then, we use it to get the relationship between $\widetilde{\nabla}^{2}{d_{t}}(y)$ and $\widetilde{\nabla}^{2}{d_{t_{+}}}(y)$. Recall that \nabla^{2}{\psi_{t}}(x):=\frac{M_{t}^{2}}{4}\big{[}\nabla^{2}{g}(x)+t\nabla^{2}{f}(x)\big{]}. Moreover, if $t\in[0,1]$, then $\frac{M_{t}^{2}}{4}:=\max\left\{1,\frac{1}{t}\right\}=\frac{1}{t}$. Therefore, we can write

[TABLE]

For any $0\leq t_{+}\leq t\leq 1$, we have

[TABLE]

In addition, using the self-concordance of $\psi_{t}$ and (72), we also have

[TABLE]

If we take the inverse of both sides of (73), then we get

[TABLE]

Since $\widetilde{\nabla}^{2}{d_{t}}(y)=\frac{M_{t}^{4}}{16}A\nabla^{2}{\psi_{t}}(\widetilde{x}^{\ast}_{t}(y))^{-1}A^{\top}$ and $\widetilde{\nabla}^{2}{d_{t_{+}}}(y)=\frac{M_{t_{+}}^{4}}{16}A\nabla^{2}{\psi_{t_{+}}}(\widetilde{x}^{\ast}_{t_{+}}(y))^{-1}A^{\top}$, the last inequalities imply

[TABLE]

which are respectively equivalent to

[TABLE]

Since $\frac{M_{t}^{2}}{4}=\frac{1}{t}$ and $\frac{M_{t_{+}}^{2}}{4}=\frac{1}{t_{+}}$, we obtain from the above inequalities that

[TABLE]

Now we can estimate the first term $\mathcal{T}_{1}$ in (71) as

[TABLE]

To estimate the second term $\mathcal{T}_{2}$ of (71), by the definition of $\widetilde{\nabla}{d_{t}}$, we have

[TABLE]

Here, we use the fact that $A^{\top}(AQ^{-1}A^{\top})^{-1}A\preceq Q$ for any symmetric positive definite matrix $Q$ and any full-row rank matrix $A$.

Plugging (75) and (76) into (71), we get

[TABLE]

This inequality is equivalent to

[TABLE]

Finally, we can derive

[TABLE]

which is exactly (40) due to the update $t_{+}:=\sigma t$. $\square$

In order to prove Lemma 3 we need the following auxiliary result.

Lemma 4

Let $a\in(0,1)$ and $b\in(0,1)$ be two positive numbers such that $a+b<1$. Let

[TABLE]

Then, $\mathcal{N}(a,b)\subseteq\left\{(u,v)\in\mathbb{R}^{2}_{+}~{}\mid~{}u\leq\frac{a+b}{1-a-b},~{}v\leq\frac{a+b}{1-a-b}\right\}$.

Proof

Suppose $(u,v)\in\mathcal{N}(a,b)$ and $u>\frac{a+b}{1-a-b}$. Then, according to $\frac{u^{2}}{1+u}\leq au+bv$, we have

[TABLE]

Therefore, we can show that

[TABLE]

However, because $\frac{v^{2}}{1+v}\leq av+bu$, one can show that

[TABLE]

This contradicts (79). Consequently, we must have $u\leq\frac{a+b}{1-a-b}$. Use the symmetry between $u$ and $v$, we also have $v\leq\frac{a+b}{1-a-b}$. $\square$

A.2.3 The proof of Lemma 3: Upper bound on the solution difference $\tilde{\Delta}$.

First, by the self-concordance of $\psi_{t}$, we have

[TABLE]

Next, since $\nabla{\psi_{t}}(x)=\frac{1}{t}\nabla{g}(x)+\nabla{f}(x)$ and $\nabla{\psi_{t_{+}}}(x)=\frac{1}{t_{+}}\nabla{g}(x)+\nabla{f}(x)$, we have

[TABLE]

Therefore, we can bound

[TABLE]

Now, since $\nabla^{2}{\psi_{t_{+}}}(x)\succeq\nabla^{2}{f}(x)$, we can show that

[TABLE]

Substituting this and (82) into (81), and using the definition of $\hat{\delta}_{+}$ and $\hat{\delta}$, we get

[TABLE]

Finally, using $t_{+}=\sigma t$, we obtain the first estimate of (42) from the last inequality.

Similarly, by following the same argument as in the proof of the first estimate in (42), we can show that

[TABLE]

which is the second estimate of (42).

Assume that we choose $\hat{\delta}\leq\delta$ and $\hat{\delta}_{+}\leq\delta$ for some $\delta\in(0,1)$. Since $t_{+}=\sigma t$, if we denote by $c_{\nu}(\sigma):=\frac{\delta}{\sigma}+\frac{(1-\sigma)}{\sigma}\sqrt{\nu_{f}}\in(0,1)$. Assume further that $\delta+c_{\nu}(\sigma)<1$. Then, it is clear that $\frac{\hat{\delta}_{+}t_{+}}{t}+\frac{(t-t_{+})}{t}\sqrt{\nu_{f}}\leq c_{\nu}(\sigma)$ and $\frac{\hat{\delta}t}{t_{+}}+\frac{(t-t_{+})}{t_{+}}\sqrt{\nu_{f}}\leq c_{\nu}(\sigma)$. Applying Lemma 4, we can see that $(\tilde{\Delta},\tilde{\Delta}_{+})\in\mathcal{N}(\delta,c_{\nu}(\sigma))$. Hence, we have

[TABLE]

which proves (44). $\square$

A.3 The proof of result in Subsection 5.5: Finding an initial point.

The proof of Theorem 5.3 requires the following key lemma.

Lemma 5

Let $\left\{\hat{y}_{j}\right\}$ be the sequence generated by Algorithm 2, where the step-size $\alpha_{j}$ is chosen such that $\alpha_{j}\in(0,1]$ and $\frac{\alpha_{j}\hat{\lambda}_{j}}{1-\delta_{j}}<1$. Then, the following estimate holds

[TABLE]

where $D_{t}$ is defined in (14) and $\omega_{\ast}(\tau):=-\tau-\ln(1-\tau)$. The optimal step-size $\alpha_{j}$ that minimizes the right-hand side of (83) is

[TABLE]

The corresponding estimate from (83) with this step-size is

[TABLE]

In particular, if we set $\delta_{j}=\epsilon_{j}=0$, then we get the original damped-step proximal-Newton step-size $\alpha_{j}=\frac{1}{1+\lambda_{j}}$ and the estimate $D_{t_{0}}(\hat{y}^{j+1})\leq D_{t_{0}}(\hat{y}^{j})-\omega(\lambda_{j})$ for $\omega(\tau):=\tau-\ln(1+\tau)$.

Proof

Firstly, from the self concordance of $d_{t_{0}}$ defined in (14) and $\hat{y}_{j+1}=(1-\alpha)\hat{y}_{j}+\alpha\hat{s}_{j}$, we can show that

[TABLE]

Next, we will prove that

[TABLE]

where $\lambda_{j}:=|\!\|\hat{y}^{j}-s^{j}|\!\|_{\hat{y}_{j},t_{0}}$.

Indeed, by the Cauchy-Schwarz inequality, we have

[TABLE]

Since

[TABLE]

we have

[TABLE]

Using $0\in\widetilde{\nabla}d_{t_{0}}(\hat{y}^{j})+\widetilde{\nabla}^{2}d_{t_{0}}(s^{j}-\hat{y}^{j})+\partial h_{t_{0}}(s^{j})$, we can further estimate

[TABLE]

where $\nabla{h}_{t_{0}}(s^{j})\in\partial h_{t_{0}}(s^{j})$. Combining (89) and (90), we get

[TABLE]

Now, we can prove (87) as follows:

[TABLE]

Combining (86) and (87), and notice that $\omega_{\ast}(\alpha\|\hat{s}^{j}-\hat{y}^{j}\|_{\hat{y}_{j},t_{0}})\leq\omega_{\ast}(\frac{\alpha\hat{\lambda}_{j}}{1-\delta_{j}})$ we can deduce

[TABLE]

Using the fact that $\lambda_{j}\geq\hat{\lambda}_{j}-\epsilon_{j}$ and the definition $D_{t_{0}}:=d_{t_{0}}+h_{t_{0}}$, we obtain (83).

Next, if we maximize \zeta(\alpha):=\alpha\big{[}\hat{\lambda}_{j}^{2}-(\epsilon_{j}+\delta_{j})\hat{\lambda}_{j}\big{]}-\omega_{\ast}\left(\frac{\alpha\hat{\lambda}_{j}}{1-\delta_{j}}\right), we have $\alpha^{\star}:=\frac{(\hat{\lambda}_{j}-\epsilon_{j}-\delta_{j})(1-\delta_{j})^{2}}{\left[1+(1-\delta_{j})(\hat{\lambda}_{j}-\epsilon_{j}-\delta_{j})\right]\hat{\lambda}_{j}}$ as defined by (84). Plugging $\alpha^{\star}$ into $\zeta(\alpha)$, we get

[TABLE]

Since $\frac{x^{2}}{1+x}-\omega_{\ast}\left(\frac{x}{1+x}\right)=\omega(x)$, we finally have $\zeta(\alpha^{\star})=\omega\left((\hat{\lambda}_{j}-\epsilon_{j}-\delta_{j})(1-\delta_{j})\right)$, which proves (85). If $\delta_{j}=\epsilon_{j}=0$, then $\alpha_{j}$ reduces to $\frac{1}{1+\lambda_{j}}$ and we obtain $D_{t_{0}}(\hat{y}^{j+1})\leq D_{t_{0}}(\hat{y}^{j})-\omega(\lambda_{j})$ from (85). $\square$