A Universal Primal-Dual Convex Optimization Framework

Introduction

This paper constructs an algorithmic framework for the following convex optimization template:

When the set $\mathbf{A}\mathbf{x}\!-\!\mathbf{b}\!\in\!\mathcal{K}$ is absent in (1), other methods can be preferable to primal-dual algorithms. For instance, if $f$ has Lipschitz gradient, then we can use the accelerated proximal gradient methods by applying the proximal operator for the indicator function of the set $\mathcal{X}$ . However, as the problem dimensions become increasingly larger, the proximal tractability assumption can be restrictive. This fact increased the popularity of the generalized conditional gradient (GCG) methods (or Frank-Wolfe-type algorithms), which instead leverage the following Fenchel-type oracles

To this end, we propose a new primal-dual algorithmic framework that can exploit the sharp-operator of $f$ in lieu of its proximal operator. Our aim is to combine the flexibility of proximal primal-dual methods in addressing the general template (1) while leveraging the computational advantages of the GCG-type methods. As a result, we trade off the computational difficulty per iteration with the overall rate of convergence. While we obtain optimal rates based on the sharp-operator oracles, we note that the rates reduce to $\mathcal{O}\left(1/\epsilon^{2}\right)$ with the sharp operator vs. $\mathcal{O}\left(1/\epsilon\right)$ with the proximal operator when $f$ is completely non-smooth (cf. Definition 1.1). Intriguingly, the convergence rates are the same when $f$ is strongly convex. Unlike GCG-type methods, our approach can now handle nonsmooth objectives in addition to complex constraint structures as in (1).

Our algorithmic framework features a gradient method and its accelerated variant that operates on the dual formulation of (1). For the accelerated variant, we study an alternative to the universal accelerated method of based on FISTA since it requires less proximal operators in the dual. While the FISTA scheme is classical, our analysis of it with the Hölder continuous assumption is new. Given the dual iterates, we then use a new averaging scheme to construct the primal-iterates for the constrained template (1). In contrast to the non-adaptive weighting schemes of GCG-type algorithms, our weights explicitly depend on the local estimates of the Hölder constants $M_{\nu}$ at each iteration. Finally, we derive the worst-case complexity results. Our results are optimal since they match the computational lowerbounds in the sense of first-order black-box methods .

Section 2 briefly recalls primal-dual formulation of problem (1) with some standard assumptions. Section 3 defines the universal gradient mapping and its properties. Section 4 presents the primal-dual universal gradient methods (both the standard and accelerated variants), and analyzes their convergence. Section 5 provides numerical illustrations, followed by our conclusions. The supplementary material includes the technical proofs and additional implementation details.

Given an accuracy level $\epsilon>0$, a point $\mathbf{x}_{\epsilon}\in\mathcal{X}$ is said to be an $\epsilon$-solution of (1) if

Primal-dual preliminaries

In this section, we briefly summarise the primal-dual formulation with some standard assumptions. For the ease of presentation, we reformulate (1) by introducing a slack variable $\mathbf{r}$ as follows:

Let $\mathbf{z}\!:=\![\mathbf{x},\mathbf{r}]$ and $\mathcal{Z}\!:=\!\mathcal{X}\!\times\!\mathcal{K}$. Then, we have $\mathcal{D}\!:=\!\left\{\mathbf{z}\in\mathcal{Z}:\mathbf{A}\mathbf{x}\!-\!\mathbf{r}\!=\!\mathbf{b}\right\}$ as the feasible set of (3).

The Lagrange function associated with the linear constraint $\mathbf{A}\mathbf{x}-\mathbf{r}=\mathbf{b}$ is defined as $\mathcal{L}(\mathbf{x},\mathbf{r},\boldsymbol{\lambda}):=f(\mathbf{x})+\langle\boldsymbol{\lambda},\mathbf{A}\mathbf{x}-\mathbf{r}-\mathbf{b}\rangle$, and the dual function $d$ of (3) can be defined and decomposed as follows:

To characterize the primal-dual relation between (1) and (4), we require the following assumptions :

Universal gradient mappings

This section defines the universal gradient mapping and its properties.

We first adopt the composite convex minimization formulation of (4) in convex optimization for better interpretability as

where $G^{\star}=-d^{\star}$, and the correspondence between $(g,h)$ and $(d_{x},d_{r})$ is as follows:

Since $g$ and $h$ are generally non-smooth, FISTA and its proximal-based analysis are not directly applicable. Recall the sharp operator defined in (2), then $g$ can be expressed as

and we define the optimal solution to the $g$ subproblem above as follows:

The second term, $h$, depends on the structure of $\mathcal{K}$. We consider three special cases:

2 Hölder continuity of the dual universal gradient

Let $\nabla{g}(\cdot)$ be a subgradient of $g$, which can be computed as $\nabla{g}(\boldsymbol{\lambda})=\mathbf{b}-\mathbf{A}\mathbf{x}^{*}(\boldsymbol{\lambda})$. Next, we define

where $\nu\geq 0$ is the Hölder smoothness order. Note that the parameter $M_{\nu}$ explicitly depends on $\nu$ . We are interested in the case $\nu\in$, and especially the two extremal cases, where we either have the Lipschitz gradient that corresponds to $\nu=1$, or the bounded subgradient that corresponds to $\nu=0$.

We require the following condition in the sequel:

$\hat{M}(g):=\displaystyle\inf_{0\leq\nu\leq 1}M_{\nu}(g)<+\infty$.

Assumption A.2 is reasonable. We explain this claim with the following two examples. First, if $g$ is subdifferentiable and $\mathcal{X}$ is bounded, then $\nabla{g}(\cdot)$ is also bounded. Indeed, we have

Hence, we can choose $\nu=0$ and $\hat{M}_{\nu}(g)=2D_{\mathcal{X}}^{\mathbf{A}}<\infty$.

3 The proximal-gradient step for the dual problem

as an approximate quadratic surrogate of $g$. Then, we consider the following update rule:

For a given accuracy $\epsilon>0$, we define

Universal primal-dual gradient methods

We apply the universal gradient mappings to the dual problem (5), and propose an averaging scheme to construct $\{\bar{\mathbf{x}}_{k}\}$ for approximating $\mathbf{x}^{\star}$. Then, we develop an accelerated variant based on the FISTA scheme , and construct another primal sequence $\{\bar{\bar{\mathbf{x}}}_{k}\}$ for approximating $\mathbf{x}^{\star}$.

Our algorithm is shown in Algorithm 1. The dual steps are simply the universal gradient method in , while the new primal step allows to approximate the solution of (1).

Complexity-per-iteration: First, computing $\mathbf{x}^{*}(\boldsymbol{\lambda}_{k})$ at Step 1 requires the solution $\mathbf{x}^{*}(\boldsymbol{\lambda}_{k})\in[-\mathbf{A}^{T}\boldsymbol{\lambda}_{k}]^{\sharp}_{\mathcal{X},f}$. For many $\mathcal{X}$ and $f$, we can compute $\mathbf{x}^{*}(\boldsymbol{\lambda}_{k})$ efficiently and often in a closed form. Second, in the line-search procedure, we require the solution $\boldsymbol{\lambda}_{k,i}$ at Step 3.a, and the evaluation of $g(\boldsymbol{\lambda}_{k,i})$. The total computational cost depends on the proximal operator of $h$ and the evaluations of $g$. We prove below that our algorithm requires two oracle queries of $g$ on average.

The primal sequence $\left\{\bar{\mathbf{x}}_{k}\right\}$ generated by the Algorithm 1 satisfies

where $\widebar{M}_{\epsilon}$ is defined by (10), $\boldsymbol{\lambda}^{\star}\in\boldsymbol{\Lambda}^{\star}$ is an arbitrary dual solution, and $\epsilon$ is the desired accuracy.

The worst-case analytical complexity: We establish the total number of iterations $k_{\max}$ to achieve an $\epsilon$-solution $\bar{\mathbf{x}}_{k}$ of (1). The supplementary material proves that

where $\|\boldsymbol{\lambda}^{\star}\|_{}=\max{\{\|\boldsymbol{\lambda}^{\star}\|,1\}}$. This complexity is optimal for $\nu=0$, but not for $\nu>0$ .

At each iteration $k$, the linesearch procedure at Step 3 requires the evaluations of $g$. The supplementary material bounds the total number $N_{1}(k)$ of oracle queries, including the function $G$ and its gradient evaluations, up to the $k$th iteration as follows:

Hence, we have $N_{1}(k)\approx 2(k+1)$, i.e., we require approximately two oracle queries at each iteration on the average.

2 Accelerated universal primal-dual gradient method

Complexity per-iteration: The per-iteration complexity of Algorithm 2 remains essentially the same as that of Algorithm 1.

The primal sequence $\left\{\bar{\bar{\mathbf{x}}}_{k}\right\}$ generated by the Algorithm 2 satisfies

where $\widebar{M}_{\epsilon}$ is defined by (10), $\boldsymbol{\lambda}^{\star}\in\boldsymbol{\Lambda}^{\star}$ is an arbitrary dual solution, and $\epsilon$ is the desired accuracy.

The worst-case analytical complexity: The supplementary material proves the following worst-case complexity of Algorithm 2 to achieve an $\epsilon$-solution $\bar{\bar{\mathbf{x}}}_{k}$:

This worst-case complexity is optimal in the sense of first-order black box models .

The line-search procedure at Step 3 of Algorithm 2 also terminates after a finite number of iterations. Similar to Algorithm 1, Algorithm 2 requires $1$ gradient query and $i_{k}$ function evaluations of $g$ at each iteration. The supplementary material proves that the number of oracle queries in Algorithm 2 is upperbounded as follows:

Roughly speaking, Algorithm 2 requires approximately two oracle query per iteration on average.

Numerical experiments

This section illustrates the scalability and the flexibility of our primal-dual framework using some applications in the quantum tomography (QT) and the matrix completion (MC).

We consider the QT problem which aims to extract information from a physical quantum system. A $q$-qubit quantum system is mathematically characterized by its density matrix, which is a complex $p\times p$ positive semidefinite Hermitian matrix $\mathbf{X}^{\natural}\in\mathcal{S}^{p}_{+}$, where $p=2^{q}$. Surprisingly, we can provably deduce the state from performing compressive linear measurements $\mathbf{b}=\mathcal{A}(\mathbf{X})\in\mathcal{C}^{n}$ based on Pauli operators $\mathcal{A}$ . While the size of the density matrix grows exponentially in $q$, a significantly fewer compressive measurements (i.e., $n\!=\!\mathcal{O}(p\log p)$) suffices to recover a pure state $q$-qubit density matrix as a result of the following convex optimization problem:​​

where the constraint ensures that $\mathbf{X}^{\star}$ is a density matrix. The recovery is also robust to noise .

Since the objective function has Lipschitz gradient and the constraint (i.e., the Spectrahedron) is tuning-free, the QT problem provides an ideal scalability test for both our framework and GCG-type algorithms. To verify the performance of the algorithms with respect to the optimal solution in large-scale, we remain within the noiseless setting. However, the timing and the convergence behavior of the algorithms remain qualitatively the same under polarization and additive Gaussian noise.

To this end, we generate a random pure quantum state (e.g., rank-1 $\mathbf{X}^{\natural}$), and we take $n=2p\log p$ random Pauli measurements. For $q=14$ qubits system, this corresponds to a $268^{\prime}435^{\prime}456$ dimensional problem with $n=138^{\prime}099$ measurements. We recast (19) into (1) by introducing the slack variable $\mathbf{r}=\mathcal{A}(\mathbf{X})-\mathbf{b}$.

We compare our algorithms vs. the Frank-Wolfe method, which has optimal convergence rate guarantees for this problem, and its line-search variant. Computing the sharp-operator $[\mathbf{x}]^{\sharp}$ requires a top-eigenvector $\mathbf{e}_{1}$ of $\mathcal{A}^{*}(\boldsymbol{\lambda})$, while evaluating $g$ corresponds to just computing the top-eigenvalue $\sigma_{1}$ of $\mathcal{A}^{*}(\boldsymbol{\lambda})$ via a power method. All methods use the same power method subroutine, which is implemented in MATLAB’s eigs function. We set $\epsilon=2\times 10^{-4}$ for our methods and have a wall-time $2\times 10^{4}$s in order to stop the algorithms. However, our algorithms seems insensitive to the choice of $\epsilon$ for the QT problem.

Figure 1 illustrates the iteration and the timing complexities of the algorithms. UniPDGrad algorithm, with an average of $1.978$ line-search steps per iteration, has similar iteration and timing performance as compared to the standard Frank-Wolfe scheme with step-size $\gamma_{k}=2/(k+2)$. The line-search variant of Frank-Wolfe improves over the standard one; however, our accelerated variant, with an average of $1.057$ line-search steps, is the clear winner in terms of both iterations and time. We can empirically improve the performance of our algorithms even further by adapting a similar line-search strategy in the weighting step as Frank-Wolfe, i.e., by choosing the weights $w_{k}$ in a greedy fashion to minimize the objective function. The practical improvements due to line-search appear quite significant.

2 Matrix completion with MovieLens dataset

Convex formulations involving the nuclear norm have been shown to be quite effective in estimating low-rank matrices from limited number of measurements . For instance, we can solve

with Frank-Wolfe-type methods, where $\kappa$ is a tuning parameter, which may not be available a priori. We can also solve the following parameter-free version

While the nonsmooth objective of (21) prevents the tuning parameter, it clearly burdens the computational efficiency of the convex optimization algorithms.

We apply our algorithms to (20) and (21) using the MovieLens 100K dataset. Frank-Wolfe algorithms cannot handle (21) and only solve (20). For this experiment, we did not pre-process the data and took the default ub test and training data partition. We start out algorithms form $\boldsymbol{\lambda}_{0}=\mathbf{0}^{n}$, we set the target accuracy $\epsilon=10^{-3}$, and we choose the tuning parameter $\kappa=9975/2$ as in . We use lansvd function (MATLAB version) from PROPACK to compute the top singular vectors, and a simple implementation of the power method to find the top singular value in the line-search, both with $10^{-5}$ relative error tolerance.

The first two plots in Figure 2 show the performance of the algorithms for (20). Our metrics are the normalized objective residual and the root mean squared error (RMSE) calculated for the test data. Since we do not have access to the optimal solutions, we approximated the optimal values, $\varphi^{\star}$ and RMSE⋆, by $5000$ iterations of AccUniPDGrad. Other two plots in Figure 2 compare the performance of the formulations (20) and (21) which are represented by the empty and the filled markers, respectively. Note that, the dashed line for AccUniPDGrad corresponds to the line-search variant, where the weights $w_{k}$ are chosen to minimize the feasibility gap. Additional details about the numerical experiments can be found in the supplementary material.

Conclusions

This paper proposes a new primal-dual algorithmic framework that combines the flexibility of proximal primal-dual methods in addressing the general template (1) while leveraging the computational advantages of the GCG-type methods. The algorithmic instances of our framework are universal since they can automatically adapt to the unknown Hölder continuity properties implied by the template. Our analysis technique unifies Nesterov’s universal gradient methods and GCG-type methods to address the more broadly applicable primal-dual setting. The hallmarks of our approach includes the optimal worst-case complexity and its flexibility to handle nonsmooth objectives and complex constraints, compared to existing primal-dual algorithm as well as GCG-type algorithms, while essentially preserving their low cost iteration complexity.

This work was supported in part by ERC Future Proof, SNF 200021-146750 and SNF CRSII2-147633. We would like to thank Dr. Stephen Becker of University of Colorado at Boulder for his support in preparing the numerical experiments.

References

References

Appendix A The key estimate of the proximal-gradient step

Lemma 2 in , which we present below as Lemma A.1, provides key properties for constructing universal gradient algorithms. We refer to for the proof of this lemma.

This lemma provides an approximate quadratic upper bound for $g$. However, it depends on the choice of the inexactness parameter $\delta$ and the smoothness parameter $\nu$. If $\nu=1$, then $M$ can be set to the Lipschitz constant $M_{1}$, and it becomes independent of $\delta$.

The algorithms that we develop in this paper are based on the proximal-gradient step (9) on the dual objective function $G$. This update rule guarantees the following estimate:

Let $Q_{M}$ be the quadratic model of $g$. If $\boldsymbol{\lambda}_{k+1}$, which is defined by (9), satisfies

We note that the optimality condition of (9) is

which can be written as $\hat{\boldsymbol{\lambda}}_{k}-\boldsymbol{\lambda}_{k+1}\in M_{k}^{-1}(\nabla{g}_{k}(\hat{\boldsymbol{\lambda}}_{k})+\partial{h}(\boldsymbol{\lambda}_{k+1}))$. Let $\nabla{h}(\hat{\boldsymbol{\lambda}}_{k+1})\in\partial{h}(\boldsymbol{\lambda}_{k+1})$ be a subgradient of $h$ at $\boldsymbol{\lambda}_{k+1}$. Then, we have

where the last inequality directly follows the convexity of $h$. ∎

Clearly, (22) holds if $M_{k}\geq\widebar{M}_{\epsilon}$, which is defined by (10), due to Lemma A.1, whenever $\delta_{k}\!=\!\epsilon>0$.

If $\nu$ and $M_{\nu}$ are known, we can set $M_{k}=\widebar{M}_{\epsilon}$, then the condition (22) is automatically satisfied. However, we do not know $\nu$ and $M_{\nu}$ a priori in general. In this case, $M_{k}$ can be determined via a line-search procedure on the condition (22).

The following lemma guarantees that the line-search procedure in Algorithms 1 and 2 terminates after a finite number of line-search iterations.

The line-search procedure in Algorithm 1 terminates after at most

Similarly, the line-search procedure in Algorithm 2 terminates after at most

Now, we show that the line-search procedure in Algorithm 2 is also finite. By the updating rule of $t_{k}$, we have $t_{k+1}:=0.5(1+\sqrt{1+4t_{k}^{2}})\leq 0.5(1+(1+2t_{k}))=t_{k}+1$. By induction and $t_{0}=1$, we have $t_{k}\leq k+1$. Using the definition (10) of $\widebar{M}_{\delta_{k}}$ with $\delta_{k}=\frac{\epsilon}{t_{k}}$ and $t_{k}\leq k+1$, we can show that

Next, we note that the condition (22) holds whenever $M_{k,i}\geq\widebar{M}_{\delta_{k}}$. However, since $M_{k,i}=2^{i}M_{k,0}\geq 2^{i}{M}_{-1}$, by using (24), it is sufficient to show that the following condition holds for a finite $i$:

This condition leads to $i\geq\log_{2}\left(\left[\frac{k+1}{\epsilon}\right]^{\frac{1-\nu}{1+\nu}}M_{\nu}^{\frac{2}{1+\nu}}\right)-\log_{2}({M}_{-1})$. Hence, at the $k$th iteration, we require at most $i_{k}=\left\lfloor\log_{2}\left(\frac{k+1}{\epsilon}\right)+\log_{2}\left(\frac{M_{\nu}^{\frac{2}{1+\nu}}}{{M}_{-1}}\right)\right\rfloor+1$ line-search iterations, which is finite. ∎

Appendix B Convergence analysis of the universal primal-dual gradient algorithm

In this section, we analyze the convergence of the Algorithm 1 (UniPDGrad). We first provide the convergence guarantee of the dual function in Theorem B.1. Then, we prove the convergence rate and the worst-case complexity given in Theorem 4.1.

Let $\left\{\boldsymbol{\lambda}_{k}\right\}$ be the sequence generated by UniPDGrad. Then,

For $\widebar{M}_{\epsilon}$ defined by (10), since the line-search is successful as shown in Lemma A.1, the condition (22) is satisfied at iteration $i$ with $M_{i}\leq 2\widebar{M}_{\epsilon}$. The following inequality directly follows Lemma A.2 considering the convexity of $g$:

Taking the weighted sum of this inequality over $i$, we get

B.2 The proof of Theorem 4.1: Convergence rate of the primal sequence

We use the following three expressions to relate the convergence in the dual sequence to the convergence in the primal sequence:

Taking the weighted sum of this inequality over $i$ and considering the convexity of $f$, we get

Setting $\boldsymbol{\lambda}=\mathbf{0}^{n}$, we get the bound on the right hand side of (15),

The inequality on the left hand side of (11) follows the following saddle point formulation:

$\forall\mathbf{r}\in\mathcal{K}$ and $\forall\mathbf{x}\in\mathcal{X}$, where the last inequality holds due to Cauchy-Schwarz inequality. The proof of the convergence rate in the objective residual (11) follows by setting $\mathbf{x}=\bar{\mathbf{x}}$ in (32).

Next, we prove the convergence rate of the feasibility gap (12). We start from the following saddle point formulation:

Substituting this estimate with $\mathbf{x}=\bar{\mathbf{x}}_{k}$ into (31), we get the following inequality:

Using Cauchy-Schwarz inequality, this implies

Solving this inequality for $\|\mathbf{A}\bar{\mathbf{x}}_{k}-\mathbf{b}-\bar{\mathbf{r}}\|$, we get

We note that $S_{k}\geq\frac{k+1}{2\widebar{M}_{\epsilon}}$, and this completes the proof. ∎

Hence, the worst-case complexity to obtain an $\epsilon$-solution of (1) in the sense of Definition 1.1 is

Next, we estimate the total number of oracle quires in UniPDGrad, as in . The total number of oracle quires up to the iteration $k$ is given by $N_{1}(k)=\sum_{j=0}^{k}(i_{j}+1)$. However, since $i_{j}-1=\log_{2}(M_{j}/M_{j-1})$, we have

It remains to use $M_{k}\leq 2\widebar{M}_{\epsilon}$ to obtain (14).

Appendix C Convergence analysis of the accelerated universal primal-dual algorithm

We now analyze the convergence of AccUniPDGrad (Algorithm 2) in terms of the objective residual and the feasibility gap.

The parameter $M_{k}$ is determined based on the following line-search condition:

Next, we simplify the scheme (33) in the following lemma:

The scheme (33) can be restated as follows:

where $\hat{\boldsymbol{\lambda}}_{0}=\boldsymbol{\lambda}_{0}$ and $t_{0}=1$, and $M_{k}$ is determined based on the line-search condition (35).

This dual scheme is of the FISTA form , except for the line-search step.

Hence $\hat{\boldsymbol{\lambda}}_{k+1}=\boldsymbol{\lambda}_{k+1}+\frac{t_{k}-1}{t_{k+1}}(\boldsymbol{\lambda}_{k+1}-\boldsymbol{\lambda}_{k})$, which is the third step of (36).

Next, from the condition (34), we have $t_{k+1}^{2}-t_{k+1}-t_{k}^{2}=0$. Hence, $t_{k+1}=\frac{1}{2}\left[1+\sqrt{1+4t_{k}^{2}}\right]$, which is exactly the second step of (36). ∎

and we set $\boldsymbol{\lambda}=\boldsymbol{\lambda}_{k}$ in (38), and then subtract $G^{\star}$ from the both sides, that results in the following inequality:

We obtain the following estimate by summing the two inequalities that we get by multiplying (39) by $\tau_{k}$ and (40) by $(1-\tau_{k})$, and then dividing the resulting estimate by $M_{k}\tau_{k}^{2}$:

Next, we sum this inequality over $k$ as follows:

where the second inequality follows $\tau_{0}=1$ and $\frac{(1-\tau_{k})}{M_{k}\tau_{k}^{2}}\leq\frac{1}{M_{k-1}\tau_{k-1}^{2}}$ for $k=1,2,\dots$, which holds since $M_{k}\geq M_{k-1}$. This implies the followings:

Now, we use the following expressions to map this estimate into the primal sequence:

Then, considering the convexity of $f$, we get

where we obtain the second inequality by setting $\boldsymbol{\lambda}=\mathbf{0}^{n}$.

We can reformulate (34) as $\frac{1}{\tau_{k}}=\frac{1}{\tau_{k}^{2}}-\frac{1}{\tau_{k-1}^{2}}$. Using this relation, $M_{0}\leq M_{i}\leq M_{k}\leq 2\widebar{M}_{\epsilon\tau_{k}}=2t_{k}^{\frac{1-\nu}{1+\nu}}\widebar{M}_{\epsilon}\leq 2(k+2)^{\frac{1-\nu}{1+\nu}}\widebar{M}_{\epsilon}$ and $\frac{k+2}{2}\leq t_{k}<k+2$ for $i=0,1,\dots,k$, we can show that

We get the bound on the right hand side of (15) by substituting (43) into (42). The inequality on the left hand side of (15) follows the saddle point formulation (32) by setting $\mathbf{x}=\bar{\bar{\mathbf{x}}}_{k}$.

Finally, we prove the convergence rate in the feasibility gap (16). By the same arguments as in the proof of Theorem 4.1, we have

We complete the proof by substituting (43) into this estimate. ∎

C.2 The worst-case complexity analysis

We analyze the worst-case complexity of AccUniPDGrad algorithm to achieve an $\epsilon$-solution $\bar{\bar{\mathbf{x}}}_{k}$. For simplicity, we consider the case $\boldsymbol{\lambda}_{0}=\mathbf{0}^{n}$ without loss of generality. Then, we require

due to the Theorem 4.2, where $\|\boldsymbol{\lambda}^{\star}\|_{}:=\max{\{\|\boldsymbol{\lambda}^{\star}\|,1\}}$. By solving this inequality, we get

Using the definition (10) of $\widebar{M}_{\epsilon}$ and considering the fact that $\left[\frac{1-\nu}{1+\nu}\right]^{\frac{1-\nu}{1+\nu}}\leq 1$ for $\nu\in$, we find the maximum number of iterations that satisfies the above inequality as follows:

Hence, the worst-case complexity to obtain an $\epsilon$-solution of (1) in the sense of Definition 1.1 is

which is optimal in the sense of first-order black box models .

Next, we consider the number of oracle quires in AccUniPDGrad. At iteration $k$, the algorithm requires $i_{k}\!+\!2$ function evaluations of $g$, as we need $i_{k}\!+\!1$ in the line-search and one for $g(\hat{\boldsymbol{\lambda}}_{k})$. Hence, the total number of oracle quires up to the iteration $k$ is $N_{2}(k)=\sum_{j=0}^{k}(i_{j}+2)$. Since $i_{j}=\log_{2}(M_{j}/M_{j-1})$, we have

Using the same argument as in the proof of Lemma A.3, we have $M_{k}\leq 2\widebar{M}_{\epsilon\tau_{k}}\leq 2\left[\frac{k+1}{\epsilon}\right]^{\frac{1-\nu}{1+\nu}}M_{\nu}^{\frac{2}{1+\nu}}$. Hence, we obtain (18) as

Appendix D The implementation details

In this section, we specify key steps of UniPDGrad and AccUniPDGrad for two important applications that we used in Section 5. We also provide an analytic step-size that guarantees the line-search condition without function evaluation.

We performed the experiments in MATLAB, using a computational resource with 4 CPUs of 2.40 GHz and 16 GB memory space for the matrix completion, and 16 CPUs of 2.40 GHz and 512 GB memory space for the quantum tomography problem.

In both quantum tomography and the matrix completion problems, we consider some problem formulations from the following convex optimization template that involves a quadratic cost:

Evaluation of the sharp-operator corresponding to the objective function $1/2\|\mathcal{A}(\mathbf{x})-\mathbf{b}\|^{2}$ requires a significant computational effort. Yet, by introducing the slack variable $\mathbf{r}=\mathcal{A}(\mathbf{x})-\mathbf{b}$, we can write an equivalent problem as

We can write the Lagrange function associated with the linear constraint as

from which we can derive the (negation of the) dual function

where $\mathbf{x}^{\ast}(\boldsymbol{\lambda})\in[\mathcal{A}^{T}(\boldsymbol{\lambda})]^{\sharp}_{\mathcal{X}}\equiv\arg\max_{\mathbf{x}\in\mathcal{X}}~{}\langle\mathcal{A}^{T}(\boldsymbol{\lambda}),\mathbf{x}\rangle$.

For the special case, $\mathcal{X}$ is a norm ball, i.e., $\mathcal{X}\equiv\left\{\mathbf{x}:\|\mathbf{x}\|\leq\kappa\right\}$, we can simplify (44) as follows:

Computing an analytical step-size: Now, we consider the line-search procedure in UniPDGrad and AccUniPDGrad. Since $h(\boldsymbol{\lambda})$ term is absent in these problems, the line-search condition (22) can be simplified as

where we use the notational convention $\hat{\boldsymbol{\lambda}}_{k}=\boldsymbol{\lambda}_{k}$ and $\delta_{k}=\epsilon$ for UniPDGrad, and $\delta_{k}=\epsilon/t_{k}$ for AccUniPDGrad. Using the definition (45), we can upper bound $g(\hat{\boldsymbol{\lambda}}_{k}-\alpha_{k}\nabla{g}(\hat{\boldsymbol{\lambda}}_{k}))$ by

The condition (46) holds if $U(\alpha_{k})=g(\hat{\boldsymbol{\lambda}}_{k})-\frac{\alpha_{k}}{2}\|\nabla{g}(\hat{\boldsymbol{\lambda}}_{k})\|^{2}+\delta_{k}/2$. Solving this second order equation, we obtain $\alpha_{k}$ explicitly as

where $P:=\|\nabla{g}(\hat{\boldsymbol{\lambda}}_{k})\|^{2}+2\kappa\|\mathcal{A}^{T}(\nabla{g}(\hat{\boldsymbol{\lambda}}_{k}))\|-2\langle\hat{\boldsymbol{\lambda}}_{k}-\mathbf{b},\nabla{g}(\hat{\boldsymbol{\lambda}}_{k})\rangle$. Note that, we can use this method to find a good estimate for the initial smoothness constant $M_{-1}$ in the initialization step.

D.2 Constrained convex optimization involving a norm cost

Now, we consider the second application, which is reformulated as

Clearly, the dual components $g$ and $h$ defined in (6) can be expressed as:

where $\|\cdot\|$ represents the Euclidean norm for vectors and the spectral norm for matrices. In (21), we consider a special case where $\mathcal{K}\equiv\left\{\mathbf{0}^{n}\right\}$, hence $h(\boldsymbol{\lambda})=0$.

Clearly, $\mathbf{X}^{*}(\boldsymbol{\lambda})=\sigma_{1}\mathbf{e}_{1}\mathbf{e}_{1}^{T}\in[\mathcal{A}^{T}(\boldsymbol{\lambda})]^{\sharp}_{\psi}$, where $\sigma_{1}=\|\mathcal{A}^{T}(\boldsymbol{\lambda})\|$ is the top singular value of $\mathcal{A}^{T}(\boldsymbol{\lambda})$ and $\mathbf{e}_{1}$ is the associated left singular vector. Hence, we can write the (sub)gradient of g as

We can compute both $\sigma_{1}$ and $\mathbf{e}_{1}$ efficiently by using the power method or the Lanczos algorithm.

References