A New Homotopy Proximal Variable-Metric Framework for Composite Convex Minimization

TL;DR

This paper introduces a novel homotopy proximal variable-metric framework for composite convex minimization, achieving finite iteration bounds and computational efficiency, especially in covariance estimation problems without matrix inversion.

Contribution

It develops a new parameterization for optimality conditions enabling homotopy proximal variable-metric methods and introduces a primal-dual-primal framework for proximal-Newton methods.

Findings

Achieves finite global iteration-complexity bounds under certain conditions.

Develops a primal-dual-primal framework that improves computational features.

Applies the method to covariance estimation without matrix inversion.

Abstract

This paper suggests two novel ideas to develop new proximal variable-metric methods for solving a class of composite convex optimization problems. The first idea is a new parameterization of the optimality condition which allows us to develop a class of homotopy proximal variable-metric methods. We show that under appropriate assumptions such as strong convexity-type and smoothness, or self-concordance, our new schemes can achieve finite global iteration-complexity bounds. Our second idea is a primal-dual-primal framework for proximal-Newton methods which can lead to some useful computational features for a subclass of nonsmooth composite convex optimization problems. Starting from the primal problem, we formulate its dual problem, and use our homotopy proximal Newton method to solve this dual problem. Instead of solving the subproblem directly in the dual space, we suggest to dualize this subproblem to go back to the primal space. The resulting subproblem shares some similarity promoted by the regularizer of the original problem and leads to some computational advantages. As a byproduct, we specialize the proposed algorithm to solve covariance estimation problems. Surprisingly, our new algorithm does not require any matrix inversion or Cholesky factorization, and function evaluation, while it works in the primal space with sparsity structures that are promoted by the regularizer. Numerical examples on several applications are given to illustrate our theoretical development and to compare with state-of-the-arts.

Tables

Table 1. Table 1: The information of binary classification datasets used in our experiments.

Dataset​	​​a1a​	​​a9a​	​​covtype​	​​news20​	​​mnist17​	​​mnist38​	​​rcv1​	​​real-sim​	​​w1a​	​​w8a​
#samples[$n$]​	​​30956​	​​16281​	​​581012​	​​19996​	​​317402​	​​292363​	​​677399​	​​72309​	​​47272​	​​14951​
#features[$p$]​	​​123​	​​122​	​​54​	​​1355191​	​​784​	​​784​	​​47236​	​​20958​	​​300​	​​300​

Table 2. Table 2: The performance and results of four algorithms on the logistic regression problem ( 47 ).

​​Datasets​	​​HomoPN​			​​HomoQuasiPN​			​​PG​			​​Ls-Rs-APG​			​​Sparsity​	​​$\rho$​
	​​iter	​​time[s]	​​rgap	​​iter	​​time[s]	​​rgap	​​iter	​​time[s]	​​rgap	​​iter	​​time[s]	​​rgap
​​a1a​​	​​7	​​ 0.247	​​0.3	​​ 13	​​ 0.190	​​0.3	​​1049	​​ 12.264	​​10.0	​​1088	​​ 16.272	​​7.6	​​11.38%	​​1e-2
​​a9a​​	​​6	​​ 0.153	​​0.4	​​ 13	​​ 0.169	​​1.1	​​1101	​​ 6.849	​​10.0	​​ 651	​​ 4.959	​​3.5	​​12.30%	​​1e-2
​​covtype​​	​​12	​​ 3.621	​​2.1	​​121	​​ 5.938	​​7.0	—	—	—	​​10508	​​3997.766	​​9.8	9.26%	​​1e-5
​​mnist17​​	​​6	​​ 101.088	​​0.3	​​ 24	​​25.030	​​6.2	—	—	—	​​ 912	​​2030.214	​​8.0	10.08%	​​2.5e-3
​​mnist38​​	​​4	​​ 115.304	​​6.1	​​ 14	​​36.109	​​9.7	—	—	—	​​13392	​​4151.542	​​10.0	11.86%	​​5.5e-3
​​news20​​	​​5	​​ 9.149	​​0.1	​​ 13	​​104.030	​​1.2	​​1307	​​ 71.579	​​10.0	​​ 90	​​ 57.225	​​3.4	​​11.04%	​​2.5e-6
​​rcv1​​	​​11	​​ 244.589	​​7.7	​​ 41	​​782.247	​​3.8	​​10291	​​3254.890	​​10.0	​​ 409	​​ 1222.426	​​5.0	​​10.14%	​​2.5e-6
​​real-sim​​	​​6	​​ 3.966	​​7.7	​​ 22	​​ 24.055	​​6.5	​​ 730	​​ 14.182	​​10.0	​​ 103	​​ 14.255	​​9.0	​​11.25%	​​2e-5
​​w1a​​	​​7	​​ 0.620	​​1.6	​​ 19	​​ 0.426	​​6.7	​​ 910	​​ 31.084	​​10.0	​​ 666	​​ 13.710	​​7.6	​​10.00%	​​1e-3
​​w8a​​	​​7	​​ 0.245	​​1.4	​​ 19	​​ 0.380	​​2.1	​​ 891	​​ 9.172	​​10.0	​​ 297	​​ 1.856	​​6.6	​​10.00%	​​1e-3

Table 3. Table 3: The performance of 4 4 4 algorithms on the D 𝐷 D -optimal experimental design problem ( 48 ).

Problem		(#Iterations) Time[s]				Objective value $F\left(x^k\right)$
${\chi }_i$	$p$	HomoPN	MUL[24]	IP[39]	SDPT3	HomoPN	MUL[24]	IP[39]	SDPT3
1	10000	(7)0.088	(2509)0.407	(127)0.262	0.399	20.51196	20.51254	20.51195	20.51241
1	50000	(7)0.368	(2510)1.050	(122)0.998	1.127	20.50907	20.50981	20.50907	20.50908
1	100000	(7)0.956	(2855)3.073	(120)1.962	2.008	20.50871	20.50943	20.50872	20.50883
2	10000	(7)0.057	(2493)0.571	(102)0.211	0.368	0.410236	0.410745	0.410221	0.410220
2	50000	(6)0.266	(3278)3.122	(101)0.867	1.280	0.409263	0.409964	0.409267	0.409260
2	100000	(5)0.525	(3910)9.361	(100)1.935	2.701	0.409143	0.409795	0.409154	0.409145
3	10000	(5)0.048	(1619)0.428	(97)0.236	0.415	5.142670	5.142919	5.142671	5.142670
3	40000	(5)0.214	(2208)2.378	(97)0.784	1.578	5.082114	5.082363	5.082119	5.082113
3	90000	(5)0.535	(2407)7.178	(95)1.774	3.827	5.062011	5.062261	5.062024	5.062011
4	10000	(6)0.069	(2512)0.488	(135)0.271	0.325	7.251897	7.252565	7.251890	7.251888
4	50000	(6)0.402	(3413)2.870	(130)1.045	1.064	7.251892	7.252527	7.251895	7.251889
4	100000	(6)1.032	(4013)8.480	(128)2.554	2.294	7.251891	7.252461	7.251902	7.251888

Table 4. Table 4: The performance of 3 3 3 algorithms on ( 48 ) for χ 1 ​ ( 8 ) subscript 𝜒 1 8 \chi_{1}(8) , χ 2 ​ ( 10 ) subscript 𝜒 2 10 \chi_{2}(10) and χ 3 ​ ( 10 ) subscript 𝜒 3 10 \chi_{3}(10) .

Problem		Time [s]			Objective value $F\left(x^k\right)$
$\chi$	$p$	HomoPN	MUL[24]	IP[39]	HomoPN	MUL[24]	IP[39]
1	10000	0.203	9.172	57.830	92.552013	92.552097	96.377758
1	50000	1.134	48.288	264.658	92.541322	92.541593	96.373413
1	100000	2.589	97.433	468.258	92.539162	92.540231	96.372813
2	10000	0.303	6.950	71.1327	18.424093	18.424563	22.387641
2	50000	1.555	43.340	229.223	18.416971	18.417293	22.388450
2	100000	2.862	84.389	543.588	18.416223	18.416356	22.389739
3	10000	0.188	17.657	2.140	30.158301	30.158329	30.158325
3	50000	0.822	72.969	240.569	29.956255	29.956513	37.689778
3	100000	2.217	165.105	556.789	29.889216	29.889565	37.692640

Table 5. Table 5: Poisson datasets information.

Dataset	vegas	games	news	911	facbook	wine
# samples	504	4940	39644	33059	99003	150930
# features	90	156	59	21	254	1154

Table 6. Table 6: The performance and results of three algorithms on the Poisson model ( 49 ).

Datasets	Time [s]			Number of iterations			$F\left(x^k\right)$
	HomoPN	BFGS	L-BFGS-LS-R	HomoPN	BFGS	L-BFGS-LS-R
vegas	0.04	0.23	0.21	5	11	11	4.2695e+00
games	0.12	1.34	0.99	6	22	22	4.5610e+00
news	1.82	52.94	48.57	9	220	215	1.1447e+02
911	0.10	0.87	0.55	6	15	15	8.2290e+00
facebook	1.54	62.36	64.68	8	39	48	4.6651e+01
wine	2.89	455.28	459.81	6	40	40	1.9819e+01

Equations

$${F^{\star}:=\min_{x\in\mathbb{R}^{p}}\Big{\{}{F(x)}:=f(x)+g(x)\Big{\}},}\vspace{-0.75ex}$$

$${\phi^{\star}:=\min_{X\succ 0}\Big{\{}\phi(X):=\mathrm{trace}\left(\Sigma X\right)-\log\det(X)+g(X)\Big{\}},}\vspace{-0.75ex}$$

$${\Vert u\Vert }_x:=⟨{\nabla }^{2}f(x)u,u{⟩}^{1/2}\text{and}{\Vert v\Vert }_x^{\ast }:=⟨{\nabla }^{2}f(x{)}^{-1}v,v{⟩}^{1/2},\text{vspace}-0.75ex$$

$${\mathrm{prox}^{H}_{g}(x):=\mathrm{arg}\!\displaystyle\min_{u\in\mathbb{R}^{p}}\big{\{}g(u)+\tfrac{1}{2}\|u-x\|_{H}^{2}\big{\}}.}\vspace{-0.75ex}$$

$$\Vert {\mathrm{prox}}_g^H(x)-{\mathrm{prox}}_g^H(y){\Vert }_H\le \Vert x-y{\Vert }_H,\text{vspace}-0.75ex$$

$${\begin{array}[]{ll}\|\nabla{g}(y)\|_{x}^{\ast}&=\max\left\{\langle\nabla{g}(y),z-y\rangle\mid\|z-y\|_{x}\leq 1\right\}\vspace{1ex}\\ &\leq\max\left\{|g(z)-g(y)|\mid\left\|z-y\right\|_{x}\leq 1\right\}\vspace{1ex}\\ &\leq L_{g}\max\left\{\left\|z-y\right\|_{x}\mid\left\|z-y\right\|_{x}\leq 1\right\}=L_{g}.\end{array}}\vspace{-0.75ex}$$

$$0\in \nabla f(x^{\star })+\partial g(x^{\star }).\text{vspace}-0.75ex$$

$$x^{\star }={\mathrm{prox}}_g^H\left(x^{\star }-H^{-1}\nabla f(x^{\star })\right).\text{vspace}-0.75ex$$

$$f_{\tau }(x):=\tau f(x)-(1-\tau )⟨{\xi }^0,x⟩,\text{vspace}-0.75ex$$

$${x_{\tau}^{\ast}=\mathrm{arg}\!\displaystyle\min_{x\in\mathbb{R}^{p}}\Big{\{}F_{\tau}(x):=f_{\tau}(x)+g(x)\Big{\}}.}\vspace{-0.75ex}$$

$$0\in \nabla f_{\tau }(x_{\tau }^{\ast })+\partial g(x_{\tau }^{\ast })\equiv \tau \nabla f(x_{\tau }^{\ast })-(1-\tau ){\xi }^0+\partial g(x_{\tau }^{\ast }),\text{vspace}-0.75ex$$

$$0\in \nabla f(x_{\tau }^{\ast })-\left(\frac{1}{\tau }-1\right){\xi }^0+\frac{1}{\tau }\partial g(x_{\tau }^{\ast }).\text{vspace}-0.75ex$$

$${x^{\ast}_{\tau}=\mathrm{prox}^{H}_{\frac{1}{\tau}g}\left(x^{\ast}_{\tau}-H^{-1}\big{(}\nabla{f}(x^{\ast}_{\tau})-(\tfrac{1}{\tau}-1)\xi^{0}\big{)}\right),}\vspace{-0.75ex}$$

$${G^{H}_{\tau}(x)=H\Big{(}x-\mathrm{prox}^{H}_{\tfrac{1}{\tau}g}\left(x-H^{-1}\big{(}\nabla{f}(x)-(\tfrac{1}{\tau}-1)\xi^{0}\big{)}\right)\Big{)}.}\vspace{-0.75ex}$$

$${x^{k+1}:\approx\mathrm{prox}_{\frac{1}{\tau_{k+1}}g}^{H_{k}}\Big{(}x^{k}-H_{k}^{-1}\left(\nabla{f}(x^{k})-\big{(}\tfrac{1}{\tau_{k+1}}-1\big{)}\xi^{0}\right)\Big{)}.}\vspace{-0.75ex}$$

$${\bar{x}^{k+1}:=\mathrm{arg}\!\displaystyle\min_{x\in\mathbb{R}^{p}}\Big{\{}\mathcal{P}_{k}(x):=\langle\nabla{f_{k\!+\!1}}(x^{k}),x-x^{k}\rangle+\tfrac{1}{2}\langle\nabla^{2}{f}(x^{k})(x-x^{k}),x-x^{k}\rangle+\tfrac{1}{\tau_{k+1}}g(x)\Big{\}},{\!\!\!}}\vspace{-0.75ex}$$

$${\bar{x}}^{k+1}={\mathrm{prox}}_{\frac{1}{{\tau }_{k+1}}g}^{{\nabla }^{2}f(x^k)}\left(x^k-{\nabla }^{2}f(x^k{)}^{-1}\nabla f_{k+1}(x^k)\right),\text{vspace}-0.75ex$$

$$x^{k+1}:\approx {\mathrm{prox}}_{\frac{1}{{\tau }_{k+1}}g}^{{\nabla }^{2}f(x^k)}\left(x^k-{\nabla }^{2}f(x^k{)}^{-1}\nabla f_{k+1}(x^k)\right).\text{vspace}-0.75ex$$

$${\mathcal{P}}_k(x^{k+1})-{\mathcal{P}}_k({\bar{x}}^{k+1})\le \frac{{\delta }_k^2}{2}.\text{vspace}-0.75ex$$

$$\delta (x^k):=\Vert x^{k+1}-{\bar{x}}^{k+1}{\Vert }_{x^k}\le {\delta }_k.\text{vspace}-0.75ex$$

$${\lambda }_{k+1}:=\Vert x^{k+1}-x_{{\tau }_{k+1}}^{\ast }{\Vert }_{x_{{\tau }_{k+1}}^{\ast }}\text{and}{\hat{\lambda }}_k:=\Vert x^k-x_{{\tau }_{k+1}}^{\ast }{\Vert }_{x_{{\tau }_{k+1}}^{\ast }}.\text{vspace}-0.75ex$$

$$C:=\frac{\Vert \nabla f(x^0)+{\xi }^0{\Vert }_2}{{\mu }_f}\text{and}\sigma :=\frac{1-{\tau }_0+{\tau }_0\omega \Gamma }{1-{\tau }_0+{\tau }_0\Gamma }\in (\omega ,1),\text{where}\Gamma :=\frac{\Vert \nabla f(x^0)+{\xi }^0{\Vert }_2}{\omega (L_g+\Vert {\xi }^0{\Vert }_2)}.\text{vspace}-0.75ex$$

$${\tau }_k:=1-\frac{(1-{\tau }_0){\sigma }^k}{{\tau }_0+(1-{\tau }_0){\sigma }^k}.\text{vspace}-0.75ex$$

$${x^{k+1}:=\mathrm{prox}_{\frac{\gamma}{\tau_{k+1}}g}\left(x^{k}-\gamma\big{(}\nabla{f}(x^{k})-(\tfrac{1}{\tau_{k+1}}-1)\xi^{0}\big{)}\right),}\vspace{-0.75ex}$$

$$\frac{2L_g(1-{\tau }_0)(1-\sqrt{\sigma })}{{\tau }_0-2L_g(1-{\tau }_0)(1-\sqrt{\sigma })}\le \frac{\sqrt{\sigma -0.01}}{10}-\frac{1}{18}.\text{vspace}-0.75ex$$

$${\tau }_{k+1}:=\left[1+\frac{{\Delta }_k{\tau }_k}{2L_g(1+{\Delta }_k)-{\Delta }_k{\tau }_k}\right]{\tau }_k\text{where}{\Delta }_k:=\left(\frac{1}{10}\sqrt{\sigma -0.01}-\frac{1}{18}{\sqrt{\sigma }}^k\right){\sqrt{\sigma }}^k.\text{vspace}-0.75ex$$

$$C_0:=\frac{\sqrt{\sigma -0.01}}{10}-\frac{1}{18}>0\text{and}\sigma \ge {\left(1-\frac{{\tau }_0{C}_0}{(1-{\tau }_0)(1+C_0)(\sqrt{{\nu }_f}+{\bar{c}}_0)}\right)}^2.\text{vspace}-0.75ex$$

$${\tau }_{k+1}:=\left(1+\frac{{\Delta }_k}{(1+{\Delta }_k)(\sqrt{{\nu }_f}+{\bar{c}}_0)}\right){\tau }_k,\text{with}{\Delta }_k:=\left(\frac{\sqrt{\sigma -0.01}}{10}-\frac{{\sqrt{\sigma }}^k}{18}\right){\sqrt{\sigma }}^k,\text{vspace}-0.75ex$$

$${\left\{\begin{array}[]{ll}\hat{s}^{j+1}&:\approx\mathrm{prox}_{\frac{1}{\tau_{0}}g}^{\nabla^{2}{f}{(\hat{x}^{j})}}\Big{(}\hat{x}^{j}-\nabla^{2}{f}(\hat{x}^{j})^{-1}\big{(}\nabla{f}(\hat{x}^{j})-(\frac{1}{\tau_{0}}-1)\hat{\xi}^{0}\big{)}\Big{)}\\ \hat{x}^{k+1}&:=(1-\alpha_{j})\hat{x}^{j}+\alpha_{j}\hat{s}^{j+1},~{}~{}\text{with}~{}~{}\hat{\zeta}_{j}:=\|\hat{s}^{j+1}-\hat{x}^{j}\|_{\hat{x}^{j}}~{}\text{and}~{}\alpha_{j}:=\frac{\hat{\zeta}_{j}-\hat{\delta}_{j}}{(1+\hat{\zeta}_{j}-\hat{\delta}_{j})\hat{\zeta}_{j}}.\end{array}\right.}\vspace{-0.75ex}$$

$$\lfloor \frac{F_{{\tau }_0}({\hat{x}}^0)-F_{{\tau }_0}(x_{{\tau }_0}^{\ast })}{\omega (0.9\beta )}\rfloor \text{vspace}-0.75ex$$

Full text

Keywords: Homotopy method; proximal variable-metric algorithm; global convergence rate; finite iteration-complexity; primal-dual-primal framework; composite convex minimization.

{AMS}

90C25, 90C06, 90-08

1 Introduction

Problem statement

We are interested in the following composite convex minimization template that covers various of applications in different fields including statistics, machine learning, image and signal processing, and engineering [3, 4, 8, 9, 15, 49]:

[TABLE]

where $f:\mathbb{R}^{p}\to\mathbb{R}\cup\{+\infty\}$ and $g:\mathbb{R}^{p}\to\mathbb{R}\cup\{+\infty\}$ are proper, closed, and convex functions. Here, $f$ often represents a loss function or a data fidelity term, while $g$ is considered as a regularizer or a penalty to promote some desired structures of the final solutions.

Motivation

This paper aims at addressing two questions arisen from numerical methods for solving (1). The first question concerns the global iteration-complexity of second-order-type methods. It is well-known that second-order methods such as Newton-type algorithms have fast local convergence rates under certain assumptions. In particular, the classical Newton method can achieve a local quadratic convergence rate under the local Lipschitz continuity of the Hessian around an optimal solution and the regularity of such an optimal solution [14]. However, global convergence behaviors as well as global convergence rates and iteration-complexity estimates of second-order-type methods have not yet been well understood. Recent attempts to address the aforementioned issues have been made for Newton-type methods [42, 44, 45], but they are still limited to some subclasses of problems such as self-concordant and global Lipschitz Hessian functions. In the first part of this paper, we address the following question.

• •

When can we design second-order-type methods that achieve global iteration-complexity?

Unfortunately, we do not have a complete answer for this question. However, we identify three different subclasses of (1) where we can develop new proximal variable-metric methods to achieve a global iteration-complexity. Our algorithms can solve nonsmooth instances of (1), but require $f$ to be smooth and satisfy other additional mild conditions that are different from existing methods.

In the second part of this paper, we address another situation of (1). We observe that existing methods for solving nonsmooth instances of (1) can be classified into the following categories:

1. (a)

If $f$ is smooth with the Lipschitz gradient and $g$ is nonsmooth but proximally tractable as defined in Subsection 2.2, then accelerated proximal-gradient methods achieve optimal convergence rate of $\mathcal{O}\left(\frac{1}{k^{2}}\right)$, where $k$ is the iteration counter. If $f$ is twice differentiable, and its Hessian is Lipschitz continuous [34] or self-concordant [41], then we can apply proximal-Newton methods [34, 59, 61] to efficiently solve (1). 2. (b)

If both $f$ and $g$ are proximally tractable, then operator splitting schemes such as Douglas-Rachford’s methods can be used to efficiently solve (1) but with a sublinear rate [3, 12]. 3. (c)

If $g(x)=\psi(Dx)$ for a given linear operator $D$, and both $f$ and $\psi$ are proximally tractable, then primal-dual methods such as Chambolle-Pock’s and primal-dual hybrid gradient methods, and alternating direction methods of multipliers (ADMM) can be applied to (1). These methods also achieve a sublinear rate in general.

We instead consider the following subclass of (1), where

1. (d)

$f$ is self-concordant as defined in Definition A.1; and $g$ is given by $g(x)=\psi(Dx)$, where $D$ is a linear operator, and $\psi$ is nonsmooth and convex, but proximally tractable.

Under this setting, existing methods such as proximal-gradient-type schemes are often not efficient for solving (1) due to the expensive evaluation of the proximal operator of $g$. We address the following research question:

• •

What is an appropriate solution method to solve (1) under the conditions stated in the subclass $(\mathrm{d})$?

This question may have multiple answers. One can apply some primal-dual methods to solve it. However, these methods only have a sublinear convergence rate. We instead propose a primal-dual-primal approach to solve (1) which consists of the following steps:

Construct the Fenchel dual problem of (1) when $g(x)=\psi(Dx)$. 2. 2.

Apply our homotopy proximal-Newton method in the first part to solve the dual problem. 3. 3.

Instead of solving the dual subproblem, we dualize it to go back to the primal space. 4. 4.

Construct an approximate primal solution of (1) from its dual approximate solution.

The idea of using primal-dual approach is classical, but our primal-dual-primal method has various computational advantages as well as a linear convergence rate when it is applied to the subclass (d) of (1). As a motivating example, we will show in Section 6.6 that this approach is very suitable for covariance estimation problem (2) below.

Examples

Apart from two research questions above, our paper is also motivated by several prominent applications. Let us recall a few concrete examples of (1):

*Covariance estimation models: * If $f(X):=-\log\det(X)+\mathrm{trace}\left(\Sigma X\right)$ in (1), where $\Sigma$ is a given symmetric matrix, then (1) covers both covariance and inverse covariance estimation problems in the literature depending on the choice of $g$ [2, 18, 32]:

[TABLE] 2. 2.

Poisson log-likelihood models: If we choose $f(x):=\sum_{i=1}^{n}\left(a_{i}^{\top}x-y_{i}\log(a_{i}^{\top}x)\right)$, where $\left\{(a_{i},y_{i})\right\}_{i=1}^{n}$ is a given dataset, then (1) covers Poisson log-likelihood models used in medical imaging, see, e.g., [35]. 3. 3.

Regularized logistic regression: If we choose $f(x):=\frac{1}{n}\sum_{i=1}^{n}\log(1+\exp(-y_{i}(a_{i}^{\top}x)))+\frac{\mu_{f}}{2}\|x\|_{2}^{2}$, where $\left\{(a_{i},y_{i})\right\}_{i=1}^{n}$ is a given dataset, and $\mu_{f}>0$ is a regularization parameter, then (1) covers the well-known logistic models including both sparse and group sparse settings under an appropriate choice of $g$. 4. 4.

Poisson regression: If $f(x):=\frac{1}{n}\sum_{i=1}^{n}\left(-y_{i}\exp(-\frac{1}{2}a_{i}^{\top}x)+\exp(\frac{1}{2}a_{i}^{\top}x)\right)+\frac{\mu_{f}}{2}\|x\|^{2}$, where $\left\{(a_{i},y_{i})\right\}_{i=1}^{n}$ is given, then we obtain a Poisson regression problem as studied in [27, 29]. 5. 5.

Distance-weighted discrimination (DWD): If $f(x):=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{(a_{i}^{\top}x+\mu_{i})^{q}}+\frac{\mu_{f}}{2}\|x\|_{2}^{2}$, for some fixed order $q>0$, then this model can be considered as a slight modification of the distance-weighted discrimination (DWD) for binary classification studied in [33, 40].

Many other applications of (1) that fit our assumptions can be found, for example in [11, 48, 61].

Literature review

Problem (1) is well studied in the literature under different assumptions on $f$ and $g$. Hitherto, several methods have been proposed for solving (1). Such methods include disciplined convex programming [22, 64], proximal gradient and accelerated proximal gradient [4, 41, 43], proximal Newton-type [6, 34, 61], splitting and alternating optimization [8, 15, 20, 68], primal-dual [9, 58], coordinate descent [16, 50, 46, 65], conditional gradient [23, 28], stochastic gradient-type methods [1, 13, 47, 53, 66], and incremental proximal gradient schemes [7].

Existing first-order methods for solving (1) heavily rely on the assumption that $f$ has Lipschitz gradient [41] and the proximal tractability of $g$ [3, 49] as defined in Subsection 2.2. Another common subclass of (1) is that $g(x)=\psi(Dx)$ for a given linear operator $D$, and both $f$ and $\psi$ are proximally tractable. Under this setting, operator splitting and primal-dual approaches can be applied to solve (1). Notable works in this direction include primal-dual hybrid gradient schemes, Chambolle-Pock’s methods, Douglas-Rachford and Vu-Condat splitting algorithms, and alternating direction methods of multipliers [3, 8, 9, 15, 21]. While first-order methods offer a low per-iteration computational complexity, they often require a large number of iterations and have a sublinear convergence rate. In addition, their efficiency also depends sensitively on the scaling and conditioning of the problem [19].

Proximal second-order methods such as proximal quasi-Newton [6, 30, 54] and proximal-Newton methods [34, 61] often achieve a high accuracy solution and have good local convergence rate but they usually have high per-iteration computational complexity. In proximal second-order-type methods, the trade-off between iteration-complexity and per-iteration computational complexity is crucial to obtain a good performance. Some existing works such as [6, 25, 26, 30, 59, 61] have provided evidence showing that second-order methods outperform first-order methods for some important subclasses of (1). The recent work [31] also studied the global linear convergence of Newton methods, but using a different concept called “$c$-Hessian stable”. Nevertheless, it is completely different from our approach.

Our approach

Our approach here relies on a combination of different ideas. The first idea is the homotopy method, which has been used in interior-point methods [41] and recently in path-following proximal Newton algorithms [60], where the main iterations rely on a scaled proximal Newton scheme [60]. The second idea is a new parameterization of the optimality condition of (1) as presented in Subsection 3.1. Our third idea is inspired by the generalized self-concordance concept introduced in [56]. The last one is a primal-dual-primal framework that we have mentioned above.

Our contribution

Our contribution can be summarized as follows.

• $\mathrm{(a)}$

We suggest a new parameterization for the optimality condition of (1) as a framework to study homotopy proximal variable-metric methods for solving different subclasses of (1). This framework covers homotopy proximal-gradient, proximal quasi-Newton, and proximal-Newton methods, and their inexact variants as special cases.

• $\mathrm{(b)}$

We propose a homotopy proximal variable-metric scheme, Algorithm 1, to solve (1) based on our new parameterization strategy. We show that this scheme achieves a global linear convergence rate under the strong convexity and Lipschitz gradient assumptions w.r.t. a local norm, and the Lipschitz continuity of $g$. We also propose an inexact homotopy proximal-Newton method to solve (1). Under the self-concordant property of $f$, and either the Lipschitz continuity of $g$ or the barrier property of $f$, our algorithm can also achieve a finite global iteration-complexity estimate. With an appropriate choice of initial points or suitable assumptions on $f$ and/or $g$, our method can achieve a linear convergence rate.

• $\mathrm{(c)}$

We propose a primal-dual-primal approach for a subclass of (1) where $f$ is self-concordant. This approach produces a new homotopy primal-dual proximal-Newton algorithm which can also achieve a linear convergence rate under given assumptions.

• $\mathrm{(d)}$

We specialize our algorithm to solve a special case (2) of (1) known as a regularized covariance estimation problems studied in [2, 18, 25, 26, 32, 59, 61]. This algorithmic variant possesses the following new features compared to existing works [25, 26, 59, 61]. First, it is applicable to any regularizer $g$ instead of just the $\ell_{1}$-norm as in [25, 26]. Second, it deals with the dual form of (2), while allowing one to reconstruct an approximate primal solution for (2). Third, it does not require any Cholesky factorization or matrix inversion as in [59] by working on the dual form. Fourth, the subproblem for computing proximal Newton directions is in the primal space of $X$ which has some special structures as promoted by the regularizer $g$ instead of in the dual space where the dual variable has structures that are correspondingly promoted by the conjugate $g^{\ast}$. The last point is important computationally when $g$ promotes the sparsity or low-rankness of the solutions.

Let us emphasize the following aspects of our contribution. Firstly, our new parameterization strategy can potentially be used to develop new numerical methods for different subclasses of (1) instead of only the four cases studied in this paper. Secondly, our path-following scheme for finding an appropriate initial point of Algorithm 1 is independent of a starting point as shown in Theorem 5.2. Thirdly, even for a strongly convex and Lipschitz gradient function $f$, our homotopy scheme has advantages in sparse optimization as discussed in Section 4.1. Fourthly, (1) is different from the barrier formulation considered in [60], where we do not use any penalty parameter for $f$ in (1) as compared to [60]. In addition, [60] is aimed at solving constrained convex optimization problems where the barrier is induced from the feasible set. Finally, for the covariance estimation problem (2), our method shares some similarity with [25, 26, 59, 60], but it is still fundamentally different. While [25, 26] focused on the sparse instance of (2), we consider a more general form in (1) that covers this example as a special case. Our algorithm and its convergence guarantee are completely different and rely on a different approach compared to [25, 26]. It has an iteration-complexity analysis for a genera $g$, while the analysis in [25, 26] rely critically on the special structure of the $\ell_{1}$-norm for $g$.

Paper organization

In Section 2, we recall some preliminary results used in this paper. Section 3 presents a new parameterization for the optimality condition of (1) and a conceptual three-stage proximal variable-metric framework, Algorithm 1, for solving (1). Section 4 analyzes the convergence of Algorithm 1 under three sets of assumptions. Section 5 proposes some procedures to find an appropriate starting point for Algorithm 1. Section 6 proposes a primal-dual-primal method for solving a nonsmooth subclass of (1), and its application to the covariance estimation problem (2). Section 7 provides several numerical experiments to illustrate our theoretical results. All technical proofs are deferred to the appendices.

2 Preliminaries: Scaled proximal operators and optimality condition

In this section, we recall some basic concepts which will be used in the sequel.

2.1 Basic notation and concepts

We work on the vector space $\mathbb{R}^{p}$ equipped with the standard inner product $\langle\cdot,\cdot\rangle$ and the corresponding Euclidean norm $\left\|\cdot\right\|_{2}$. We use $\mathbb{S}^{p}_{++}$ to denote the set of all symmetric positive definite matrices in $\mathbb{R}^{p\times p}$. For a given $H\in\mathbb{S}^{p}_{++}$, we use $\left\|x\right\|_{H}:=\langle Hx,x\rangle^{1/2}$ to denote the weighted norm. The corresponding dual norm is $\|y\|_{H}^{\ast}=\langle H^{-1}y,y\rangle^{1/2}$.

For a subset $\mathcal{X}$, $\mathrm{int}\left(\mathcal{X}\right)$ denotes the interior of $\mathcal{X}$, and $\partial{\mathcal{X}}$ denotes its boundary. Let $f:\mathbb{R}^{p}\to\mathbb{R}\cup\{+\infty\}$ be a convex function. As usual, $\mathrm{dom}(f)$ denotes the effective domain of $f$, and $\partial{f}$ denotes its subdifferential [51]. If $f$ is twice differentiable, then $\nabla{f}$ and $\nabla^{2}{f}$ denote its gradient and Hessian, respectively. For a given twice differentiable convex function $f$, if $x\in\mathrm{dom}(f)$ such that $\nabla^{2}{f}(x)\succ 0$, we define a local norm and its dual norm associated with $f$ as in [44]:

[TABLE]

for any $u,v\in\mathbb{R}^{p}$. Clearly, $\langle u,v\rangle\leq\left\|u\right\|_{x}\left\|v\right\|_{x}^{\ast}$. This is the weighted norm with $H=\nabla^{2}{f}(x)$. For a real number $a$, we use $\lfloor a\rfloor$ to denote the integer less than or equal to $a$. We use $[a]_{+}:=\max\left\{0,a\right\}$ for any real number $a$.

Given a nonempty convex set $\mathcal{X}$ in $\mathbb{R}^{p}$, and a point $x\in\mathbb{R}^{p}$, the distance from $x$ to $\mathcal{X}$ corresponding to the weighted-norm $\left\|\cdot\right\|_{H}$ is defined as $\mathrm{dist}\left(x,\mathcal{X}\right):=\inf_{y\in\mathcal{X}}\left\|x-y\right\|_{H}$. For a given convex function $f$, we say that $f$ is $\mu_{f}$-strongly convex if $f(\cdot)-\tfrac{\mu_{f}}{2}\left\|\cdot\right\|_{2}^{2}$ remains convex, where $\mu_{f}>0$ is called the strong convexity parameter of $f$. We say that $f$ is $L_{f}$-smooth (i.e., Lipschitz gradient continuous) if $f$ is differentiable on $\mathrm{dom}(f)$ and $\nabla{f}$ is Lipschitz continuous with a Lipschitz constant $L_{f}\in[0,+\infty)$, i.e., $\|\nabla{f}(x)-\nabla{f}(y)\|_{2}\leq L_{f}\|x-y\|_{2}$ for all $x,y\in\mathrm{dom}(f)$. We denote the class of $\mu_{f}$-strongly convex and $L_{f}$-smooth functions by $\mathcal{F}_{L,\mu}^{1,1}$. A convex function $g$ is Lipschitz continuous on $\mathrm{dom}(g)$ with a Lipschitz constant $L_{g}\in[0,+\infty)$ if $|g(x)-g(y)|\leq L_{g}\|y-x\|_{2}$ for all $x,y\in\mathrm{dom}(g)$.

2.2 Scaled proximal operators

Let $g:\mathbb{R}^{p}\to\mathbb{R}\cup\{+\infty\}$ be a proper, closed, and convex function, and $H\in\mathbb{S}^{p}_{++}$. We define the following scaled proximal operator [17] of $g$:

[TABLE]

The optimality condition of this minimization problem is $0\in H(\mathrm{prox}^{H}_{g}(x)-x)+\partial{g}(\mathrm{prox}^{H}_{g}(x))$, which can be written as $x\in(\mathbb{I}+H^{-1}\partial{g})(\mathrm{prox}^{H}_{g}(x))$, or $\mathrm{prox}^{H}_{g}(x)=(\mathbb{I}+H^{-1}\partial{g})^{-1}(x)$. When $H=\frac{1}{\gamma}\mathbb{I}$, where $\gamma>0$ and $\mathbb{I}$ is the identity matrix, $\mathrm{prox}^{H}_{g}(\cdot)$ becomes a classical proximal operator [3, 49], and is usually denoted by $\mathrm{prox}_{\gamma g}(\cdot)$. An important property of $\mathrm{prox}_{g}^{H}$ is its nonexpansiveness

[TABLE]

for any $x,y$ in $\mathbb{R}^{p}$. We say that $g$ is proximally tractable if $\mathrm{prox}^{H}_{g}(\cdot)$ can be efficiently evaluated, e.g., in a closed form or by a low-order polynomial time algorithm (e.g., $\mathcal{O}\left(p\log(p)\right)$). Computational methods for evaluating this scaled proximal operator and its classical forms can be easily found in the literature including [17, 49].

2.3 Lipschitz continuity w.r.t. local norm

Let $\|\cdot\|_{x}$ and its dual norm $\left\|\cdot\right\|_{x}^{\ast}$ be defined by a strictly smooth convex function $f:\mathbb{R}^{p}\to\mathbb{R}$, and $g:\mathbb{R}^{p}\to\mathbb{R}\cup\{+\infty\}$ be a proper, closed, and convex function.

Definition 2.1*.*

We say that $g$ is $L_{g}$-Lipschitz continuous w.r.t. $\|\cdot\|_{x}$ with a Lipschitz constant $L_{g}\in[0,+\infty)$, if for any $x,y,z\in\mathrm{dom}(f)\cap\mathrm{dom}(g)$, we have $|g(y)-g(z)|\leq L_{g}\|y-z\|_{x}$.

As a concrete example, assuming that $f(x)=\frac{1}{2}x^{\top}Qx+q^{\top}x$ is a strongly convex quadratic function, then $g$ is Lipschitz continuous in $\ell_{2}$-norm if and only if $g$ is Lipschitz continuous w.r.t. the local norm defined by $f$.

Lemma 2.2*.*

A proper, closed, and convex function $g$ is $L_{g}$-Lipschitz continuous w.r.t. $\|\cdot\|_{x}$ with a Lipschitz constant $L_{g}$ on $\mathrm{dom}(f)\cap\mathrm{dom}(g)$ if and only if $\|\nabla{g}(y)\|^{\ast}_{x}\leq L_{g}$ for any $x,y\in\mathrm{dom}(f)\cap\mathrm{dom}(g)$ and $\nabla{g}(y)\in\partial{g}(y)$.

In particular, if $f$ is strongly convex with a strong convexity parameter $\mu_{f}>0$ and $g$ is Lipschitz continuous in $\ell_{2}$-norm with a Lipschitz constant $\bar{L}_{g}\geq 0$ $($i.e. $|{g(y)-g(z)}|\leq\bar{L}_{g}\|y-z\|_{2}$ for any y,z\in\mathrm{dom}(f)\cap\mathrm{dom}(g)$$), then $g$ is $L_{g}$-Lipschitz continuous w.r.t. $\|\cdot\|_{x}$ on $\mathrm{dom}(f)\cap\mathrm{dom}(g)$ with the Lipschitz constant $L_{g}:=\frac{\bar{L}_{g}}{\sqrt{\mu_{f}}}$. However, the converse statement does not hold in general.

Proof 2.3*.*

For any $x,y\in\mathrm{dom}(f)\cap\mathrm{dom}(g)$ and $\nabla{g}(y)\in\partial{g}(y)$, we have

[TABLE]

Conversely, by convexity of $g$, we have $g(y)-g(z)\leq\langle\nabla{g}(y),y-z\rangle\leq\left\|\nabla{g}(y)\right\|_{x}^{\ast}\|y-z\|_{x}\leq L_{g}\|y-z\|_{x}$. By exchanging $y$ and $z$, we finally get $|g(z)-g(y)|\leq L_{g}\|z-y\|_{x}$.

If $f$ is strongly convex with a strong convexity parameter $\mu_{f}>0$, then we have $\nabla^{2}{f}(x)\succeq\mu_{f}\mathbb{I}$ for any $x\in\mathrm{dom}(f)$. Therefore, we have $\|y-z\|_{2}\leq\frac{1}{\sqrt{\mu_{f}}}\left\|y-z\right\|_{x}$. This shows that $|g(y)-g(z)|\leq\bar{L}_{g}\|y-z\|_{2}\leq\frac{\bar{L}_{g}}{\sqrt{\mu_{f}}}\left\|y-z\right\|_{x}$. Hence, $g$ is $L_{g}$-Lipschitz continuous w.r.t. $\|\cdot\|_{x}$ with $L_{g}:=\frac{\bar{L}_{g}}{\sqrt{\mu_{f}}}$.

As an example, if $\mathrm{dom}(g)$ is contained in an affine subspace defined by $\mathcal{L}:=\left\{x\in\mathbb{R}^{p}\mid Ax=b\right\}$, and $\nabla^{2}f(x)$ is uniformly positive definite on $\mathcal{L}$, then $g$ is $L_{g}$-Lipschitz and Lemma 2.2 still holds, but $f$ is still non-strongly convex. In this case, we say that $f$ is restricted strongly convex. Lemma 2.2 shows that the Lipschitz continuity w.r.t. the local norm $\|\cdot\|_{x}$ of $f$ is weaker than the global Lipschitz continuity of $g$ since we only require the condition to hold on $\mathrm{dom}(f)\cap\mathrm{dom}(g)$.

2.4 Fundamental assumption and optimality condition

Throughout this paper, we rely on the following fundamental assumption:

Assumption ​​ 1**.**

${\mathrm{dom}(F)}\!:=\!\mathrm{dom}(f)\cap\mathrm{dom}(g)\!\neq\!\emptyset$*. The solution set $\mathcal{X}^{\star}$ of (1) is nonempty. *

Assumption 1 is a standard one that is required in any solution method. Throughout this paper, we assume that Assumption 1 holds without recalling it.

The optimality condition associated with (1) becomes

[TABLE]

This condition is necessary and sufficient for $x^{\star}$ to be an optimal solution of (1). For any $H\in\mathbb{S}^{p}_{++}$, we can reformulate this optimality condition as a fixed-point condition:

[TABLE]

This formulation shows that $x^{\star}$ is a fixed point of $\mathbb{T}_{g}^{H}(\cdot):=\mathrm{prox}_{g}^{H}(\cdot-H^{-1}\nabla{f}(\cdot))$.

3 A Conceptual Homotopy Proximal Variable-Metric Framework

In this section, we introduce a novel parameterization of the optimality condition (6) and propose a conceptual framework for designing homotopy proximal variable-metric methods for solving (1).

3.1 Parametrization of the optimality condition

Given $x^{0}\in\mathrm{dom}(F)$, we compute a subgradient $\xi^{0}\in\partial{g}(x^{0})$. Then, we parameterize $f$ as follows:

[TABLE]

where $\tau\in[0,1]$. Clearly, $f_{1}(x)=f(x)$, $\nabla{f_{\tau}}(x)=\tau\nabla{f}(x)-(1-\tau)\xi^{0}$, and $\nabla^{2}{f_{\tau}}(x)=\tau\nabla^{2}{f}(x)$. In addition, $\mathrm{dom}(f_{\tau})=\mathrm{dom}(f)$ for any $\tau\in(0,1]$. Note that if we can choose $\xi^{0}\in\partial{g}(x^{0})$ such that $\xi^{0}=\boldsymbol{0}^{p}$, then $f_{\tau}(x)$ in (7) reduces to $f_{\tau}(x)=\tau f(x)$.

Next, we consider the following composite convex optimization problem derived from (1):

[TABLE]

This problem is similar to (1) and can be considered as a parametric perturbation instance of (1). The optimality condition of this problem is given by

[TABLE]

which is necessary and sufficient for $x^{\ast}_{\tau}$ to be an optimal solution of (8). We call this condition a parametric optimality condition of (1).

From the optimality condition (9), we can show that

• •

If $\tau=1$, then (9) becomes $0\in\nabla{f}(x^{\ast}_{1})+\partial{g}(x^{\ast}_{1})$, which is exactly the original optimality condition (6) of (1). Hence, $x^{\ast}_{1}=x^{\star}$ is an exact optimal solution of (1).

• •

If $\tau=0$, then (9) reduces to $\xi^{0}\in\partial{g}(x^{\ast}_{0})$. Therefore, we can choose $x^{\ast}_{0}=x_{0}$, the initial point, as an optimal solution of (8) at $\tau=0$.

Our main idea is to start from a small value $\tau_{0}\approx 0$ and follow a homotopy path on $\tau$ to find an approximate solution of $x^{\ast}_{\tau}$ at $\tau\approx 1$. As we will show later, we do not start from $\tau_{0}=0$, but from a sufficiently small value $\tau_{0}>0$.

Note that, when $\tau>0$, we can write (9) as

[TABLE]

If $g(\cdot)=\rho\|\cdot\|_{1}$, an $\ell_{1}$-regularizer, for a given regularization parameter $\rho>0$, then when $\tau$ is close to zero, the weight $\frac{1}{\tau}$ on $g$ is large and the solution of (8) is expected to be very sparse. This potentially can reduce the computational complexity of the underlying optimization method by working on sparse vectors or matrices. This property of the regularizer $g$ is also expected in other applications such as low-rank and group sparsity models.

Remark 3.1*.*

The formulation (9) is new and it does not reduce to any existing homotopy formulation including [41, Formula 4.2.26] to the best of our knowledge. This formulation is expected to lead to more efficient homotopy-type algorithms for solving sparse and low-rank convex optimization as explained above.

3.2 A fixed-point interpretation of the parametric optimality condition

Recall the optimality condition (9) given as $0\in\nabla{f_{\tau}}(x^{\ast}_{\tau})+\partial{g}(x^{\ast}_{\tau})$. By using the scaled proximal operator $\mathrm{prox}_{g}^{H}$, we can reformulate (9) into a fixed-point problem:

[TABLE]

for any $H\in\mathbb{S}^{p}_{++}$. Let us define the following mapping for any $x\in\mathrm{dom}(F)$:

[TABLE]

Clearly, (10) is equivalent to $G_{\tau}^{H}(x^{\ast}_{\tau})=0$. We call $G_{\tau}^{H}$ the scaled generalized gradient mapping of the parametric problem (8). The most common case is $H=\frac{1}{\gamma}\mathbb{I}$ as mentioned above for some $\gamma>0$. Then, $G_{\tau}^{H}$ reduces to the standard generalized gradient mapping [41].

3.3 Conceptual framework of homotopy proximal variable-metric methods

We first describe our conceptual three-stage proximal variable-metric algorithm as in Algorithm 1.

We will provide the details of each stage in the sequel based on some appropriate assumptions for (1). The main step of Algorithm 1 is (12), where we need to evaluate the scaled proximal operator $\mathrm{prox}_{\frac{1}{\tau}g}^{H}(\cdot)$. Depending on the choice of the variable matrix $H_{k}$, we obtain different methods:

• •

If $H_{k}$ is diagonal, then we obtain a homotopy proximal gradient method.

• •

If $H_{k}$ approximates $\nabla^{2}{f}(x^{k})$, then we obtain a homotopy proximal quasi-Newton method.

• •

If $H_{k}=\nabla^{2}{f}(x^{k})$, then we obtain a homotopy proximal Newton method.

The choice of an initial point $x^{0}$, the initial value $\tau_{0}$ of the parameter $\tau$, the update rule of $\tau_{k}$, and the approximation rule of (12) in Algorithm 1 will be specified in the sequel.

3.4 Inexact proximal Newton scheme

Let $H_{k}=\nabla^{2}{f}(x^{k})$. The exact evaluation of the scaled proximal operator in (12) is equivalent to solving the following convex subproblem:

[TABLE]

where \nabla{f}_{k+1}(x^{k}):=f(x^{k})-\big{(}\tfrac{1}{\tau_{k+1}}-1\big{)}\xi^{0}. In this case, we can write $\bar{x}^{k+1}$ as

[TABLE]

the exact solution of (13).

When $g$ is nontrivial (e.g., not a linear function), we can only approximate the true solution $\bar{x}^{k+1}$ of (13) by an approximation $x^{k+1}$ such that

[TABLE]

Here the approximation “$:\approx$” is defined explicitly next in Definition 3.2.

Definition 3.2*.*

Let $\bar{x}^{k+1}$ be the exact solution of (13), and $\delta_{k}\geq 0$ be a given accuracy. We say that $x^{k+1}$ is a $\delta_{k}$-approximate solution to $\bar{x}^{k+1}$, denoted by $x^{k+1}:\approx\bar{x}^{k+1}$ as in (14), if

[TABLE]

Using this definition, we have the following result, see [60, Lemma 3.2.].

Fact 1*.*

$\frac{1}{2}\|x^{k+1}-\bar{x}^{k+1}\|_{x^{k}}^{2}\leq\mathcal{P}_{k}(x^{k+1})-\mathcal{P}_{k}(\bar{x}^{k+1})$. Consequently, combining this inequality and (15), we can show that if (15) holds, then

[TABLE]

The condition (15) can be guaranteed by using several optimization methods in the literature such as accelerated proximal-gradient [4, 41, 52], ADMM [8], or semi-smooth Newton-CG augmented Lagrangian methods [67]. In Subsection 6.3, we approximately compute (13) via solving its dual.

We define the following local distances to measure the distance of approximations $x^{k+1}$ and ${x}^{k}$ to the true solution $x^{\ast}_{\tau_{k+1}}$ of the parameterized problem (9):

[TABLE]

These metrics will be used to analyze the convergence of our methods.

4 Convergence and iteration-complexity analysis

We analyze the convergence and iteration-complexity of Algorithm 1 for solving (1) under three different subclasses of $f$ and $g$.

4.1 Linear convergence for the smooth and strongly convex case

The first class of models is when $f$ and $g$ in (1) satisfies the following assumption.

Assumption ​​ 2**.**

*Assume that $f$ is $\mu_{f}$-strongly convex and $L_{f}$-smooth. The function $g$ is $L_{g}$-Lipschitz continuous on $\mathrm{dom}(g)$. *

Under Assumption 2, Algorithm 1 only has Stage 2 and we do not need to perform Stage 1 and Stage 3 of Algorithm 1. We can start from any starting point $x^{0}\in\mathrm{dom}(F)$. We show that Algorithm 1 achieves a global linear convergence rate. This result is stated in the following theorem, whose proof can be found in Apppendix B.1.

Theorem 4.1*.*

Under Assumption 2, let $m$ and $L$ be two constants such that $0<m\leq L<+\infty$ and $\omega:=\frac{1}{m}\sqrt{(L-2\mu_{f})m+L_{f}^{2}}<1$. For any given $\tau_{0}\in(0,1)$ and $x^{0}\in\mathrm{dom}(F)$, we define

[TABLE]

Let $\left\{(x^{k},\tau_{k})\right\}$ be the sequence generated by the exact scheme (12) in Algorithm 1, where $H_{k}\in\mathbb{S}^{p}_{++}$ is chosen such that $m\mathbb{I}\preceq H_{k}\preceq L\mathbb{I}$ and $\tau_{k}$ is updated by

[TABLE]

Then, for any $k\geq 0$, we have $\|x^{k}-x^{\ast}_{\tau_{k}}\|_{2}\leq C\sigma^{k}$ and $0<1-\tau_{k}\leq\frac{(1-\tau_{0})}{\tau_{0}}\sigma^{k}$. Consequently, both sequences $\left\{\|x^{k}-x^{\ast}_{\tau_{k}}\|_{2}\right\}$ and $\left\{1-\tau_{k}\right\}$ globally converge to zero at a linear rate.

The sequence $\left\{x^{k}\right\}$ also satisfies $\|x^{k}-x^{\star}\|_{2}\leq\hat{C}\sigma^{k}$, where $\hat{C}:=C+\frac{(1-\tau_{0})\omega(L_{g}+\|\xi^{0}\|_{2})}{\tau_{0}^{3}\mu_{f}}$. Hence, $\left\{x^{k}\right\}$ globally converges to a solution $x^{\star}$ of (1) at a linear rate.

Let us make some remarks on the result of Theorem 4.1. First, the condition $\omega<1$ is equivalent to $(L-2\mu_{f})m+L_{f}^{2}<m^{2}$. If we choose $m=L>0$, then we have $L_{f}^{2}<2\mu_{f}L$ which leads to $L>\tfrac{L_{f}^{2}}{2\mu_{f}}$. In this case, if we define $\gamma:=\frac{2}{L+m}=\frac{1}{L}$, then (12) becomes

[TABLE]

which reduces to a homotopy proximal gradient method.

To optimize the contraction factor, we need to minimize $1-2\mu_{f}t+L_{f}^{2}t^{2}$ over $t$. This gives us $t=\frac{\mu_{f}}{L_{f}^{2}}$ showing that $m=L=\frac{L_{f}^{2}}{\mu_{f}}$. Hence, we must choose $H_{k}=\tfrac{L_{f}^{2}}{\mu_{f}}\mathbb{I}$, and we obtain $\omega=1-\tfrac{\mu_{f}}{L_{f}}$. Another simple choice of $H_{k}$ is H_{k}=\big{(}\frac{L+m}{2}\big{)}\mathbb{I}.

Next, note that the convergence rate of $\left\{\|x^{k}-x^{\star}\|_{2}\right\}$ in (12) is slower than in the standard proximal variable-metric method. Its contraction factor is $\sigma$ defined in (18). However, (12) possesses some computational advantages that the standard proximal variable-metric method does not have as we will discuss in Section 7.

The linear convergence rate under Assumption 2 is known from the literature for both gradient and Newton-type methods. Nevertheless, our method is new, which works on the parameterized function $f_{\tau}$ instead of $f$. Our method also allows us to flexibly choose the variable matrix $H_{k}$ as long as it satisfies the condition of Theorem 4.1. Another appropriate choice of $H_{k}$ is a rank-one update as proposed in [6].

Remark 4.2* (Nesterov’s accelerated variant).*

Note that we can develop Nesterov’s accelerated variant for (12) under Assumption 2. In this case, the convergence factor in Theorem 4.1 will be improved from $1-\frac{\mu_{f}}{L_{f}}$ to $1-\sqrt{\frac{\mu_{f}}{L_{f}}}$. However, we skip this modification in this paper.

4.2 Linear convergence for self-concordant function $f$ without barrier

We consider the second case where $f$ and g satisfy the following assumption.

Assumption ​​ 3**.**

*The function $f$ in (1) is standard self-concordant $($see Definition A.1$)$. The function $g$ is $L_{g}$-Lipschitz continuous w.r.t. the local norm $\|\cdot\|_{x}$ defined by $f$ with a Lipschitz constant $L_{g}\in[0,+\infty)$ in $\mathrm{dom}(f)\cap\mathrm{dom}(g)$. *

Note that, in Assumption 3, we only require $f$ to be self-concordant, but not necessary a self-concordant barrier. The class of self-concordant functions is much larger than the class of self-concordant barriers. As indicated in Proposition A.2, any generalized self-concordant and strongly convex function is self-concordant. In particular, it covers a few representative applications presented in the introduction. For other examples, we refer the reader to [56, 61].

Under Assumption 3, Algorithm 1 has two stages: Stage 1 finds an initial point, and Stage 2 performs a homotopy scheme. We can skip Stage 3. For any initial value $\tau_{0}\in(0,1)$ of $\tau$, let us choose $\sigma\in(0.318642,1)$ that solves the following inequation:

[TABLE]

Note that there always exists $\sigma\in(0.318642,1]$ that solves (20). Theorem 4.3 below states the convergence of (14) under the Assumption 3. Its proof can be found in Appendix B.2.

Theorem 4.3*.*

Suppose that Assumption 3 holds for (1). Let $\tau_{0}\in(0,1)$ and $\sigma\in(0,1]$ be two constants satisfying (20) and $\left\{(x^{k},\tau_{k})\right\}$ be the sequence generated by (14). Moreover, $x^{0}$ and $\tau_{0}\in(0,1)$ are chosen such that $\lambda_{0}:=\|x^{0}-x^{\ast}_{\tau_{0}}\|_{x^{\ast}_{\tau_{0}}}\leq\beta$ with $\beta:=0.05$. Let us choose $0\leq\delta_{k}\leq\frac{\lambda_{k}}{113}$, and update the parameter $\tau_{k}$ as

[TABLE]

Then, $\|x^{k}-x^{\ast}_{\tau_{k}}\|_{x^{\ast}_{\tau_{k}}}\leq\beta\sigma^{k}$ for $k\geq 0$ and $0<1-\tau_{k}\leq\frac{(1-\tau_{0})}{\tau_{0}}\sqrt{\sigma}^{k}$. Therefore, the sequences $\{\|x^{k}-x^{\ast}_{\tau_{k}}\|_{x^{\ast}_{\tau_{k}}}\}$ and $\left\{1-\tau_{k}\right\}$ both globally converge to zero at a linear rate.

Moreover, there exists $\hat{C}>0$ such that $\|x^{k}-x^{\star}\|_{x^{\star}}\leq\hat{C}\sqrt{\sigma}^{k}$ for all $k\geq 0$. Hence, the sequence $\left\{\|x^{k}-x^{\star}\|_{x^{\star}}\right\}$ also globally converges to zero at a linear rate.

The following result is a direct consequence of Theorem 4.3 and Lemma 2.2 when $g$ is $\bar{L}_{g}$-Lipschitz continuous in $\ell_{2}$-norm, and $f$ is strongly convex and generalized self-concordant.

Corollary 4.4*.*

Assume that $f$ is generalized self-concordant as defined in Definition A.1 and $g$ is Lipschitz continuous with a Lipschitz constant $\bar{L}_{g}\geq 0$ in $\ell_{2}$-norm instead of $g$ being $L_{g}$-Lipschitz continuous w.r.t. $\left\|\cdot\right\|_{x}$. Assume additionally that $f$ is strongly convex with a strong convexity parameter $\mu_{f}>0$. Then, the conclusion of Theorem 4.3 still holds with $L_{g}:=\frac{\bar{L}_{g}}{\sqrt{\mu_{f}}}$.

Remark 4.5*.*

(a) The $L_{g}$-Lipschitz continuity of $g$ w.r.t. a local norm $\left\|\cdot\right\|_{x}$ in Assumption 3 can be replaced by assuming that $\|\nabla{g}(x)\|_{x}^{\ast}\leq L_{g}$ for some $\nabla{g}(x)\in\partial{g}(x)$ for any $x\in\mathrm{dom}(f)$ and $\xi^{0}=\boldsymbol{0}^{p}\in\partial{g}(x^{0})$. By Lemma 2.2, we can easily see that this condition is weaker than the $L_{g}$-Lipschitz continuity of $g$ w.r.t. $\left\|\cdot\right\|_{x}$. For example, if $f(x)=-\ln(x)$, and $g(x)=x^{2}$, then $\|\nabla{g}(x)\|_{x}^{2}=1$ for all $x>0$. In this case, the conclusions of Theorem 4.3 still hold.

(b) Observe that in (21), if the rate of change from $\Delta_{k}$ to $\Delta_{k+1}$ is slow so that $\Delta_{k+1}\approx\Delta_{k}$, then the rate of increment from $\tau_{k}$ to $\tau_{k+1}$ will become faster when $k$ increases.

4.3 Linear convergence under the self-concordant barrier of $f$

When $f$ is a self-concordant barrier, we use a different analysis, and no longer require $g$ to be Lipschitz continuous as stated in the following assumption:

Assumption ​​ 4**.**

*The function $f$ is a $\nu_{f}$-self-concordant barrier as defined in Definition A.1, and $g$ is proper, closed, and convex. For a given $x^{0}\in\mathrm{dom}(F)$, either the analytic center $x^{\star}_{f}$ of $f$ defined by (50) on the interior of the level set $\mathcal{L}_{F}(x^{0}):=\left\{x\in{\mathrm{dom}(F)\mid F(x)}\leq F(x^{0})\right\}$ exists or $\xi^{0}=\boldsymbol{0}^{p}\in\partial{g}(x^{0})$. *

Under Assumption 4, Algorithm 1 also requires Stage 1 and Stage 2, while we can skip Stage 3. For any $\tau_{0}\in(0,1)$, we choose $\sigma\in(0.318642,1]$ such that

[TABLE]

Here, $\bar{c}_{0}:=\theta_{f}\|\xi^{0}\|_{x^{\star}_{f}}$ ($\theta_{f}$ is defined in Appendix A) if $x^{\star}_{f}$ defined by (50) exists, and $\bar{c}_{0}:=0$, otherwise. These two constants $\tau_{0}$ and $C_{0}$ always exist, and $C_{0}\leq\frac{9}{20\sqrt{10}}\approx 0.14230$. The following theorem states the convergence of Algorithm 1 under the self-concordant barrier assumption on $f$, whose proof can be found in Appendix B.3.

Theorem 4.6*.*

Let us choose $\sigma\in(0.318642,1]$ such that (22) holds. Let $\left\{(x^{k},\tau_{k})\right\}$ be the sequence generated by Algorithm 1 using (14) with $\tau_{0}\in(0,1)$ such that $\|x^{0}-x^{\ast}_{\tau_{0}}\|_{x^{\ast}_{\tau_{0}}}\leq\beta$ with $\beta:=0.05$. Let us choose $0\leq\delta_{k}\leq\frac{\lambda_{k}}{113}$, and update $\tau_{k}$ as

[TABLE]

Then, $\|x^{k}-x^{\ast}_{\tau_{k}}\|_{x^{\ast}_{\tau_{k}}}\leq\beta\sigma^{k}$ and $1-\tau_{k}\leq\frac{C_{0}\sqrt{\sigma}^{k}}{(1~{}+~{}C_{0})(\sqrt{\nu_{f}}+\bar{c}_{0})(1-\sqrt{\sigma})}$. Therefore, the sequences $\{\|x^{k}-x^{\ast}_{\tau_{k}}\|_{x^{\ast}_{\tau_{k}}}\}$ and $\left\{1-\tau_{k}\right\}$ both globally converge to zero at a linear rate.

Moreover, if we choose $\sigma\in(0.318642,1]$ such that (22) holds and $\widetilde{C}_{1}:=\frac{C_{0}}{\tau_{0}(1+C_{0})(1-\sqrt{\sigma})-C_{0}}>0$ $($always exist$)$, then $\|x^{k}-x^{\star}\|_{x^{\star}}\leq\frac{\beta\sigma^{k}}{1-\widetilde{C}_{1}\sqrt{\sigma}^{k}}+\widetilde{C}_{1}\sqrt{\sigma}^{k}$. Hence, the sequence $\left\{x^{k}\right\}$ converges to an optimal solution $x^{\star}$ of (1) at a linear rate.

Theorem 4.6 shows the global linear convergence of our inexact proximal-Newton method for solving (1) under Assumption 4. Here, the constants $\sigma$ and $\beta$ balance between the contraction factor and the step-size of the homotopy parameter $\tau$. The choice of $\sigma$ from (22) is conservative due to several rough estimates in our proof. In practice, $\sigma$ can be chosen to be much smaller than one as observed in our numerical experiments.

5 Stage 1: Finding an appropriate initial point

While the variant of Algorithm 1 in Theorem 4.1 can start from any initial point $x^{0}\in\mathrm{dom}(F)$, the variants in Theorem 4.3 and Theorem 4.6 require an appropriate initial point $x^{0}$. More precisely, we need to choose an initial value $\tau_{0}\in(0,1)$ and $x^{0}$ such that $\|x^{0}-x^{\ast}_{\tau_{0}}\|_{x^{\ast}_{\tau_{0}}}\leq\beta$ for a given $\beta:=0.05$. We consider two cases: $g$ is $\mu_{g}$-strongly convex and $g$ is non-strongly convex.

5.1 Inexact damped-step proximal-Newton scheme

We can apply the following inexact damped-step proximal-Newton scheme proposed in [60] to find $x^{0}$. Let us start from any initial point $\hat{x}^{0}\in\mathrm{dom}(F)$, compute a subgradient $\hat{\xi}^{0}\in\partial{g}(\hat{x}^{0})$, and update:

[TABLE]

Here, $0\leq\hat{\delta}_{j}<\hat{\zeta}_{j}$ is the accuracy level defined as in Definition 3.2 and we use the “hat” notation for the iterates to distinguish this procedure from Algorithm 1. The following proposition provides an estimation on the number of iterations needed to find the initial point $x^{0}$, whose proof can be found in [60, Lemma 4.3.].

Proposition 5.1*.*

Let $\left\{\hat{x}^{j}\right\}$ be generated by (24) with $\hat{\delta}_{j}:=\frac{\hat{\zeta}_{j}}{10}$, then after at most

[TABLE]

iterations, we obtain $\hat{x}^{j_{\max}}$ such that $\|\hat{x}^{j_{\max}}-x^{\ast}_{\tau_{0}}\|_{x^{\ast}_{\tau_{0}}}\leq\beta$, where $F_{\tau_{0}}(x):=f(x)-(\frac{1}{\tau_{0}}-1)\langle\hat{\xi}^{0},x\rangle+\tfrac{1}{\tau_{0}}g(x)$, and $\omega(t):=t-\ln(1+t)$.

Proposition 5.1 suggests that we can perform a finite number of damped-step proximal-Newton scheme (24) to find $x^{0}:=\hat{x}^{j_{\max}}$ such that $\|x^{0}-x^{\ast}_{\tau_{0}}\|_{x^{\ast}_{\tau_{0}}}\leq\beta$. Hence, $x^{0}$ is an initial point that satisfies the conditions of Theorem 4.3 and Theorem 4.6.

5.2 Strong convexity of $g$

If $g$ is strongly convex with a strong convexity parameter $\mu_{g}>0$, and $x^{0}$ is not an optimal solution of (1), then we can choose $\tau_{0}$ as

[TABLE]

Here, $\lambda_{\max}(\nabla^{2}{f}(x^{0}))$ is the maximum eigenvalue of $\nabla^{2}{f}(x^{0})$. We will show in Appendix B.4 that $x^{0}$ satisfies $\|x^{0}-x^{\ast}_{\tau_{0}}\|_{x^{\ast}_{\tau_{0}}}\leq\beta$. Hence, Algorithm 1 can start from an arbitrary point $x^{0}\in\mathrm{dom}(F)$.

5.3 Non-strong convexity of $g$ and strong convexity of $f$

We adopt our recent idea in [62] to develop a homotopy scheme to find this initial point $x^{0}$ in a finite number of iterations.

Starting from any $\hat{x}^{0}\in\mathrm{dom}(F)$, we consider the following auxiliary optimality condition depending on a new homotopy parameter $t>0$ and a fixed value $\tau_{0}\in(0,1)$:

[TABLE]

for any $\hat{\xi}^{0}\in\partial{g}(\hat{x}^{0})$. Clearly, when $t=0$, (26) reduces to (9) at $\tau=\tau_{0}$ and $x^{0}\equiv\hat{x}^{0}$. When $t=1$, (26) becomes $0\in\nabla{f}(x^{\ast}_{0})-\nabla{f}(\hat{x}^{0})+\frac{1}{\tau_{0}}(\partial{g}(x^{\ast}_{0})-\hat{\xi}^{0})$, which shows that ${x^{*}_{0}=}\hat{x}^{0}$ is a solution of (26). By applying the homotopy method starting from $t_{0}\approx 1$, and decreases $t_{j}$ to zero, we obtain an approximation $\hat{x}^{j}$ to $x^{\ast}_{\tau_{0}}$. The main step of this scheme is given as follows:

[TABLE]

where the approximation “$:\approx$” is defined as in Definition 3.2, and $t_{0}>0$ is a starting value of $t$. We also use the “hat” notation for the iterates to distinguish this procedure from Algorithm 1.

This scheme is slightly different from (14) with the additional term $-t_{j+1}(\nabla{f}(\hat{x}^{0})+\hat{\xi}^{0})$. The following theorem shows us how to choose $t_{0}$ and update $t$ to guarantee $\|\hat{x}^{j}-x^{\ast}_{\tau_{0}}\|_{x^{\ast}_{\tau_{0}}}\leq\beta$, whose proof is given in Appendix B.5.

Theorem 5.2*.*

Assume that $f$ is self-concordant and $\mu_{f}$-strongly convex with $\mu_{f}>0$. For any given $\beta\in(0,0.05]$, we defined $\Theta:=\sqrt{\frac{99\beta}{500}}-\frac{10\beta}{9}$ $>0.$ Let $\hat{x}^{0}\in\mathrm{dom}(F)$ be an arbitrary starting point, $\hat{\xi}^{0}\in\partial{g}(\hat{x}^{0})$, and $t_{0}$ be chosen such that

[TABLE]

Let $\left\{(\hat{x}^{j},t_{j})\right\}$ be the sequence generated by (27) starting from this $\hat{x}^{0}$ and $t_{0}$. Support further that $t_{j}$ is updated by t_{j+1}:=\big{[}t_{j}-\frac{\Theta}{L_{g}(1+\Theta)}\big{]}_{+} and $\delta_{j}$ satisfies $\hat{\delta}_{j}\leq\frac{\lambda_{j}}{113}$. Then, after at most $j_{\max}:=\left\lfloor\frac{t_{0}M_{0}(1+\Theta)}{\Theta}\right\rfloor$ iterations with $M_{0}:=\frac{\|\nabla{f}(\hat{x}^{0})+\hat{\xi}^{0}\|_{2}}{\sqrt{\mu_{f}}}$, we have $t_{j_{\max}}=0$, and $\|\hat{x}^{j_{\max}}-x^{\ast}_{\tau_{0}}\|_{x^{\ast}_{\tau_{0}}}\leq\beta$.

Theorem 5.2 shows that to find an initial point $x^{0}:=\hat{x}^{j_{\max}}$ for Algorithm 1 such that $\|x^{0}-x^{\ast}_{\tau_{0}}\|_{x^{\ast}_{\tau_{0}}}\leq\beta$, we only need a finite number of iterations $j_{\max}$ as defined in Theorem 5.2. Moreover, in this case, we can take $\xi^{0}:=\hat{\xi}^{0}$ in Algorithm 1.

5.4 Implementation remarks for Algorithm 1

Theoretically, the variants of Algorithm 1 stated in Theorem 4.3 and Theorem 4.6 require a good starting point $x^{0}$ such that $\|x^{0}-x^{\ast}_{\tau_{0}}\|_{x^{\ast}_{\tau_{0}}}\leq\beta$. To find this point, we can use either (24) or (27). However, since we know that when $\tau_{0}=0$, $x^{\ast}_{\tau_{0}}\equiv x^{\ast}_{0}=x^{0}$, in practice we can choose $\tau_{0}>0$ to be sufficiently small such that $x^{\ast}_{\tau_{0}}\approx x^{0}$, and skip Stage 1.

Practically, we only perform two stages as follows:

• •

Skip Stage 1 and choose $\tau_{0}>0$ sufficiently small such that $\|x^{0}-x^{\ast}_{\tau_{0}}\|_{x^{\ast}_{\tau_{0}}}$ is small.

• •

In Stage 2, we choose $\sigma=1$ to guarantee that $\|x^{k}-x^{\ast}_{\tau_{k}}\|_{x^{\ast}_{\tau_{k}}}\leq\beta$ instead of $\|x^{k}-x^{\ast}_{\tau_{k}}\|_{x^{\ast}_{\tau_{k}}}\leq\beta\sigma^{k}$. Then we update $\tau_{k}$ from $\tau_{0}$ to $\tau_{k}\approx 1$.

• •

In Stage 3, we fix $\tau_{k}$ and perform a couple of iterations to reach $\|x^{k}-x^{\ast}_{\tau_{k}}\|_{x^{\ast}_{\tau_{k}}}\leq\varepsilon$.

We only perform Stage 3 if we choose $\sigma=1$. In this case, we only have $\|x^{k}-x^{\ast}_{\tau_{k}}\|_{x^{\ast}_{\tau_{k}}}\leq\beta$. To achieve $\|x^{k}-x^{\ast}_{\tau_{k}}\|_{x^{\ast}_{\tau_{k}}}\leq\varepsilon$, we need to perform a few proximal-Newton iterations with fixed $\tau_{k}$.

6 Primal-Dual-Primal Method

Our second idea is a primal-dual-primal approach to solve (1). We propose a primal-dual-primal method which consists of the following steps:

• •

Construct the Fenchel dual problem (29) of (1).

• •

Apply Algorithm 1 to solve the dual problem (29).

• •

Instead of solving the dual subproblem (13), we dualize it to go back to the primal space.

• •

Construct an approximate primal solution of (1) from its dual approximate solution.

We will show in Section 6.6 that this approach is useful for the well-known model (2). Now, we present this method in detail as follows.

6.1 The dual problem

We assume that $g(x):=\psi(Dx)$, where $\psi$ is a proper, closed, and convex function from $\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$, and $D:\mathbb{R}^{p}\to\mathbb{R}^{n}$ is a linear operator such that $n\leq p$. The dual problem of (1) in this case becomes

[TABLE]

where $f^{\ast}$ and $\psi^{\ast}$ are the Fenchel conjugates of $f$ and $\psi$, respectively.

Let us define $\varphi(y):=f^{\ast}(-D^{\top}y)$. Then, we can compute the gradient and Hessian of $\varphi$ as

[TABLE]

We impose the following assumption.

Assumption ​​ 5**.**

*The function $f$ in (1) is self-concordant as defined in Definition A.1, and $g(x):=\psi(Dx)$, where $\psi:\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$ is a proper, closed, and convex function and $D:\mathbb{R}^{p}\to\mathbb{R}^{n}$ is a linear operator such that $n\leq p$. In addition, $D$ has full-row rank. *

Under Assumption 5, the function $\varphi$ is still a self-concordant function as stated in [44, Theorem 2.4.1] for $\mathrm{dom}(\varphi)$ defined as

[TABLE]

We define the local norm with respect to $\varphi$ as $\|u\|_{y}:=(u^{\top}\nabla^{2}{\varphi}(y)u)^{1/2}$ and its dual norm $\|v\|_{y}^{\ast}:=(v^{\top}\nabla^{2}{\varphi}(y)^{-1}v)^{1/2}$. The optimality condition of the dual problem (29) becomes

[TABLE]

which is necessary and sufficient for $y^{\star}$ to be an optimal solution of (29) if $\mathrm{dom}(\varphi)\cap{\mathrm{dom}(\psi^{*})}\neq\emptyset$.

Let $y^{\star}$ be an optimal solution of (29). Then, from (31), if we define

[TABLE]

then $-D^{\top}y^{\star}\in\partial{f}(x^{\star})$, which leads to $0\in D^{\top}y^{\star}+\partial{f}(x^{\star})$. On the other hand, we have $Dx^{\star}\in\partial{\psi^{\ast}}(y^{\star})$, which leads to $y^{\star}\in\partial{\psi}(Dx^{\star})$. Combining both expressions, we have $0\in D^{\top}\partial{\psi}(Dx^{\star})+\partial{f}(x^{\star})$. Therefore, $x^{\star}$ given by (32) is an exact solution of the primal problem (1).

6.2 The homotopy proximal Newton method methods for the dual problem

To fulfill the assumptions of Theorems 4.3 and 4.6, we assume that one of the following conditions holds:

• •

$f$ satisfied Assumption 5 and $\psi^{\ast}$ is $L_{\psi^{\ast}}$-Lipschitz continuous w.r.t. $\|\cdot\|_{y}$ defined by $\varphi$.

• •

$f$ satisfied Assumption 5 and is $\nu_{f}$-self-concordant barrier.

One can show that $\psi^{\ast}$ is $L_{\psi^{\ast}}$-Lipschitz continuous if $\mathrm{dom}(\partial{\psi})$ is bounded w.r.t. the local norm defined by $f$, i.e. there exists $L_{\psi^{\ast}}>0$ such that $\|u\|_{x}\leq L_{\psi^{\ast}}$ for any $u\in\mathrm{dom}(\partial{\psi})$. Since the dual problem (29) has the same property as the primal one (1) under the above assumptions, let us apply Algorithm 1 with $H_{k}:=\nabla^{2}{\varphi}(y^{k})$ to solve this problem, which leads to

[TABLE]

where \nabla{\varphi_{\tau_{k+1}}}(y^{k}):=\nabla{\varphi}(y^{k})-\big{(}\frac{1}{\tau_{k+1}}-1\big{)}\xi^{0} with $\xi^{0}\in\partial{\psi^{\ast}}(y^{0})$. Here, $y^{k+1}$ is an approximation to the true solution $\bar{y}^{k+1}$ as defined in Definition 3.2, where $\bar{y}^{k+1}$ is given by

[TABLE]

and both the gradient mapping $\nabla{\varphi}$ and the Hessian mapping $\nabla^{2}{\varphi}$ of $\varphi$ are given in (30), respectively. This problem in general does not have a closed form solution. But observe that (34) is a convex composite quadratic programming problem for which highly advanced algorithms such as the semismooth Newton augmented Lagrangian method developed in [37, 69] can be designed to solve it efficiently, as we shall demonstrate later in the numerical experiments.

6.3 The dualization of the subproblem (34)

Instead of solving the dual subproblem (34) directly, we dualize it to obtain the following subproblem in the primal space of $Dx$:

[TABLE]

where

[TABLE]

Clearly, this problem is again a composite strongly convex quadratic program of the same form as (34), but in the primal space of $Dx$. Specially, if $D=\mathbb{I}$, the identity matrix, then (35) is in the primal space of $x$ as in (13).

6.4 Solution reconstruction for (34)

Recall that $\bar{z}^{k+1}$ denotes the exact solution of (35), then we can construct

[TABLE]

as an exact solution of (34).

Assume that we can only solve (35) up to a given accuracy $\delta\geq 0$. In this case, we say that $z^{k+1}$ is a $\delta$-approximate solution to $\bar{z}^{k+1}$ of (35) if for any $\tilde{e}_{k}$ such that $\left\|\tilde{e}_{k}\right\|_{y^{k}}\leq\delta$, we have

[TABLE]

To guarantee (37), we can apply inexact first-order methods to solve (35), see, e.g., in [52, 63].

If $z^{k+1}$ satisfies (37), then we can construct an approximate solution $y^{k+1}$ to $\bar{y}^{k+1}$ as

[TABLE]

The following lemma shows a relation between $z^{k+1}$ of (35) and the approximate solution $y^{k+1}$ of (34), whose proof is given in Appendix B.6.

Lemma 6.1*.*

Let $z^{k+1}$ be a $\delta$-approximate solution to $\bar{z}^{k+1}$ of (35) in the sense of (37). Then, $y^{k+1}$ constructed by (38) is also a $\delta$-approximate solution to the true solution $\bar{y}^{k+1}$ of (34) such that $\mathcal{P}_{k}(y^{k+1})-\mathcal{P}_{k}(\bar{y}^{k+1})\leq\frac{\delta^{2}}{2}$.

6.5 Primal solution recovery

Finally, we show how to recover an approximate primal solution $x^{k}$ of the original problem (1) from its dual approximate solution $y^{k}$. Based on (32), we show below that for an approximate solution $y^{k}$ to $y^{\star}$, the following point

[TABLE]

is an approximate solution to the true solution $x^{\star}$ of (1) as stated in the following theorem whose proof is given in Appendix B.7. In particular, if $D$ is an invertible matrix, then one can show that we can construct an approximate solution $x^{k+1}$ to $x^{\star}$ of (1) from $z^{k+1}$.

Theorem 6.2*.*

Let $y^{\star}$ be an exact solution of the dual problem (29). Then

• $(\mathrm{a})$

$x^{\star}$ constructed by (32) is an exact solution of (1).

• $(\mathrm{b})$

Let $\left\{y^{k}\right\}$ be computed by (38) and $\left\{x^{k}\right\}$ be given by (39) such that $\|y^{k}-y^{\star}\|_{y^{\star}}\!<\!1$. Then

[TABLE]

Consequently, under the conditions of Theorem 4.3 or Theorem 4.6, the sequence $\left\{x^{k}\right\}$ converges linearly to the optimal solution $x^{\star}$ of (1).

• $(\mathrm{c})$

Let $z^{k+1}$ be an approximate solution of (35). If $\|y^{k}-y^{\star}\|_{y^{k}}<1$, then

[TABLE]

Assume that we apply Algorithm 1 to solve the dual problem (29) under the assumptions of Theorem 4.3 or Theorem 4.6 and the choice $\delta_{k}\leq\frac{\lambda_{k}}{113}$. If, in addition, $D$ is invertible, then $x^{k+1}:=D^{-1}z^{k+1}$ is an approximate solution to $x^{\star}$ of (1). Moreover, \big{\{}\|x^{k+1}-x^{\star}\|^{\ast}_{y^{k}}\big{\}} converges linearly to zero.

From Theorem 6.2, we can see that if $D$ is invertible, then we can directly use $x^{k+1}:=D^{-1}z^{k+1}$ to approximate the solution $x^{\star}$ of (1). Otherwise, we can construct an approximate solution $x^{k}$ to $x^{\star}$ by using (39), which requires one evaluation of $\nabla{f^{\ast}}$.

6.6 Applications to covariance estimation

In this section, we apply Algorithm 1 and the primal-dual-primal method in Section 6 to solve the regularized covariance estimation problem (2) as in [18] and its least-squares extension in [32].

We recall the primal regularized covariance estimation problem given in (2). Associated with (2), we can also consider its dual form:

[TABLE]

Here, $\psi^{\ast}$ is the Fenchel conjugate of $\psi(X):=g(X)$. This problem again has the same form as (1). Instead of solving the primal problem (2), we apply Algorithm 1 to solve the dual problem (42) and reconstruct a solution of (2) from its dual.

6.6.1 The main steps of the algorithm

Given $Y_{k}$ such that $Y_{k}+\Sigma\succ 0$, we define $X_{k}:=(Y_{k}+\Sigma)^{-1}$. The main step of the algorithm is to solve the following subproblem

[TABLE]

where \widehat{X}_{k}:=X_{k}-\Xi_{k}\equiv X_{k}-\big{(}\tfrac{1}{\tau_{k+1}}-1\big{)}\Xi_{0} for a fixed $\Xi_{0}\in\partial{\psi^{\ast}(Y_{0})}$. As discussed in Section 6, instead of solving (43), we look at its dual form

[TABLE]

where $C_{k}:=2Y_{k}-\Xi_{k}+\Sigma$. Once $Z_{k+1}$ is computed from (44), we can reconstruct $Y_{k+1}$ as follows:

[TABLE]

and compute an inexact Newton decrement

[TABLE]

Finally, when an $\varepsilon$-solution $\widetilde{Y}^{\star}$ of (42) is computed (i.e. $\widetilde{Y}^{\star}:=Y_{k_{\max}}$), we can reconstruct an approximate solution $\widetilde{X}^{\star}$ of the primal problem (2) by taking $\widetilde{X}^{\star}:=(\Sigma+\widetilde{Y}^{\star})^{-1}$. This computation requires the inverse of a symmetric positive definite matrix, which can be done efficiently by Cholesky decomposition. However, as shown in Theorem 6.2, we can use $Z_{k+1}$ computed by (44) to approximate the true solution $X^{\star}$. This allows us to avoid the matrix inversion $(\Sigma+\widetilde{Y}^{\star})^{-1}$.

6.6.2 The algorithm

Putting together these steps, we obtain a new algorithmic variant for solving (2) as presented in Algorithm 2.

Let us highlight some new features of Algorithm 2 as compared to existing methods in the literature, e.g., [18, 25, 26, 59, 61].

• (a)

Firstly, Algorithm 2 deals with a general regularizer compared to [18, 25, 26]. When $g$ is the $\ell_{1}$-norm regularizer, we can apply coordinate descent methods as in [18, 25, 26] for solving (44) to improve its practical performance.

• (b)

Secondly, Algorithm 2 relies on Algorithm 1 to solve the dual problem (42) instead of standard proximal-Newton methods. It has a linear convergence rate compared to the damped-step scheme which only has a sublinear convergence rate as shown in [59, 61].

• (c)

Thirdly, it does not require any linesearch or any additional assumption in our analysis to achieve a linear convergence rate.

• (d)

Fourthly, the whole algorithm does not require any matrix inversion or Cholesky decomposition as long as we can solve the subproblem (44) with a first order method. This is an important feature for designing parallel and distributed variants of Algorithm 2 as compared to [26].

• (e)

Finally, the subproblem (44) works on the original regularizer $g$ instead of the dual problem as in [59], which preserves the structure such as sparsity on the iterates as promoted by the regularizer $g$.

7 Numerical experiments

We provide some numerical experiments to illustrate our theoretical development. Our experiments are implemented in Matlab 2018a running on a Dell Optiplex 9010, 3.4 GHz Intel Core i7-3770 with 16GB 1600 MHz DDR3 memory.

7.1 Lipschitz gradient and strongly convex models

Now we evaluate the performance of the homotopy proximal-Newton scheme (12) by applying it to solve the following logistic regression problem with an elastic-net regularizer:

[TABLE]

where $\mu_{f}>0$ and $\rho>0$ are two regularization parameters, and $(a_{i},y_{i})\in\mathbb{R}^{p}\times\left\{-1,1\right\}$, $i=1,\ldots,n$ is a given dataset. As shown in [70], the elastic-net regularizer helps to remove variable limitation with more freedom than the classical LASSO model, and it can also carter for groups of nonzero variables. Clearly, f(x):=\frac{1}{n}\sum_{i=1}^{n}\log\big{(}1+\exp(-y_{i}a_{i}^{\top}x)\big{)}+{\frac{\mu_{f}}{2}}\|x\|^{2} is $\mu_{f}$-strongly convex, and $L_{f}$-Lipschitz gradient continuous with $L_{f}:={\frac{1}{2n}}\|A\|^{2}+\mu_{f}$, where $A^{\top}=[a_{1},\ldots,a_{n}]\in\mathbb{R}^{p\times n}.$ Moreover, the function $g(x):=\rho\|x\|_{1}$ is $L_{g}$-Lipschitz continuous with $L_{g}:=\rho$. Hence, Assumption 2 of Theorem 4.1 is satisfied.

We implement Algorithm 1 to solve (47) and compare it with homotopy quasi-Newton variant, standard proximal-gradient scheme [4], and the accelerated proximal-gradient method with line-search and restart [4, 5, 55]. These methods are abbreviated as “HomoPN”, “HomoQuasiPN”, “PG”, and “Ls-Rs-APG”, respectively. We test these algorithms on several binary classification datasets a1a, a9a, w1a, w8a, covtype.binary, news20.binary, rcv1.binary and real-sim from [10], and mnist17 and mnist38 from the mnist dataset where we choose the digits $1$ and $7$, and the digits 3 and 8. The details of these dataset set is given in Table 1.

Following [13], we set $\mu_{f}=\frac{1}{n}$. The parameter $\rho$ for $\ell_{1}$-regularization is selected to produce about $10$ percent of nonzeros coefficients. We should mention here that the subpropblem in the model (47) is an elastic-net regularized least-squares problem. For Algorithm 1 to be numerically efficient, it is crucial for us to solve those subproblems efficiently. Fortunately we can adapt the highly efficient semismooth Newton augmented Lagrangian method in [37] to solve the subproblems (12). We terminate the experiments when the relative gaps are less than a given tolerance $\varepsilon=10^{-6}$, based on the KKT system of (47). Moreover, for PG and Ls-Rs-APG methods, we set the maximum number of iterations at $2\times 10^{4}$. If any method does not achieve our desired accuracy after at most $2\times 10^{4}$ iterations, we use ”—” to represent the result. Our final results are reported in Table 2, where iter is the number of iterations, time[s] is the computational time in second, and rgap is the relative gap times $10^{-7}$. We highlight that (47) can be solved effectively by a stochastic method to a modest level of accuracy when the number of data points $n$ is large. However, in our experiments, we only focus on relatively moderate datasets as we are interested in solving the problems accurately to evaluate the performance of Algorithm 1, and ignore the comparison with stochastic methods.

Table 2 shows that since $f$ is strongly convex and has continuous Lipschitz gradient, the standard proximal-gradient method PG can solve some of the problems slightly better than the homotopy quasi-Newton method HomoQuasiPN, particularly for some of the datasets where the number of features is large (e.g., news20 and real-sim). Note that since the problem (47) is strongly convex, both PG and Ls-Rs-APG have linear convergence rate. Ls-Rs-APG outperforms PG in terms of iteration numbers as stated by the theory. But since we adopt the line search scheme, we found that the times taken by these two methods in solving some of the problems are roughly at the same level. However, HomoQuasiPN is promising when the number of features is moderate but the number of observations is large, for instance the datasets mnist17 and mnist38. The reason is that when the number of observations is large, the cost of computing the $H_{k}$ for HomoPN will be more expensive than that in HomoQuasiPN. Hence, with warm start, HomoQuasiPN will only need a few more iterations to converge but with lower computing cost in each iteration. In other situations, our HomoPN outperforms all other methods in terms of iteration numbers and computation time. Furthermore, our homotopy methods can often solve the problems more accurately. To conclude, our homotopy methods are highly efficient for solving a varieties of large-scale datasets in logistic regression.

7.2 Self-concordant barrier models

We illustrate the variant of Algorithm 1 in Theorem 4.6 for solving the following constrained convex problem with a self-concordant barrier function $f$ arising from D-optimal experimental design, see, e.g. [24, 39]:

[TABLE]

where $A_{i}$ for $i=1,\cdots,p$ are $m\times m$ symmetric positive semidefinite matrices. If we define f(x):=-\log\det\Big{(}\sum_{i=1}^{p}x_{i}A_{i}\Big{)}, and $g(x):=\delta_{\Delta_{p}}(x)$, where $\Delta_{p}$ is the standard simplex in $\mathbb{R}^{p}$, then we can reformulate (48) into (1), and $f$ satisfies the assumptions of Theorem 4.6.

We implement Algorithm 1 to solve (48), and compare it with the interior-point method in [39], where their code is available online at http://www.mypolyuweb.hk/ tkpong/OD_final_codes/. To approximately solve the proximal-Newton subproblem (44), we adapt the semi-smooth Newton-CG augmented Lagrangian method in [37, 69] to solve the composite convex QP problem where the nonsmooth term is given by $g(x)=\delta_{\Delta_{p}}(x)$. We also compare our method with the multiplicative method proposed in [24] and the interior-point method implemented in SDPT3-v.4.0 [57]. We abbreviate these solvers by HomoPN, MUL, IP, and SDPT3, respectively. Unlike other well-established IP solvers, SDPT3 allows us to handle directly the log-determinant functions without reformulation or approximation.

We follow the same procedure as in [39] to generate the data and use the implementation of MUL from [39]. More precisely, we consider the following four design spaces:

[TABLE]

where $s_{i}=\frac{3i}{p},r_{i}=\frac{2i}{p}-1$ and $t_{i}=\frac{i}{p}$.

For each design space, we set $A_{i}:=x_{i}x_{i}^{\top}$ for $i=1,\cdots,p$, with $p=10000,50000,100000$ for $\chi_{1},\chi_{2},\chi_{4}$, and $p=10000,40000,90000$ for $\chi_{3}$. In this case, the problem dimension $m$ is $m=4$ in $\chi_{1}$, $\chi_{2}$ and $\chi_{4}$ and $m=5$ in $\chi_{3}$. The performance of these four methods on $9$ problems of different sizes are reported in Table 3, where (#Iterations) and Time[s] denote the number of iterations and computational time taken, respectively, and $F(x^{k})$ is the approximate optimal objective value attained for (48).

In Table 3, the objective value $F(x^{k})$ is rounded off to seven significant digits. We can see that our homotopy method, HomoPN outperforms the multiplicative algorithm MUL in terms of computational time, and achieving much smaller objective values in all the instances. Our HomoPN also outperforms the interior point method (IP) and SDPT3 in terms of time and also gives slightly better objective values in most instances.

To see the performance of our method compared to MUL and IP, in the following test, we extend the dimension of datasets in $\chi_{1}$, $\chi_{2}$, and $\chi_{3}$ as follows:

[TABLE]

We first run MUL up to $10,000$ iterations, and check whether the HomoPN and IP methods can solve the problem in terms of the function value. We terminate the algorithms when the objective value is less than the value procedured by MUL. When the IP method cannot achieve this objective value, we terminate it and record the time and objective value after at most $10,000$ iterations (this is a large number of iterations for interior point methods). The computation results are presented in Table 4. For these larger dimensional problems, SDPT3 fails to solve them to the required accuracy, and we do not add this method to our experiments.

From Table  4, we can see that our algorithm outperforms both the MUL and IP methods in terms of computation time, while the IP method cannot achieve the desired accuracy when the dimension of the matrix variable $m$ is increased. It also illustrates the ability of our method to solve the problems to an intermediate accuracy. The reason why the IP method cannot solve the problem is that when the dimension is large, the Newton system is very ill-condition, and the solver cannot handle the ill-conditioning difficultly satisfactorily.

7.3 Non-Lipschitz gradient models

Next, we illustrate the variant of Algorithm 1 stated in Theorem 4.3 under Assumption 3 through a non-Lipschitz gradient model. We consider the following problem from Poisson regression [27, 29]:

[TABLE]

where $a_{i}\in\mathbb{R}^{p}$, $y_{i}\in\mathbb{R}$ are given, and $\mu_{f}>0$ and $\rho>0$ are given regularization parameters.

Let $f(x):=\frac{1}{n}\sum_{i=1}^{n}\left(y_{i}\exp\left(-\tfrac{a_{i}^{\top}x}{2}\right)+\exp\left(\tfrac{a_{i}^{\top}x}{2}\right)\right)+\frac{\mu_{f}}{2}\|x\|^{2}_{2}$. Note that if $\mu_{f}>0$, then by Proposition A.2 in Appendix A, $f$ is $M_{f}$-self-concordant with $M_{f}:=\frac{\max_{i}\left\{\|a_{i}\|_{2}\right\}}{\sqrt{\mu_{f}}}$. Moreover, it is also $\mu_{f}$-strongly convex, and $g(x):=\rho\|x\|_{1}$ is $L_{g}$-Lipschitz continuous with $L_{g}:=\rho$. Hence, problem (49) satisfies Assumption 3. Consequently, the results of Theorem 4.3 hold for this problem.

We implement Algorithm 1 to solve (49) and compare it with a quasi-Newton method using BFGS with linesearch scheme and a limited-memory quasi-Newton method (L-BFGS) with both linesearch and restarting scheme. We name these three schemes HomoPN, BFGS, and L-BFGS-LS-R, respectively. Similar to the above example, we adapt the semi-smooth Newton-CG augmented Lagrangian method in [37, 69] to approximately solve the convex composite quadratic programming subproblem (13) in HomoPN.

Note that since the gradient $\nabla{f}$ is not Lipschitz continuous, first-order methods such as proximal gradient-type methods are not applicable. Our experiment reveals that first-order methods such as accelerated Barzilai-Borwein step-size proximal gradient algorithms often failed to converge due to the explosion of the estimated Lipschitz constant of $\nabla{f}$.

We test three algorithms on $6$ datasets downloaded from the UCI dataset repository [38] and Kaggle (https://www.kaggle.com/datasets). Such datasets are used to predict the number of communications for a social post (news and facebook), the review score of a hotel or wine (vegas and wine), the number of goals scored in a game season or the number of 911 calls in a period of time. The sizes of these datasets are given in Table 5. Since these are raw datasets, to use them in our model (49), we perform a sequence of data-preprocessing routines to make our feature matrices well scaled by using methods such as one-hot encoding and min-max scaling.

We test three algorithms: HomoPN, BFGS, and L-BFGS-LS-R on the datasets in Table 5. The computation results are reported in Table 6, where $F(x^{k})$ is the objective value of (49). We can see from the results that while all methods we have tested can solve the problem using real datasets (they all reached almost the same optimal function values), the performance of the two quasi-Newton methods are similar and our homotopy method outperforms them both in terms of computational times and number of iterations. Moreover, based on our observation, when the problem is hard (i.e. the linear system is very ill-conditioned), the number of iterations in the quasi-Newton-type methods increases rapidly. However, our homotopy method is relatively stable in terms of the number of iterations, as was indicated in Section 7.4. To conclude, our method is highly efficient for real datasets.

7.4 The initial point independence of Algorithm 2

Our goal in this example is to observe the independence of performance w.r.t. to the initial point $X_{0}$ in our inexact primal-dual-primal homotopy PN scheme, i.e. Algorithm 2, compared to standard proximal-Newton methods for sparse optimization. We demonstrate this observation on the sparse inverse covariance estimation problem covered by (2) which perfectly fits our assumptions.

As mentioned, in theory, we have shown that Algorithm 2 can achieve a global linear convergence rate. In practice, however, this rate may still be slow. We instead update the homotopy parameter $\tau$ in Algorithm 1 by a longer step based on a linesearch such that the new iterate $X_{k}$ remains in $\mathrm{dom}(\phi)$. This trick allows us to go faster from $\tau\approx 0$ to $\tau\approx 1$ within a few iterations instead of using the worst-case factor. Then, we fix the homotopy parameter $\tau_{k}\approx 1$ and apply a few full-step proximal-Newton iterations to compute the desired solution (Stage 3 of Algorithm 1).

We implement Algorithm 2 and compare it with the primal proximal-Newton method in [61]. We abbreviate these methods by HomoPN and Primal PN, respectively. We use a restarting accelerated proximal-gradient algorithm to approximate the proximal Newton direction. We generate two problem instances using the same procedure as in [36] with $p=500$ and $p=1000$, respectively. In order to obtain a desired sparse solution of (2), we choose $\rho:=0.01$. To see the affect of the initial point $X_{0}$ on the methods we generate two different initial points $X_{0}$.

• •

Case 1 (dense): $X_{0}:=\Sigma^{\dagger}+10^{-6}\times\mathbb{I}$, where $\Sigma^{\dagger}$ is the pseudo-inverse of $\Sigma$.

• •

Case 2 (sparse): $X_{0}:=\mathrm{diag}(1./\mathrm{diag}(\Sigma))$ a diagonal matrix.

Figure 1 shows the convergence of HomoPN and Primal PN for these two cases: the sparse initial point $X_{0}$ and the dense initial point $X_{0}$.

Figure 1 shows that our new algorithmic variant (HomoPN) depends weakly on the initial point $X_{0}$, while the primal proximal Newton method (Primal PN) in [61] strongly depends on the choice of $X_{0}$. The convergence rate of HomoPN is divided into two parts: the homotopy part with a linear rate, and the refining part with a quadratic rate. Figure 1 shows the quadratic convergence on the relative objective residual $\frac{\phi(X_{k})-\phi^{\star}}{|\phi^{\star}|}$ of both algorithms when the iterates approach the optimal solution, but Primal PN (dense) requires a large number of iterations to reach its quadratic convergence region.

Acknowledgments

This work is supported in part by the NSF grant, no. DMS-16-2044 (USA).

Appendix A Self-concordant and generalized self-concordant functions

Our central concept is the self-concordance and self-concordant barrier introduced by Nesterov and Nemirovskii [41, 44], and extensions to generalized self-concordance studied in [56].

For a three-time continuously differentiable and convex function $f$, let $\nabla^{3}{f}(x)[u]$ denote the third-order derivative along the direction $u$. We denote this class of functions by $\mathcal{C}^{3}(\mathrm{dom}(f))$.

Definition A.1*.*

For $f:\mathrm{dom}(f)\subseteq\mathbb{R}^{p}\to\mathbb{R}$ in $\mathcal{C}^{3}(\mathrm{dom}(f))$, we say that:

• $\mathrm{(a)}$

$f$ is $M_{f}$-self-concordant ($$M_{f}\geq 0$$) if for any $x\in\mathrm{dom}(f)$ and $u\in\mathbb{R}^{p}$, we have

[TABLE]

If $M_{f}=2$, then we say that $f$ is standard self-concordant.

• $\mathrm{(b)}$

$f$ is $(M_{f},\kappa)$-generalized self-concordant if for any $x\in\mathrm{dom}(f)$ and $u,v\in\mathbb{R}^{p}$, we have

[TABLE]

where if $0\leq\kappa<2$ or $\kappa>3$, then we use the convention $\frac{0}{0}=0$.

• $\mathrm{(c)}$

$f$ is a $\nu_{f}$-self-concordant barrier $(\nu_{f}\geq 1)$ if $f$ is standard self-concordant with $\mathrm{dom}(f)=\mathrm{int}\left(\mathcal{X}\right)$, $f(x)$ tends to $+\infty$ as $x$ approaches the boundary $\partial{\mathcal{X}}$ of $\mathcal{X}$, and

[TABLE]

• $\mathrm{(d)}$

$f$ is a $\nu_{f}$-self-concordant logarithmically homogeneous barrier function of $\mathcal{X}$ if it is self-concordant barrier and $f(\tau x)=f(x)-\nu_{f}\log(\tau)$ for all $x\in\mathrm{int}\left(\mathcal{X}\right)$ and $\tau>0$.

It is clear that affine and convex quadratic functions are standard self-concordant but not self-concordant barrier. The logistic, Poisson, and DWD models discussed in the introduction are generalized self-concordant with $\kappa\in[2,3]$, but not self-concordant. Self-concordant barriers are often associated with convex sets such as cones and convex bodies. Several simple sets are equipped with a self-concordant barrier. For instance, $f(x):=-\sum_{i=1}^{p}\log(x_{i})$ is a $p$-self-concordant barrier of $\mathbb{R}^{p}_{+}$, $f(X):=-\log\det(X)$ is a $p$-self-concordant barrier of $\mathbb{S}^{p}_{+}$, and $f(x,t)=-\log(t^{2}-\langle x,x\rangle)$ is a $2$-self-concordant barrier of the Lorentz cone $\mathcal{L}_{p+1}:=\left\{(x,t)\in\mathbb{R}^{p}\times\mathbb{R}_{+}\mid\left\|x\right\|\leq t\right\}$. The relation between self-concordant and generalized self-concordant is stated in the following proposition.

Proposition A.2*.*

[56] Let $f\in\mathcal{C}^{3}(\mathrm{dom}(f))$ be an $(M_{f},\kappa)$-generalized self-concordant and $\kappa\in(0,3]$ as defined in Definition A.1. Then, if $f$ is strongly convex with a strong convexity parameter $\mu_{f}>0$, then $f$ is $\widehat{M}_{f}$-self-concordant with $\widehat{M}_{f}:=\frac{M_{f}}{(\sqrt{\mu_{f}})^{3-\kappa}}$.

Proposition A.2 shows that the regularized logistic, Poisson, and DWD models in the introduction become self-concordant due to the aid of a regularization term $\frac{\mu_{f}}{2}\|x\|^{2}$.

Let $f$ be a $\nu_{f}$-self-concordant barrier of $\mathcal{X}$. Then

[TABLE]

is referred to as the analytical center of $\mathcal{X}$. If $\mathcal{X}$ is bounded, then $x^{\star}_{f}$ exists and is unique. Otherwise, we often add an artificial bound or consider $f$ on a sublevel set of $F$ in (1) to guarantee the existence of $x_{f}^{\star}$. Let $\theta_{f}:=\nu_{f}+2\sqrt{\nu_{f}}$ be defined for a general self-concordant barrier, and $\theta_{f}:=1$ be defined for a self-concordant logarithmically homogeneous barrier. Then, $\|v\|_{x}^{*}\leq\theta_{f}\|v\|_{x_{f}^{\star}}^{*}$ and $\|x-x^{\star}_{f}\|_{x^{\star}_{f}}\leq\theta_{f}$ for $x\in\mathrm{dom}(f)$ and $v\in\mathbb{R}^{p}$ if $x^{\star}_{f}$ exists.

Let $\mathcal{K}$ be a proper, closed, pointed, and convex cone. If $\mathcal{K}$ is endowed with a $\nu_{f}$-self-concordant logarithmically homogeneous barrier function $f$, then its Legendre transformation [44]:

[TABLE]

is also a $\nu_{f}$-self-concordant logarithmically homogeneous barrier of the anti-dual cone $-\mathcal{K}^{*}$ of $\mathcal{K}$. For instance, if $\mathcal{K}=\mathbb{S}_{+}^{p}$, then $\mathcal{K}^{*}=S_{+}^{p}=\mathcal{K}$ (self-dual cone). A barrier function of $\mathbb{S}_{+}^{p}$ is $f(X):=-\log\det(X)$. Hence, $f^{*}(Y)=-p-\log\det(-Y)$ is a barrier function of $-\mathcal{K}^{*}$.

Appendix B Convergence and iteration-complexity analysis of Algorithm 1

We break up our analysis into several lemmas and theorems. The following lemmas provide several key estimates for our proofs.

Lemma B.1*.*

Let $x^{0}\in\mathrm{dom}(F)$ be given, $\xi^{0}\in\partial{g}(x^{0})$, and $L_{g}$ be the Lipschitz constant of $g$. Let $x^{\ast}_{\tau}$ and $x^{\ast}_{\hat{\tau}}$ be two solutions of (9) at $\tau\in(0,1]$ and $\hat{\tau}\in(0,1]$, respectively. Then the following results hold.

• $\mathrm{(a)}$

If $f$ is $\mu_{f}$-strongly convex, then we have

[TABLE]

• $\mathrm{(b)}$

If $f$ is self-concordant, then we have

[TABLE]

• $\mathrm{(c)}$

Let $x^{\ast}_{t}$ and $x^{\ast}_{\hat{t}}$ be two solutions of (26) at $t\in(0,1]$ and $\hat{t}\in(0,1]$, respectively. If $f$ is self-concordant, then

[TABLE]

• $\mathrm{(d)}$

Let $x^{\ast}_{s}$ be a solution of either (9) or (26). If $f$ is self-concordant, and $x^{0}$ and $s\in(0,1)$ are chosen such that $\|x^{0}-x^{\ast}_{s}\|_{x^{0}}\leq\gamma<1$, then $\|x^{0}-x^{\ast}_{s}\|_{x^{\ast}_{s}}\leq\frac{\gamma}{1-\gamma}$.

• $\mathrm{(e)}$

If $g$ is $\mu_{g}$-strongly convex, and $x^{\ast}_{\tau}$ is a solution of (9), then

[TABLE]

Proof B.2*.*

(a) From (9), we have $-\tau(\nabla{f}(x^{\ast}_{\tau})+\xi^{0})\in\partial{g}(x^{\ast}_{\tau})-\xi^{0}$. Using the monotonicity of $\partial{g}$ and $\xi^{0}\in\partial{g}(x^{0})$, we have $\langle\nabla{f}(x^{\ast}_{\tau})+\xi^{0},x^{\ast}_{\tau}-x^{0}\rangle\leq 0$, which implies

[TABLE]

If $f$ is $\mu_{f}$-strongly convex, then using the Cauchy-Schwarz inequality, the last inequality leads to

[TABLE]

which can be simplified as $\|x^{0}-x^{\ast}_{\tau}\|_{2}\leq\tfrac{\|\nabla{f}(x^{0})+\xi^{0}\|_{2}}{\mu_{f}}$. This is exactly the first estimate of (51).

Using (9) again with $\hat{\tau}$ we have $-\hat{\tau}(\nabla{f}(x^{\ast}_{\hat{\tau}})+\xi^{0})\in\partial{g}(x^{\ast}_{\hat{\tau}})-\xi^{0}$. Combining this expression and $-\tau(\nabla{f}(x^{\ast}_{\tau})+\xi^{0})\in\partial{g}(x^{\ast}_{\tau})-\xi^{0}$, and using the monotonicity of $\partial{g}$, we have

[TABLE]

Rearranging this inequality, we obtain

[TABLE]

If $f$ is $\mu_{f}$-strongly convex, then we have $\langle\nabla{f}(x^{\ast}_{\tau})-\nabla{f}(x^{\ast}_{\hat{\tau}}),x^{\ast}_{\tau}-x^{\ast}_{\hat{\tau}}\rangle\geq\mu_{f}\|x^{\ast}_{\tau}-x^{\ast}_{\hat{\tau}}\|_{2}^{2}$. Combining this estimate with the last inequality, and then using the Cauchy-Schwarz inequality, we obtain the first inequality of the second line of (51).

Using again (9), we have $-\tau(\nabla{f}(x^{\ast}_{\tau})+\xi^{0})\in\partial{g}(x^{\ast}_{\tau})-\xi^{0}$. Since $g$ is $L_{g}$-Lipschitz continuous, the last expression leads to $\tau\|\nabla{f}(x^{\ast}_{\tau})+\xi^{0}\|_{2}\leq L_{g}+\|\xi^{0}\|_{2}$, which implies (51).

(b) If $f$ is self-concordant, then we have

[TABLE]

With a similar proof as in (a), we have

[TABLE]

which is the first line of (52). The second line of (52) is proved similarly as that of (51) and using the fact that $\|\xi^{0}\|_{x^{\ast}_{\tau}}^{\ast}\leq L_{g}$.

(c) The proof of (53) is very similar to the proof of (52) and we omit it here.

(d) Note that if $f$ is self-concordant, then we have $\|x^{0}-x^{\ast}_{s}\|_{x^{\ast}_{s}}\leq\frac{\|x^{0}-x^{\ast}_{s}\|_{x^{0}}}{1-\|x^{0}-x^{\ast}_{s}\|_{x^{0}}}$ due to [41, Theorem 4.1.5] as long as $\|x^{0}-x^{\ast}_{s}\|_{x^{0}}<1$. If $\|x^{0}-x^{\ast}_{s}\|_{x^{0}}\leq\gamma<1$, then the last inequality implies that $\|x^{0}-x^{\ast}_{s}\|_{x^{\ast}_{s}}\leq\frac{\gamma}{1-\gamma}$, which proves (d).

(e) By the choice of $\xi^{0}$, we have $0\in\partial{g}(x^{0})-\xi^{0}$. From (9), we also have

[TABLE]

Using the $\mu_{g}$-strong monotinicity of $\partial{g}$ and the monotonicity of $\nabla{f}$, we have

[TABLE]

By the Cauchy-Schwarz inequality, we obtain $\|x^{0}-x^{\ast}_{\tau}\|_{2}\leq\frac{\tau}{\mu_{g}}\|\nabla{f}(x^{0})+\xi^{0}\|_{2}$, which is the desired result in (54).

Lemma B.3*.*

Let $M>0$ and $q\in(0,1)$ be two given constants. Let $\left\{\tau_{k}\right\}\subset(0,1)$ be a given sequence such that $\tau_{k+1}\leq\tau_{k}+Mq^{k}\tau_{k}\tau_{k+1}$ for all $k\geq 0$. Then

[TABLE]

as long as the denominators are well-defined.

Proof B.4*.*

Define $s_{k}:=\frac{1}{\tau_{k}}$. Then, $\tau_{k+1}\leq\tau_{k}+Mq^{k}\tau_{k}\tau_{k+1}$ is equivalent to $s_{k+1}\geq s_{k}-Mq^{k}$. By induction, we have $s_{k}\geq s_{0}-M\sum_{i=0}^{k-1}q^{i}=s_{0}-\frac{M(1-q^{k})}{1-q}$. Hence, we can show that $\tau_{k}\leq\frac{\tau_{0}(1-q)}{1-q-\tau_{0}M(1-q^{k})}$. This estimate leads to

[TABLE]

which proves (55).

Lemma B.5*.*

Assume that $f$ is $\mu_{f}$-strongly convex and $L_{f}$-smooth. Suppose $0<m\leq L<+\infty$ are given parameters such that $\omega:=\frac{1}{m}\sqrt{(L-2\mu_{f})m+L_{f}^{2}}<1$. Let $\left\{x^{k}\right\}$ be the sequence generated by (12) using $H_{k}\in\mathcal{S}^{p}_{++}$ such that $m\mathbb{I}\preceq H_{k}\preceq L\mathbb{I}$. Then

[TABLE]

Proof B.6*.*

Using $\nabla{f}_{\tau_{k+1}}(\cdot):=\nabla{f}(\cdot)-(\frac{1}{\tau_{k+1}}-1)\xi^{0}$, $H=H_{k}$, and $\tau=\tau_{k+1}$ in (10), we get x^{\ast}_{\tau_{k+1}}=\mathrm{prox}_{\frac{1}{\tau_{k+1}}g}^{H_{k}}\big{(}x^{\ast}_{\tau_{k+1}}-H_{k}^{-1}\nabla{f_{\tau_{k+1}}}(x^{\ast}_{\tau_{k+1}})\big{)}. Combining this expression and (12), then using the non-expansiveness of $\mathrm{prox}_{\frac{1}{\tau_{k+1}}g}^{H_{k}}$ in (5), we can derive that

[TABLE]

Next, from the fact that $m\mathbb{I}\preceq H_{k}\preceq L\mathbb{I}$ and the strong convexity and smoothness of $f$, we have

[TABLE]

Substituting these estimates into (57), we get

[TABLE]

Using again $m\mathbb{I}\preceq H_{k}\preceq L\mathbb{I}$, the last inequality leads to

[TABLE]

This estimate can be simplified as in (56) with \omega:=\frac{1}{m}\big{(}(L-2\mu_{f})m+L_{f}^{2}\big{)}^{1/2}<1.

Lemma B.7*.*

Assume that $f$ is self-concordant. Let $x^{k+1}$ be generated by the inexact scheme (14) up to a given accuracy $\delta_{k}\geq 0$. Let $\lambda_{k+1}:=\|x^{k+1}-x^{\ast}_{\tau_{k+1}}\|_{x^{\ast}_{\tau_{k+1}}}$ and $\hat{\lambda}_{k}:=\|x^{k}-x^{\ast}_{\tau_{k+1}}\|_{x^{\ast}_{\tau_{k+1}}}$ be defined by (17). If $1-4\hat{\lambda}_{k}+2\hat{\lambda}_{k}^{2}>0$, then we have

[TABLE]

Proof B.8*.*

Note that $\nabla{f}_{\tau_{k+1}}(\cdot):=\nabla{f}(\cdot)-(\tfrac{1}{\tau_{k+1}}-1)\xi^{0}$. Using (10) with $H\equiv H^{\ast}_{k}:=\nabla^{2}f(x^{\ast}_{\tau_{k+1}})$ and $\tau=\tau_{k+1}$, we have

[TABLE]

Define $H_{k}:=\nabla^{2}{f}(x^{k})$ and $e_{k}:={H_{k}^{*}}^{-1}\left(H_{k}^{*}-H_{k}\right)(\bar{x}^{k+1}-x^{k})$. Then, we can rewrite (13) as

[TABLE]

Using (5), we can derive from (59) and (60) that

[TABLE]

First, we estimate the following term of (B.8) using definition of $\nabla{f}_{\tau_{k+1}}$ and $H_{k}^{\ast}$ as

[TABLE]

By the self-concordance of $f$, with the same proof as in [41, Theorem 4.1.14], we can derive that

[TABLE]

Next, using [41, Corollary 4.1.4], we estimate $\|e_{k}\|_{x^{\ast}_{\tau_{k+1}}}$ as follows:

[TABLE]

provided that $\|x^{k}-x^{\ast}_{\tau_{k+1}}\|_{x^{\ast}_{\tau_{k+1}}}{\!\!\!}<1$. Combining these two estimates, we can derive from (B.8) that

[TABLE]

Note that $\|\bar{x}^{k+1}-x^{k}\|_{x^{\ast}_{\tau_{k+1}}}\leq\|\bar{x}^{k+1}-x^{\ast}_{\tau_{k+1}}\|_{x^{\ast}_{\tau_{k+1}}}+~{}\|x^{k}-x^{\ast}_{\tau_{k+1}}\|_{x^{\ast}_{\tau_{k+1}}}$ by the triangle inequality. Using this in the last inequality, and rearranging the result, we obtain

[TABLE]

Moreover, using the self-concordance of $f$ in [41, Theorem 4.1.5], we can also derive that

[TABLE]

Combining the last two inequalities and using the definition of $\lambda_{k+1}$, $\hat{\lambda}_{k}$, and $\delta_{k}$, we obtain

[TABLE]

which is exactly (58). Here, the right-hand side is well-defined since $1-4\hat{\lambda}_{k}+2\hat{\lambda}_{k}^{2}>0$.

B.1 Proof of Theorem 4.1: Global linear convergence for smooth and strongly convex function $f$

From (56) of Lemma B.5, we can write

[TABLE]

Using the estimate (51) of Lemma B.1 with $\tau=\tau_{k+1}$ and $\hat{\tau}=\tau_{k}$, we obtain

[TABLE]

Combining these two estimates and using the triangle inequality, we can derive that

[TABLE]

where $\omega:=\frac{1}{m}\sqrt{(L-2\mu_{f})m+L_{f}^{2}}<1$, and $\Delta{\tau_{k}}:=|\tau_{k+1}-\tau_{k}|$. In addition, by using (51) of Lemma B.1 with $\tau=\tau_{0}$, we have

[TABLE]

Now, let us assume that $\|x^{k}-x^{\ast}_{\tau_{k}}\|_{2}\leq C\sigma^{k}$ for some $\sigma\in(\omega,1]$ and $C$ is given in (63). In order to get $\|x^{k+1}-x^{\ast}_{\tau_{k+1}}\|_{2}\leq C\sigma^{k+1}$, from (62), we impose the following condition:

[TABLE]

Since $\tau_{k}<\tau_{k+1}$ and $\sigma\in(\omega,1]$, this inequality is equivalent to

[TABLE]

Define $C_{1}:=\frac{C(\sigma-\omega)\mu_{f}}{\omega(L_{g}+\|\xi^{0}\|_{2})}=\frac{(\sigma-\omega)\|\nabla{f}(x^{0})+\xi^{0}\|_{2}}{\omega(L_{g}+\|\xi^{0}\|_{2})}>0$. Then, the last inequality leads to

[TABLE]

Using Lemma B.3 with $M:=C_{1}$ and $q:=\sigma$, we get

[TABLE]

Let us choose $\tau_{0}$ and $\sigma$ such that $(1-\tau_{0})(1-\sigma)-\tau_{0}C_{1}=0$. Using the formula of $C_{1}$ above, this condition is equivalent to

[TABLE]

Let us fix $\tau_{0}\in(0,1)$. Then (65) shows that we can choose $\sigma:=\frac{1-\tau_{0}+\tau_{0}\omega\Gamma}{1-\tau_{0}+\tau_{0}\Gamma}\in(\omega,1)$ as shown in (18) to guarantee (65), where $\Gamma:=\frac{\|\nabla{f}(x^{0})+\xi^{0}\|_{2}}{\omega(L_{g}+\|\xi^{0}\|_{2})}$.

With the choice of $\tau_{0}$ and $\sigma$, we have $0\leq 1-\tau_{k}\leq\frac{(1-\tau_{0})\sigma^{k}}{\tau_{0}+(1-\tau_{0})\sigma^{k}}$ and $0<\frac{(1-\tau_{0})\sigma^{k}}{\tau_{0}+(1-\tau_{0})\sigma^{k}}<1$. If we update $1-\tau_{k}:=\frac{(1-\tau_{0})\sigma^{k}}{\tau_{0}+(1-\tau_{0})\sigma^{k}}$ as shown in (19), then the condition (64) holds. Moreover, $0\leq 1-\tau_{k}\leq\frac{(1-\tau_{0})\sigma^{k}}{\tau_{0}}$, where the right-hand side $\frac{(1-\tau_{0})\sigma^{k}}{\tau_{0}}$ converges to zero as $k$ tends to $+\infty$.

Finally, since $x^{\star}=x^{\ast}_{1}$, by the triangle inequality, we have

[TABLE]

Hence, we obtain $\|x^{k}-x^{\star}\|_{2}\leq\hat{C}\sigma^{k}$ with $\hat{C}:=C+\frac{\omega(L_{g}+\|\xi^{0}\|_{2})(1-\tau_{0})}{\tau_{0}^{3}\mu_{f}}$. This is exactly the last conclusion of Theorem 4.1. $\square$

B.2 Proof of Theorem 4.3: Linear convergence under Assumption 3

Let $\lambda_{k}$ and $\hat{\lambda}_{k}$ be defined by (17) and let $\Delta_{k}^{\ast}:=\|x^{\ast}_{\tau_{k+1}}-x^{\ast}_{\tau_{k}}\|_{x_{\tau_{k+1}}^{\ast}}$ be the local distance between the true solutions $x^{\ast}_{\tau_{k+1}}$ and $x^{\ast}_{\tau_{k}}$ of (9). We divide the proof of Theorem 4.3 into the following steps:

Step 1: Upper bound on $\lambda_{k+1}$: Recall (58) from Lemma B.7 as follows:

[TABLE]

By the triangle inequality and [41, Theorem 4.1.5], we can also derive that

[TABLE]

Now, consider the function $\psi(t):=\frac{3-2t}{1-4t+2t^{2}}$ on $[0,0.1145]$. It is straightforward to numerically check that $3\leq\psi(t)\leq 5$. Hence, we can overestimate (66) as follows:

[TABLE]

Using the fact that $\hat{\lambda}_{k}\leq\frac{\lambda_{k}}{1-\Delta_{k}^{\ast}}+\Delta_{k}^{\ast}$ above, we overestimate again this inequality as

[TABLE]

where we require $\frac{\lambda_{k}}{1-\Delta_{k}^{\ast}}+\Delta_{k}^{\ast}\leq 0.1145$.

Now, if we impose $\lambda_{k}\leq 0.05$, then $\frac{\lambda_{k}}{1-\Delta_{k}^{\ast}}+\Delta_{k}^{\ast}\leq 0.1145$ implies that $\Delta_{k}^{\ast}\leq 0.055$. Therefore, if we choose $\delta_{k}:=\tfrac{\lambda_{k}}{113}$, then we can further simplify the above inequality to get

[TABLE]

as long as $\lambda_{k}\leq 0.05$ and $\Delta_{k}^{\ast}\leq 0.055$.

Step 2: Upper bound on $\Delta^{\ast}_{k}$: Let us assume that $\lambda_{k}\leq 0.05\sigma^{k}$ for some $\sigma\in(0,1]$. In order to guarantee $\lambda_{k+1}\leq 0.05\sigma^{k+1}$, using the last inequality (67), we have to impose the following condition

[TABLE]

Provided that $\sigma>\frac{25}{81}+0.01=0.318642$, the previous condition is equivalent to

[TABLE]

where $C_{k}:=\frac{1}{10}\sqrt{\sigma-0.01}-\frac{1}{18}\sqrt{\sigma}^{k}>0$.

Step 3: Update rule of $\tau_{k}$: Using (52) with $\tau:=\tau_{k+1}$ and $\hat{\tau}:=\tau_{k}$, we obtain

[TABLE]

If we update $\tau_{k}$ as (21), then we can see that

[TABLE]

This guarantees (68).

Step 4: Conditions on $\sigma$ and $\tau_{0}$: Define $\bar{C}_{0}:=\frac{2L_{g}(1-\tau_{0})(1-\sqrt{\sigma})}{\tau_{0}}$ and choose $\sigma\in(0.318642,1]$ such that

[TABLE]

Note that the above condition holds if

[TABLE]

In turn, the condition (72) holds if we choose $\sigma\in(0.318642,1)$ such that $\frac{\bar{C}_{0}}{1-\bar{C}_{0}}\leq\frac{\sqrt{\sigma-0.01}}{10}-\frac{1}{18}$. Using the explicit expression of $\bar{C}_{0}$, we obtain that

[TABLE]

which is exactly the condition (20). It is not hard to show that this inequality always has a solution $\sigma\in(0.318642,1]$. Moreover, since $\Delta_{k}^{\ast}\leq 0.055$, from (68), we also require $C_{k}\leq 0.055$ for all $k\geq 0$. It is easy to check numerically that this condition is always satisfied for $\sigma\in(0.318642,1]$.

Step 5: Bound on $\tau_{k}$: Next, using (71), the definition of $\bar{C}_{0}$, and Lemma B.3 with $M:=\bar{C}_{0}$ and $q:=\sqrt{\sigma}$, similar to the proof of Theorem 4.1, we can show that

[TABLE]

Hence 0\leq 1-\tau_{k}\leq\big{(}\frac{1-\tau_{0}}{\tau_{0}}\big{)}\sqrt{\sigma}^{k}.

Step 6: Linear convergence rate of $\left\{x^{k}\right\}$: Finally, note that $x^{\star}=x^{\ast}_{1}$ and $\tau_{k+1}=1$, we can show from (69) and (73) that

[TABLE]

Since in (72) we have chosen $\bar{C}_{0}$ such that $\bar{C}_{0}\sqrt{\sigma}^{k}<1$, this inequality implies that

[TABLE]

where $C_{1}:=\frac{\bar{C}_{0}}{1-\bar{C}_{0}}>0$. Using the triangle inequality, the self-concordance of $f$, and (74), we can derive that

[TABLE]

where $\hat{C}:=\frac{0.05}{1-C_{1}}+C_{1}=\frac{0.05(1-\bar{C}_{0})}{1-2\bar{C}_{0}}+\frac{\bar{C}_{0}}{1-\bar{C}_{0}}$. This inequality shows that $\left\{x^{k}\right\}$ converges to a solution $x^{\star}$ of (1) at a linear rate. $\square$

B.3 Proof of Theorem 4.6: Linear convergence

when $f$ is a self-concordant barrier

Let $\Delta_{k}^{\ast}:=\|x^{\ast}_{\tau_{k+1}}-x^{\ast}_{\tau_{k}}\|_{x_{\tau_{k+1}}^{\ast}}$ as defined in Theorem 4.1. Using (52) with $\tau:=\tau_{k+1}$ and $\hat{\tau}:=\tau_{k}$, we get

[TABLE]

Since $\|\xi^{0}\|_{x^{\ast}_{\tau_{k+1}}}^{\ast}\leq\bar{c}_{0}:=\theta_{f}\|\xi^{0}\|^{\ast}_{x^{\star}_{f}}$ and $\|\nabla{f}(x^{\ast}_{\tau_{k+1}})\|_{x^{\ast}_{\tau_{k+1}}}^{\ast}\leq\sqrt{\nu_{f}}$ due to the self-concordant barrier property of $f$, the last inequality leads to

[TABLE]

Since we want to guarantee $\|x^{k}-x^{\ast}_{\tau_{k}}\|_{x^{\ast}_{\tau_{k}}}\leq\beta\sigma^{k}$ for $\beta=0.05$ and $\sigma\in(0,1]$ as in Theorem 4.3, similar to the proof of Theorem 4.3, we can choose

[TABLE]

where $C_{k}:=\frac{1}{10}\sqrt{\sigma-0.01}-\frac{1}{18}\sqrt{\sigma}^{k}$. Here, $C_{k}>0$ for $k\geq 0$ if $\sigma\in(0.318642,1]$.

Now, from the update rule (23) of Theorem 4.6, we have

[TABLE]

This inequality shows that (76) automatically holds.

Let $a_{k}:=\frac{C_{k}\sqrt{\sigma}^{k}}{(\sqrt{\nu_{f}}+\bar{c}_{0})(1+C_{k}\sqrt{\sigma}^{k})}\in(0,0.052133]$. It is easily to show that $a_{k}\geq\tfrac{C_{0}\sqrt{\sigma}^{k}}{(\sqrt{\nu_{f}}+\bar{c}_{0})(1+C_{0}\sqrt{\sigma}^{k})}$, where $C_{0}:=\frac{\sqrt{\sigma-0.01}}{10}-\frac{1}{18}>0$. By induction, we get $\tau_{k}=\tau_{0}\prod_{i=0}^{k}(1+a_{i})$. Using an elementary inequality $\prod_{i=0}^{k}(1+a_{i})\geq 1+\sum_{i=0}^{k}a_{i}$, we have

[TABLE]

Hence, for any $\tau_{0}\in(0,1)$, if we choose $\sigma\in(0.318642,1]$ such that \sigma\geq\big{(}1-\frac{\tau_{0}C_{0}}{(1-\tau_{0})(1+C_{0})(\sqrt{\nu_{f}}+\bar{c}_{0})}\big{)}^{2}, then we have $1-\tau_{k}\leq\frac{C_{0}\sqrt{\sigma}^{k}}{(1+C_{0})(\sqrt{\nu_{f}}+\bar{c}_{0})(1-\sqrt{\sigma})}$.

Finally, note that $x^{\star}=x^{\ast}_{1}$, we can show from (75) that

[TABLE]

This shows that

[TABLE]

where ${\widetilde{C}_{1}}:=\frac{C_{0}}{\tau_{0}(1+C_{0})(1-\sqrt{\sigma})-C_{0}}>0$. Next, using the triangle inequality, the self-concordance of $f$, and (77), we can derive

[TABLE]

This inequality shows that $\left\{x^{k}\right\}$ converges to a solution $x^{\star}$ of (1) at a linear rate. $\square$

B.4 Proof of (25): The choice of $\tau_{0}$

From (54), we have $\|x^{0}-x^{\ast}_{\tau_{0}}\|_{2}\leq\frac{\tau_{0}}{\mu_{g}}\|\nabla{f}(x^{0})+\xi^{0}\|_{2}$. Let $\lambda_{\max}(\nabla^{2}{f}(x^{0}))$ be the maximum eigenvalue of $\nabla^{2}{f}(x^{0})$. Hence, if we assume that

[TABLE]

then we have

[TABLE]

Since

[TABLE]

as long as $\|x^{0}-x^{\ast}_{\tau_{0}}\|_{x^{0}}<1$, to guarantee that $\|x^{0}-x^{\ast}_{\tau_{0}}\|_{x^{\ast}_{\tau_{0}}}\leq\beta$, we impose $\frac{\gamma}{1-\gamma}\leq\beta$. The last condition is equivalent to $\gamma\leq\frac{\beta}{1+\beta}$, which shows that

[TABLE]

This condition is exactly (25). $\square$

B.5 Proof of Theorem 5.2: Finding an initial point

Let $\lambda_{j}:=\|x^{j}-x^{\ast}_{t_{j}}\|_{x^{\ast}_{t_{j}}}$ and $\Delta^{\ast}_{j}:=\|x^{\ast}_{t_{j+1}}-x^{\ast}_{t_{j}}\|_{x^{\ast}_{t_{j+1}}}$ Assume that $\lambda_{j}\leq\beta$, from (67), to guarantee $\lambda_{j+1}\leq\beta$, we impose the condition

[TABLE]

This condition holds if $\Delta_{j}^{\ast}\leq\Theta:=\sqrt{\frac{99\beta}{500}}-\frac{10\beta}{9}$. Since $\beta\leq 0.1145<0.16038$, we have $\Theta:=\sqrt{\frac{99\beta}{500}}-\frac{10\beta}{9}>0$. Using (53) with $t:=t_{j+1}$ and $\hat{t}:=t_{j}$, we obtain that

[TABLE]

Since $f$ is $\mu_{f}$-strongly convex, we have $\|\nabla{f}(x^{0})+\xi^{0}\|_{x^{\ast}_{t_{k+1}}}^{\ast}\leq\frac{1}{\sqrt{\mu_{f}}}\|\nabla{f}(\hat{x}^{0})+\hat{\xi}^{0}\|_{2}$. Hence, we get $\frac{\Delta_{j}^{\ast}}{1+\Delta_{j}^{\ast}}\leq(t_{j}-t_{j+1})\frac{\|\nabla{f}(\hat{x}^{0})+\hat{\xi}^{0}\|_{2}}{\sqrt{\mu_{f}}}$.

To guarantee this condition, we impose $\frac{\Theta}{1+\Theta}=M_{0}(t_{j}-t_{j+1})$, where $M_{0}:=\frac{\|\nabla{f}(\hat{x}^{0})+\hat{\xi}^{0}\|_{2}}{\sqrt{\mu_{f}}}$. This allows us to update $t_{j}$ to $t_{j+1}:=t_{j}-\frac{\Theta}{L_{g}(1+\Theta)}$. If we start from $t_{0}\approx 1$, then we have $t_{j}:=t_{0}-\frac{\Theta}{M_{0}(1+\Theta)}(j+1)$. Hence, $t_{j}=0$ if $j+1\geq\frac{t_{0}(1+\Theta_{0})}{\Theta}$. This shows that after $j_{\max}:=\left\lfloor\frac{t_{0}M_{0}(1+\Theta)}{\Theta}\right\rfloor$ iterations, we obtain $\hat{x}^{j_{\max}}$ such that $\|\hat{x}^{j_{\max}}-x^{\ast}_{\tau_{0}}\|_{x^{\ast}_{\tau_{0}}}\leq\beta$.

Finally, we show how to choose $t_{0}$ such that $\|\hat{x}^{0}-x^{\ast}_{t_{0}}\|_{x^{\ast}_{t_{0}}}\leq\beta$. Using (53) with $t:=t_{0}$, if $(1-t_{0})\|\nabla{f}(\hat{x}^{0})+\hat{\xi}^{0}\|_{\hat{x}^{0}}^{\ast}<\frac{1}{2}$, then we have

[TABLE]

By using Lemma B.1(e), if $\|\hat{x}^{0}-x^{\ast}_{t_{0}}\|_{\hat{x}^{0}}<1$, then we have $\|\hat{x}^{0}-x^{\ast}_{t_{0}}\|_{x^{\ast}_{t_{0}}}\leq\frac{\|\hat{x}^{0}-x^{\ast}_{t_{0}}\|_{\hat{x}^{0}}}{1-\|\hat{x}^{0}-x^{\ast}_{t_{0}}\|_{\hat{x}^{0}}}$. To guarantee that $\|\hat{x}^{0}-x^{\ast}_{t_{0}}\|_{x^{\ast}_{t_{0}}}\leq\beta$, we require $\|\hat{x}^{0}-x^{\ast}_{t_{0}}\|_{\hat{x}^{0}}\leq\frac{\beta}{1+\beta}$. Combing this estimate and (78), we can show that $\|\hat{x}^{0}-x^{\ast}_{t_{0}}\|_{x^{\ast}_{t_{0}}}\leq\beta$ if $1-t_{0}\geq\frac{\beta}{(1+2\beta)\|\nabla{f}(\hat{x}^{0})+\hat{\xi}^{0}\|_{\hat{x}^{0}}^{\ast}}$. Therefore, we can choose

[TABLE]

to guarantee that $\|\hat{x}^{0}-x^{\ast}_{t_{0}}\|_{x^{\ast}_{t_{0}}}\leq\beta$. $\square$

B.6 Proof of Lemma 6.1: Approximate solution of the dual problem

Define the quadratic function

[TABLE]

It is easy to show that

[TABLE]

By using (37) and (38), we have $y^{k+1}\in\partial{\psi}(\tau_{k+1}z^{k+1})$, which implies $z^{k+1}\in\frac{1}{\tau_{k+1}}\partial{\psi^{\ast}}(y^{k+1})$. Hence, by the convexity of $\psi^{\ast}$, we have \frac{1}{\tau_{k+1}}\big{(}\psi^{\ast}(y^{k+1})-\psi^{\ast}(\bar{y}^{k+1})\big{)}\leq\langle z^{k+1},y^{k+1}-\bar{y}^{k+1}\rangle. Using this inequality and (80), if we define $\mathcal{P}_{k}(y):=q_{k}(y)+\frac{1}{\tau_{k+1}}\psi^{\ast}(y)$, then we have

[TABLE]

Next, using (38), we have $z^{k+1}=-\nabla{\varphi_{\tau_{k+1}}}(y^{k})-\nabla^{2}{\varphi}(y^{k})(y^{k+1}-y^{k}-\tilde{e}_{k})$ for $\tilde{e}_{k}\in\mathbb{R}^{p}$. Substituting this expression into the last inequality, we obtain

[TABLE]

The last inequality implies that if $\left\|\tilde{e}_{k}\right\|_{y^{k}}\leq\delta$, then $\mathcal{P}_{k}(y^{k+1})-\mathcal{P}_{k}(\bar{y}^{k+1})\leq\frac{\delta^{2}}{2}$. This completes the proof. $\square$

B.7 Proof of Theorem 6.2: An approximate solution for $x^{\star}$ of (1)

The statement (a) is trivial. We now prove (b) and (c) as follows.

(b) From (32), (39), and the definition of $\varphi$ in (30), we can show that

[TABLE]

Let $G_{k}:=\int_{0}^{1}\nabla^{2}{f^{\ast}}(-D^{\top}{y^{\star}}-sD^{\top}(y^{k}-y^{\star})ds$. By the mean-value theorem and [41, Corollary 4.1.4], we can derive from the above expression that

[TABLE]

Here, we note that $\nabla^{2}{\varphi}(y^{\star})=D\nabla^{2}{f^{\ast}}(-D^{\top}y^{\star})D^{\top}$. This inequality leads to (40) provided that $\|y^{k}-y^{\star}\|_{y^{\star}}<1$.

Since we apply (33) to solve the dual problem (29), by Theorem 4.3 or Theorem 4.6, the sequence $\left\{y^{k}\right\}$ satisfies $\|y^{k}-y^{\star}\|_{y^{\star}}\leq\hat{C}\sqrt{\sigma}^{k}$ for given constants $\hat{C}>0$ and $\sigma\in(0,1)$. Combining this relation and (40), we can show that $\left\{x^{k}\right\}$ converges linearly to $x^{\star}$ an optimal solution of (1).

(c) First, since \nabla{\varphi}_{\tau}(y)=D\nabla{f^{\ast}}(-D^{\top}y)-\big{(}\frac{1}{\tau}-1\big{)}\xi^{0}, using (32) we can write

[TABLE]

Next, from (38) we can write

[TABLE]

Combining these expressions, we can further estimate

[TABLE]

provided that $\|y^{\star}-y^{k}\|_{y^{k}}<1$. Here, the last inequality follows from the self-concordance of $\varphi_{\tau}$ with similar proof as [62, Theorem 1]. This proves (41).

Since we apply (33) to solve the dual problem (29), by Theorem 4.3 or Theorem 4.6, we have $\|y^{k}-y^{\star}\|_{y^{k}}\leq\hat{C}\sqrt{\sigma}^{k}$ and $\|y^{k+1}-y^{\star}\|_{y^{k}}\leq\hat{C}\sqrt{\sigma}^{k+1}$ for some constants $\hat{C}>0$ and $\sigma\in(0,1)$. In addition, $0\leq\frac{1}{\tau_{k+1}}-1=\frac{1-\tau_{k+1}}{\tau_{k+1}}\leq\frac{1-\tau_{k+1}}{\tau_{0}}\leq C_{1}\sqrt{\sigma}^{k+1}$ for some constant $C_{1}>0$. Using these bounds and $\|\tilde{e}_{k}\|_{y^{k}}\leq\delta_{k}$ in (41), we can show that

[TABLE]

Under the conditions of Theorem 4.3 or Theorem 4.6, we have $\|\xi^{0}\|_{y^{k}}^{\ast}=\langle\nabla^{2}{\varphi}(y^{k})^{-1}\xi^{0},\xi^{0}\rangle^{1/2}\leq M$ for some $M>0$ and $\delta_{k}\leq\frac{\beta}{113}\sigma^{k}$ for $\beta=0.05$, then we can easily show that

[TABLE]

provided that $\hat{C}<1$. If $D$ is invertible, then we can show that $x^{k+1}:=D^{-1}z^{k+1}$ is an approximate solution to $x^{\star}$ of (1). Moreover, $\|x^{k}-x^{\star}\|_{y^{k}}$ converges to zero at a linear rate. $\square$