Analyzing Random Permutations for Cyclic Coordinate Descent

TL;DR

This paper investigates the performance of random permutation cyclic coordinate descent (RPCD) on convex quadratic problems, showing it can outperform standard cyclic methods and behave similarly to fully-random approaches, supported by theoretical analysis.

Contribution

It introduces a class of convex quadratic problems where RPCD outperforms CCD and provides a convergence analysis explaining this empirical advantage.

Findings

RPCD outperforms CCD on certain convex quadratic problems.

RPCD exhibits convergence similar to RCD in these cases.

Theoretical analysis supports empirical observations.

Abstract

We consider coordinate descent methods on convex quadratic problems, in which exact line searches are performed at each iteration. (This algorithm is identical to Gauss-Seidel on the equivalent symmetric positive definite linear system.) We describe a class of convex quadratic problems for which the random-permutations version of cyclic coordinate descent (RPCD) outperforms the standard cyclic coordinate descent (CCD) approach, yielding convergence behavior similar to the fully-random variant (RCD). A convergence analysis is developed to explain the empirical observations.

Figures

Tables

Table 1. Table 1 . Comparison of CCD, RPCD, and RCD on the matrix ( 1.13 ) with n = 100 𝑛 100 n=100 and ϵ = δ italic-ϵ 𝛿 \epsilon=\delta .

$\delta$	1.0000 (-03)	3.0000 (-03)	1.0000 (-02)	3.0000 (-02)	1.0000 (-01)
$1-{\rho }_{\text{CCD}}(\delta ,\text{observed})$	3.4122 (-04)	3.3170 (-04)	3.3527 (-04)	6.1266 (-04)	8.1036 (-04)
$1-\rho {\left(C\right)}^2$	5.9018 (-06)	1.7170 (-05)	6.0912 (-05)	1.9453 (-04)	7.5546 (-04)
$1-{\rho }_{\text{RCD}}(\delta ,\text{observed})$	2.6814 (-03)	5.8265 (-03)	2.1983 (-02)	6.8824 (-02)	1.4427 (-01)
$1-{\rho }_{\text{RCD}}(\delta ,\text{predicted})$	1.9940 (-03)	5.9466 (-03)	1.9419 (-02)	5.5047 (-02)	1.5364 (-01)
$1-{\rho }_{\text{RPCD}}(\delta ,\text{observed})$	2.7048 (-03)	6.3637 (-03)	2.1723 (-02)	6.9230 (-02)	2.0842 (-01)
Benchmark $2\delta$	2.0000 (-03)	6.0000 (-03)	2.0000 (-02)	6.0000 (-02)∗	2.0000 (-01)∗

Table 2. Table 2 . Comparison of CCD, RPCD, and RCD on the matrix ( 1.13 ) with n = 100 𝑛 100 n=100 and ϵ = δ / 10 italic-ϵ 𝛿 10 \epsilon=\sqrt{\delta/10} .

$\delta$	1.0000 (-03)	3.0000 (-03)	1.0000 (-02)	3.0000 (-02)	1.0000 (-01)
$1-{\rho }_{\text{CCD}}(\delta ,\text{observed})$	2.2372 (-04)	3.9800 (-04)	3.3538 (-04)	2.8511 (-04)	7.9319 (-04)
$1-\rho {\left(C\right)}^2$	2.7954 (-05)	5.0165 (-05)	9.7958 (-05)	2.3542 (-04)	7.8096 (-04)
$1-{\rho }_{\text{RCD}}(\delta ,\text{observed})$	2.6143 (-03)	8.6962 (-03)	1.7869 (-02)	5.8402 (-02)	1.4545 (-01)
$1-{\rho }_{\text{RCD}}(\delta ,\text{predicted})$	1.9763 (-03)	5.8634 (-03)	1.9019 (-02)	5.3824 (-02)	1.5364 (-01)
$1-{\rho }_{\text{RPCD}}(\delta ,\text{observed})$	2.8377 (-03)	7.1350 (-03)	2.1157 (-02)	6.6712 (-02)	2.0501 (-01)
Benchmark $2\delta$	2.0000 (-03)	6.0000 (-03)∗	2.0000 (-02)∗	6.0000 (-02)∗	2.0000 (-01)∗

Equations

$$\min\,f(x),\quad\mbox{where $f:\mathbb{R}^{n}\to\mathbb{R}$ is smooth and convex,}$$

$${\nabla }_{i}f(x)=[\nabla f(x){]}_i,e_i=(0,\dots ,0,1,0,\dots ,0{)}^T,$$

$$f(x)=\frac{1}{2}{x}^{T}Ax,$$

$$f(x)=\frac{1}{2}{x}^{T}Ax-b^{T}x$$

$$P^{T}AP=L_P+{\Delta }_P+L_P^T,$$

$$C_P:=-(L_P+{\Delta }_P{)}^{-1}{L}_P^T,$$

$$x^{ln}=(P_l{C}_{P_l}{P}_l^T)x^{(l-1)n},$$

$$x^{\ell n}=(P_{\ell }{C}_{P_{\ell }}{P}_{\ell }^T)(P_{\ell -1}{C}_{P_{\ell -1}}{P}_{\ell -1}^T)\dots (P_1{C}_{P_1}{P}_1^T)x^0,$$

$$f(x^{\ell n})=\frac{1}{2}(x^0{)}^T\left((P_1{C}_{P_1}^T{P}_1^T)\dots (P_{\ell }{C}_{P_{\ell }}^T{P}_{\ell }^T)A(P_{\ell }{C}_{P_{\ell }}{P}_{\ell }^T)\dots (P_1{C}_{P_1}{P}_1^T)\right)x^0.$$

$${\mathbb{E}}_{P_1,P_2,\dots ,P_{\ell }}f(x^{\ell n}).$$

$$A=\delta V\Sigma V^T+(1-\delta )11^T,$$

$$A:=\delta I+(1-\delta)\mathbf{1}\mathbf{1}^{T},\quad\mbox{where $\delta\in(0,n/(n-1))$.}$$

$$x^{ln}=(P_{l}C{P}_l^T)x^{(l-1)n},l=1,2,3,\dots .$$

$$B_u:=\delta I+(1-\delta )u{u}^T,$$

$$A_{ϵ}:=\delta I+(1-\delta )11^T+ϵD,$$

$$\delta \in (0,n/(n-1))$$

$$D=\mathrm{diag}(d)$$

$$0<\delta \le ϵ,\left|{\rho }_1\right|{ϵ}^2<\delta \ll 1,nϵ\le 1.$$

$$1-1.4\delta ,$$

$$\Vert {\mathbf{r}}_1\Vert \le 1,\Vert {\mathbf{R}}_1\Vert \le 1.$$

$$\mathbf{R_{1}}v\mathbf{1}^{T}=\rho_{1}\mathbf{r_{1}}\mathbf{1}^{T}\;\;\mbox{provided $\|v\|\leq\rho_{1}$}.$$

$$\underset{i=1,2,\dots ,n}{\min }|u_i|=\sqrt{\frac{\delta }{\delta +ϵ}},\underset{i=1,2,\dots ,n}{\max }|u_i|=1.$$

$$A_{ϵ}:=U^{-1}{B}_u{U}^{-1}=\delta U^{-2}+(1-\delta )11^T,$$

$$\frac{1}{2}({\tilde{x}}^k{)}^T{B}_u{\tilde{x}}^k=\frac{1}{2}(x^k{)}^T{A}_{ϵ}{x}^k,k=0,1,2,\dots .$$

$$\delta I+(1-\delta )11^T\le A_{ϵ}\le (1+ϵ)\left({\delta }^{\prime }I+(1-{\delta }^{\prime })11^T\right),$$

$$P^T{A}_{ϵ}P$$

$$=(1-\delta )E+(\delta I+ϵD_P)+(1-\delta )E^T,$$

$$D_P:=P^{T}DP,E:=\left[\begin{array}{cccccc}0& 0& 0& \dots & 0& 0\\ 1& 0& 0& \dots & 0& 0\\ 1& 1& 0& \dots & 0& 0\\ ⋮& ⋮& ⋮& & ⋮& ⋮\\ 1& 1& 1& \dots & 1& 0\end{array}\right].$$

$$C_P:=-(1-\delta ){\left[(1-\delta )E+(I+ϵD_P)\right]}^{-1}{E}^T.$$

$${\mathbb{E}}_{P_1,P_2,\dots ,P_{\ell }}f(x^{\ell n}),\ell =1,2,\dots ,$$

$${\bar{A}}_{ϵ}^{(0)}$$

$${\bar{A}}_{ϵ}^{(1)}$$

Full text

This work was supported by NSF Awards IIS-1447449, 1628384, 1634597, and 1740707; ONR Award N00014-13-1-0129; AFOSR Award FA9550-13-1-0138, and Subcontracts 3F-30222 and 8F-30039 from Argonne National Laboratory; and Award N660011824020 from the DARPA Lagrange Program. Ching-pei Lee was at the University of Wisconsin-Madison when this work was originally submitted.

1. Introduction

The coordinate descent (CD) approach for solving the problem

[TABLE]

follows the framework of Algorithm 1. We denote

[TABLE]

where the single nonzero in $e_{i}$ appears in position $i$. Epochs (indicated by the counter $\ell$) encompass cycles of inner iterations (indicated by $j$). At each iteration $k$, one component of $x$ is selected for updating; a steplength parameter $\alpha_{k}$ is applied to the negative gradient of $f$ with respect to this component.

The choice of coordinate $i=i(\ell,j)$ to be updated at inner iteration $j$ of epoch $\ell$ differs between variants of CD, as follows:

• •

For “cyclic CD” (CCD), we choose $i(\ell,j)=j+1$.

• •

For “fully randomized CD,” also known as “stochastic CD,” and abbreviated as RCD, we choose $i(\ell,j)$ uniformly at random from $\{1,2,\dotsc,n\}$ and independently at each iteration.

• •

For “random-permutations CD” (abbreviated as RPCD), we choose $\pi_{\ell+1}$ at the start of epoch $\ell$ to be a random permutation of the index set $\{1,2,\dotsc,n\}$ (chosen uniformly at random from the space of random permutations), then set $i(\ell,j)$ to be the $(j+1)$th entry in $\pi_{\ell+1}$, for $j=0,1,2,\dotsc,n-1$.

Note that $x^{ln}$ denotes the value of $x$ after $l$ epochs.

We consider in this paper problems in which $f$ is a strictly convex quadratic, that is

[TABLE]

with $A$ symmetric positive definite. Even this restricted class of functions reveals significant diversity in convergence behavior between the three variants of CD described above. The minimizer of (1.3) is obviously $x^{*}=0$. Although (1.3) does not contain a linear term, it is straightforward to extend our results to the case for problems of the form

[TABLE]

by replacing $x^{0}$ in several places of our analysis with $x^{0}-x^{*}$, where $x^{*}=A^{-1}b$ is the minimizer of this problem. We assume that the choice of $\alpha_{k}$ in Algorithm 1 is the exact minimizer of $f$ along the chosen coordinate direction. The resulting approach is thus equivalent to the Gauss-Seidel method applied to the linear system $Ax=0$. The variants CCD, RCD, RPCD can be interpreted as different cyclic / randomized variants of Gauss-Seidel for this system.

In the RPCD variant, we can express a single epoch as follows. Letting $P$ be the permutation matrix corresponding to the permutation $\pi$ on this epoch, we split the symmetrically permuted Hessian into strictly triangular and diagonal parts as follows:

[TABLE]

where $L_{P}$ is strictly lower triangular and $\Delta_{P}$ is diagonal. We then define

[TABLE]

so that the epoch indexed by $l-1$ can be written as follows:

[TABLE]

where $P_{l}$ denotes the matrix corresponding to permutation $\pi_{l}$. By recursing to the initial point $x^{0}$, we obtain after $\ell$ epochs that

[TABLE]

yielding a function value of

[TABLE]

We analyze convergence in terms of the expected value of $f$ after $\ell$ epochs for any given $x^{0}$, with the expectation taken over the permutations $P_{1},P_{2},\dotsc,P_{\ell}$, that is,

[TABLE]

1.1. Previous Work

Convergence of RCD is analyzed in [5], showing that when the objective is strongly convex, the method requires $O((nL_{\max}/\mu)|\log\hat{\epsilon}|)$ iterations to reach an objective function value that is within $\hat{\epsilon}$ of the optimal value, in expectation, when the coordinates are sampled in a uniform random manner, where $\mu$ is the modulus of strong convexity and $L_{\max}$ is the maximum coordinate-wise Lipschitz constant for the gradient. This rate can be improved to $O((nL_{\text{avg}}/\mu)|\log\hat{\epsilon}|)$ if the sampling probability for each coordinate is proportional to the coordinate-wise Lipschitz constants, where $L_{\text{avg}}$ is the average of these constants. On the other hand, the best known convergence rate of CCD for convex quadratic problems, given by [8], has an iteration complexity for reaching an $\hat{\epsilon}$-accurate solution deterministically that can be $O(n^{2})$ times slower than that for RCD in the worst case. The best worst-case convergence guarantees for RPCD so far are still identical to those for CCD. (The analyses for CCD assume only that each coordinate is processed exactly once per epoch, and are indifferent to the fact that the ordering of coordinates can change on each epoch, as in RPCD.) However, in practice it is sometimes observed that RPCD behaves in a manner more similar to RCD than CCD (see, for example, the experiments in [7] and the talks [11, 12]), and a rigorous explanation for the general convergence rate for the expected objective of RPCD over the random permutations has been difficult to obtain. Some trials have been conducted to tackle this problem. Recht and Ré [6] state a conjecture whose consequence is that RPCD converges faster than RCD on quadratic problems, but they prove the result only for some special random cases. Sun et al. [10] have analyzed the convergence speed of the distance between the expected iterate $\mathbb{E}[x^{k}]$ and the minimizer $x^{*}$ for convex quadratic problems, but this cannot be translated to a result for the expected squared error $\mathbb{E}\|x^{k}-x^{*}\|^{2}$ nor the expected function suboptimality $\mathbb{E}(f(x^{k})-f(x^{*}))$, which are much more informative quantities.111As an example of why a sequence $\{x^{k}\}$ for which $\mathbb{E}[x^{k}]=x^{*}$ does not give useful information about convergence rate, consider $x^{k}=x^{*}+r^{k}$, where $r^{k}$ are drawn i.i.d. from $N(0,I)$. Such a sequence has $\mathbb{E}[x^{k}]=x^{*}$, yet it has $\mathbb{E}\|x^{k}-x^{*}\|^{2}=1$, so cannot be said to converge to $x^{*}$ in expectation.

Computational experience reported in [11, 12] showed that for most convex quadratic functions (1.3), the convergence behaviors of all variants of CD are similar. For example, when $A$ is a matrix of the form $V\Sigma V^{T}$ where $V$ is random orthogonal and $\Sigma$ is a positive diagonal matrix whose diagonals (the eigenvalues of $A$) follow a log-uniform distribution, then CCD, RCD, and RPCD all converge at roughly the same rates, no matter how widely the eigenvalues are dispersed. However, these computational tests revealed a class of matrices $A$ for which the variants had radically different performance: matrices of the form

[TABLE]

for small positive values of $\delta$, and $\mathbf{1}=(1,1,\dotsc,1)^{T}$. For such matrices, the performance of RCD and RPCD is similar, but CCD converges much more slowly. In this paper, we explain much of this anomalous behavior by considering a matrix closely related to (1.10), and explaining the difference by means of a specialized analysis of RPCD.

The current paper is an extension of our paper [3] in which, motivated by the empirical observation above, we considered the special case of (1.10) in which $\Sigma=I$, that is,

[TABLE]

It was proved by [8] that this matrix achieves worst-case convergence behavior for CCD. We showed in [3] that a factor of $O(n^{2})$ fewer iterations are required by RPCD to achieve the same accuracy, and that the complexity of RPCD is similar to RCD in this case. Salient properties of the matrix (1.11) include the following.

• (a)

It has eigenvalue $\delta$ replicated $(n-1)$ times, and a single dominant eigenvalue $\delta+(1-\delta)n$, and

• (b)

it is invariant under symmetric permutations, that is, $P^{T}AP=A$ for all permutation matrices $P$.

The latter property makes the analysis of RPCD much more straightforward than for more general $A$ of the form (1.10). Specifically, it follows from (1.5) that $C_{P}\equiv C=-(L+\Delta)^{-1}L^{T}$, where $A=L+\Delta+L^{T}$, that is, $C_{P}$ is independent of $P$. For the matrix (1.11), the expression (1.6) thus simplifies to

[TABLE]

We refer to [3] for a more extensive discussion of prior related work on variants of coordinate descent. We note in particular that for general convex functions $f$, CCD has weaker convergence guarantees than for convex quadratic $f$, as analyzed in [1, 9, 4]. By contrast, the convergence results for RCD presented in [5] show no difference between quadratic and nonquadratic convex functions.

1.2. Contributions

In this work, we study the behavior of the RPCD variant of CD on problems of the form (1.3), where the coefficient matrix has the form

[TABLE]

for some $u\in\mathbb{R}^{n}$. This paper focuses on the case in which the components of $u$ are not too different in magnitude, and are all close to $1$. Rather than working directly with (1.12), we work with a diagonally scaled version that has a form more tractable for analysis. By scaling (1.12) symmetrically with the matrix $U=\operatorname{diag}(u)$, we obtain

[TABLE]

(Details are given in Section 2.) Note that both forms (1.12) and (1.13) are generalizations of (1.11). They are both closely related to the more general form (1.10), in that (1.12) can be obtained from (1.10) by a symmetric scaling with $\Sigma^{-1/2}V$, while (1.13) has the form (1.10) with $V=I$ and $\Sigma=I+(\epsilon/\delta)D$. Thus, this paper provides a significantly more complete explanation of the anomalous convergence behavior involving matrices (1.10) than our earlier work.

For matrices of the form (1.13) in (1.3), this paper proves similar convergence behavior for RPCD to what was proved in [3] for the special case (1.11), in the regime defined by the following values of the parameters $n$, $\epsilon$, and $\delta$:

[TABLE]

where $\rho_{1}$ is a positive or negative quantity of modest size, magnitude not much larger than $1$, and independent of $n$, $\epsilon$, and $\delta$. We prove that the convergence rate guarantee of RPCD for problems defined by (1.13) is similar to that of RCD, and much better than the rate bound for CCD. Specifically, we explain via analysis of a linear recurrence that captures the epoch-wise behavior of RPCD that the per-epoch objective improvement is bounded by a factor of approximately

[TABLE]

which is similar to the corresponding factors of approximately $1-\delta$ and $1-2\delta$ that are known for RCD (by different analyses), and significantly better than the factor of approximately $(1-\delta/n^{2})$ arising from worst-case theoretical guarantees for CCD. By the generalization of (1.11) to (1.13) (and thus (1.12)), we extend our understanding of the empirical behavior of RPCD, RCD, and CCD described at the beginning of Section 1.1.

1.3. Remainder of the Paper

In Section 2, we relate matrices of the forms (1.13) and (1.12), showing that the behavior of CD is similar on both. Section 3 presents our analysis for the behavior of RPCD on problem (1.3), (1.13). In particular, we define a sequence of matrices $\{\bar{A}_{\epsilon}^{(t)}\}$ such that given any initial guess $x^{0}$, the expected value of the objective $f(x^{tn})$ after the $t$th epoch is $\tfrac{1}{2}(x^{0})^{T}\bar{A}_{\epsilon}^{(t)}x^{0}$. We then define a sequence of matrices $\{\hat{A}_{\epsilon}^{(t)}\}$ that dominates $\{\bar{A}_{\epsilon}^{(t)}\}$, and that can be parametrized compactly. We analyze convergence of the sequence $\{\hat{A}_{\epsilon}^{(t)}\}$ by means of a spectral analysis of the matrix that relates its parameters at successive values of $t$, and use it to develop an estimate of the asymptotic per-epoch improvement of the objective $f(x^{tn})$, $t=0,1,2,\dotsc$. We provide an explanation in Section 3.5 for the large decrease in $f$ that is often observed in the very first iteration of CD, a phenomenon that is not explained by the asymptotic analysis. Section 4 discusses RCD and CCD variants for the problem (1.3), (1.13), while Section 5 reports computational experience with the three variants.

1.4. Notation

In addition to the notation $\rho_{1}$ mentioned above, which denotes a scalar quantity of size not much greater than $1$ and independent of $n$, $\epsilon$, and $\delta$, we make extensive use of vector quantities $\mathbf{r_{1}}\in\mathbb{R}^{n}$ and matrix quantities $\mathbf{R_{1}}\in\mathbb{R}^{n\times n}$ (symmetric in some contexts and nonsymmetric in others), which we assume are both bounded in norm by $1$, that is,

[TABLE]

In the case in which $\mathbf{R_{1}}$ is also symmetric, it follows from these assumptions that $-I\preceq\mathbf{R_{1}}\preceq I$. This notation is essential to capturing remainder terms that appear in our analysis. In particular, it allows us to keep explicit track of dependence of the remainder terms on $n$, $\epsilon$, and $\delta$. For example, a vector quantity whose size is bounded by a modest multiple of $\epsilon^{2}n^{-1}$ can be represented by $\rho_{1}\epsilon^{2}n^{-1}\mathbf{r_{1}}$. The following estimate follows immediately from this notation:

[TABLE]

Matrix and vector norms $\|\cdot\|$ signify $\|\cdot\|_{2}$ throughout, unless some other subscript is specified.

2. Quadratic functions with Hessians of the form (1.13)

We discuss here the matrix of the form (1.13), explaining its relationship to (1.12) and to (1.11), and giving some preliminaries for the analysis of RPCD on the corresponding quadratic function.

2.1. Relating (1.13) to (1.12)

Given $\epsilon>0$ and $\delta\in(0,1)$, suppose that $u\in\mathbb{R}^{n}$ satisfies

[TABLE]

Consider the matrix $B_{u}$ from (1.12). Defining $U:=\operatorname{diag}(u)$, we have

[TABLE]

and note that the diagonal elements of $U^{-2}$ are in the range $[1,\epsilon/\delta+1]$. Thus we can write $\delta U^{-2}=\delta I+\epsilon D$, where $D$ is diagonal with elements in $[0,1]$, so in fact $A_{\epsilon}$ in (2.2) has the form (1.13).

We verify in Appendix A that the iterates generated by Algorithm 1 for a given sequence of indices $i(\ell,j)$ to (1.3) with $A=B_{u}$ from (1.12), and with starting point $\tilde{x}^{0}$ and exact line search are isomorphic to the iterates generated by applying the same algorithm with the same index sequence to (1.3) with $A=A_{\epsilon}$ from (2.2), with starting point $x^{0}=U\tilde{x}^{0}$. Specifically, we have $x^{k}=U\tilde{x}^{k}$ for all $k\geq 0$, where $\{\tilde{x}^{k}\}$ is the iterate sequence corresponding to (1.12) and $\{x^{k}\}$ is the sequence corresponding to (2.2). Note that the function values coincide at each iteration, that is,

[TABLE]

Thus we expect to see similar asymptotic behavior for the quadratic objectives based on matrices (1.12) and (1.13), from starting points with the same distribution.

We note too that the matrix $A_{\epsilon}$ from (1.13) is “sandwiched” between scalar multiples of two matrices of the form (1.11). We have

[TABLE]

where $\delta^{\prime}=(\delta+\epsilon)/(1+\epsilon)$ and “$\leq$” denotes element-wise inequality. This observation suggests similar behavior for RPCD to that proved for the matrices (1.11) in [3]. Indeed, we observe similar behavior empirically, but we could not find a way to exploit the relationship (2.4) in our convergence analysis. The distinctiveness of the components of $D$ plays a key role; the effects of $D$ in (1.13) persist through the epochs. The analysis techniques in [3] make strong use of the fact that the epoch-wise iteration matrix $C_{P}$ defined in (1.5) is independent of $P$, a fact that no longer holds for matrices (1.13).

Representative numerical results for the three versions of CD on quadratics with Hessians of the form (1.11) are shown in Figure 1. We note here the nearly identical linear rates of the RPCD and RCD variants, and the much slower rate of the CCD variant. The same pattern is observed in Figure 2, which considers matrices of the forms (1.12) and (1.13). Note in particular that the latter two matrices are indistinguishable in their empirical behavior, further justifying our focus on the form (1.13) in our analysis.

2.2. RPCD Preliminaries

We now define some notation to be used in the remainder of the analysis: the matrix $C_{P}$ that defines the change in iterate $x$ over one epoch and the matrix $\bar{A}_{\epsilon}^{(\ell)}$ that defines the value $f(x^{\ell n})$ of the objective after $\ell$ epochs.

Applying to (1.13) the decomposition (1.4) into triangular and diagonal matrices, we obtain

[TABLE]

where

[TABLE]

Following (1.5), we have for $A_{\epsilon}$ that the epoch matrix is

[TABLE]

Our interest is in the quantity

[TABLE]

where $f(x^{\ell n})$ is defined by (1.8). Adapting notation from [3], we define the matrices $\bar{A}_{\epsilon}^{(t)}$, $t=0,1,2,\dotsc,\ell$ as follows:

[TABLE]

We have the following recursive relationship between successive terms in the sequence $\bar{A}_{\epsilon}^{(t)}$, $t=0,1,2,\dotsc$:

[TABLE]

where we have dropped the subscript on $P_{\ell-t+1}$ in the second equality, since the permutation matrices for each epoch are i.i.d. Using this matrix, we can compute (2.8) by

[TABLE]

3. Epoch-Wise Convergence of Expected Function Value

In this section, we analyze the behavior of the sequence of matrices $\{\bar{A}_{\epsilon}^{(t)}\}$ that govern the expected value of the objective function $f$ after $t$ epochs of RPCD. By focusing on the operation (2.9) which tracks the change from one element of this sequence to the next, we show that this sequence is bounded in norm by a quantity that decreases to zero at an asymptotic rate similar to the known rate for the fully-random variant RCD.

We show that the matrix sequence $\{\bar{A}_{\epsilon}^{(t)}\}$ is dominated222Given two symmetric matrices $F$ and $G$, we say that $F$ dominates $G$ if $F-G$ is positive semidefinite. by another sequence of positive definite matrices $\{\hat{A}_{\epsilon}^{(t)}\}$ that can be represented as a four-term recurrence

[TABLE]

where $\mathbf{r_{1}}$ is a vector such that $\|\mathbf{r_{1}}\|\leq 1$ (as defined in Section 1.4) and $(\hat{\eta}_{t},\hat{\nu}_{t},\hat{\epsilon}_{t},\hat{\tau}_{t})$ is a quadruplet of scalar coefficients for all $t=0,1,2,\cdots$. (Note that the quantities $\mathbf{r_{1}}$ in the final term are generally different for each $t$.) We set $\hat{A}_{\epsilon}^{(0)}=\bar{A}_{\epsilon}^{(0)}=A_{\epsilon}$, with

[TABLE]

and define the sequence $\{\hat{A}_{\epsilon}^{(t)}\}$ so that successive elements satisfy the same relationship as shown in (2.9) for $\{\bar{A}_{\epsilon}^{(t)}\}$, namely

[TABLE]

Our analysis consists chiefly of analyzing the convergence to zero of the sequence of quadruplets $\{(\hat{\eta}_{t},\hat{\nu}_{t},\hat{\epsilon}_{t},\hat{\tau}_{t})\}_{t=0,1,2\dotsc}$ corresponding to $\{\hat{A}_{\epsilon}^{(t)}\}_{1,2,\dotsc}$.

After several definitions and technical results in Section 3.1, we derive in Section 3.2 a tractable representation of the matrix $C_{P}$ from (1.5) that defines the transition between successive elements of the sequences $\{\bar{A}_{\epsilon}^{(t)}\}$ and $\{\hat{A}_{\epsilon}^{(t)}\}$. In Section 3.3, we examine the effect of the operation of $C_{P}$ on each of the four terms in the bounding sequence (3.1). In Section 3.4, we define the recurrence that relates successive elements of the sequence $\{(\hat{\eta}_{t},\hat{\nu}_{t},\hat{\epsilon}_{t},\hat{\tau}_{t})\}_{t=0,1,2\dotsc}$, and examine the rate at which this sequence converges to [math]. We show that the per-epoch rate is bounded by a scalar sequence that converges at a nearly linear rate of $1-1.4\delta$. (Our analysis is conservative; the true rate, observed in experiments, is often closer to $1-2\delta$.)

In most of this section, we consider the regime for parameters $n$, $\epsilon$, and $\delta$ defined by (1.14). The inequality $\delta\leq\epsilon$ is made mostly for convenience; it implies that we can replace $\delta$ by $\epsilon$ in remainder terms, and it allows wide divergence in the diagonal elements of the matrix (1.13). (We expect that the main convergence results will continue to apply in a regime in which $0\leq\epsilon<\delta$, which indeed is closer to the matrix (1.11) studied in [3], which has constant diagonals, but the remainder terms in the analysis will need to be handled differently.) In the analysis of Section 3.4, we make additional assumptions on $n$, $\epsilon$, and $\delta$.

3.1. Definitions and Technical Results

We start by defining some useful quantities, drawing on [3], and proving several elementary results. While technical, these results give an idea of the effects of applying expectations over permutations to matrices that arise in the subsequent analysis.

From (1.13) and (3.3), we have

[TABLE]

From the definition of $D$ in (1.13), we have $d_{\textrm{av}}\in(0,1)$ and $d_{\textrm{av,2}}\in(0,1)$. We use $\pi$ to denote the permutation of $\{1,2,\dotsc,n\}$ associated with the permutation matrix $P$, so that for any vector $u\in\mathbb{R}^{n}$, we have

[TABLE]

We can see immediately that

[TABLE]

A useful conditional probability is as follows:

[TABLE]

This claim follows because $Pe_{2}$ contains $n-1$ zeros and a single $1$, and the $1$ cannot appear in position $i$ (because $Pe_{1}=e_{i}$) but may appear in any other position with equal likelihood.

A quantity that appears frequently in the analysis is the matrix $F$ defined by

[TABLE]

that is, the $n\times n$ matrix of all zeros except for $1$ on the diagonal immediately above the main diagonal. We see immediately that $\|F\|=1$. Several identities follow:

[TABLE]

We also have

[TABLE]

To verify this claim, note that the diagonals of $PFP^{T}$ are zero for all permutation matrices $P$, while the off-diagonals are $1$ with equal probability. Thus the expected value of the $n(n-1)$ off-diagonal elements is obtained by distributing the $n-1$ nonzeros in $F$ with equal weight among all off-diagonal elements, giving an expected value of $1/n$ for each of these elements, as in (3.9).

We have the following results about quantities involving $F$.

Lemma 3.1**.**

[TABLE]

Proof.

For (3.10a), we see that $P^{T}DPe_{1}$ is a multiple of $e_{1}$, and that $Fe_{1}=0$.

For (3.10b), we use $\mathbb{E}_{i}$ to denote the expectation with respect to index $i$ uniformly distributed over $\{1,2,\dotsc,n\}$, and recall that $\pi$ denotes the permutation corresponding to $P$. We have

[TABLE]

as required. ∎

Finally, we make frequent use of the following trivial result about the norm of rank-1 matrices: for any vectors $v,w\in\mathbb{R}^{n}$, we have

[TABLE]

In particular, we have from $\|\mathbf{1}\|=n^{1/2}$ that

[TABLE]

and in particular, using the notation of Section 1.4, we have

[TABLE]

3.2. Properties of the Epoch Matrix $C_{P}$

As in [3], we define

[TABLE]

We noted in [3] that

[TABLE]

so by using notation (3.7), we have

[TABLE]

We have further from a standard matrix-norm inequality together with the facts that $\|\bar{L}\|_{1}\leq 2$ and $\|\bar{L}\|_{\infty}\leq 2$ that

[TABLE]

Moreover, from [3, Section 2.2], we have

[TABLE]

so that

[TABLE]

(The validity of the remainder term in this expression follows from the fact that the coefficients of $\delta^{2},\delta^{3},\dotsc,\delta^{n-1}$ in $\bar{L}E^{T}$ all have the form $\mathbf{R_{1}}+\mathbf{r_{1}}\mathbf{1}^{T}$, so we can absorb them all into a single term of order $\delta^{2}$ by summation.)

The following lemma provides a useful estimate of the epoch matrix $C_{P}$.

Lemma 3.2**.**

Suppose that (1.14) holds. Then for $C_{P}$ defined by (1.5) and (2.7), we have

[TABLE]

Proof.

Note first that for a matrix $Y$ with $\|Y\|\leq\rho_{1}$ and for $\epsilon$ satisfying (1.14), we have

[TABLE]

From (2.7), using definition (3.14), we have

[TABLE]

By substituting from (3.17) and (3.19) (noting that $\|\bar{L}D_{P}\|\leq\|\bar{L}\|\leq 2$ from (3.16)), we have

[TABLE]

where we used $\delta\leq\epsilon$ from (1.14) to absorb the term $\delta^{2}(\rho_{1}\mathbf{R_{1}}+\rho_{1}\mathbf{r_{1}}\mathbf{1}^{T})$. The result follows immediately when we use (3.15) to substitute for $\bar{L}$, and again use $\delta\leq\epsilon$ together with $\|D_{P}\|\leq 1$ and (1.17) to absorb the remainder terms. ∎

3.3. Single-Epoch Analysis

In this section we analyze the change in each term in the expression (3.1) over a single epoch. We examine in turn the following terms:

• •

the $I$ term: Lemma 3.3,

• •

the $D$ term: Lemma 3.4,

• •

the $\mathbf{1}\mathbf{1}^{T}$ and $(\mathbf{r_{1}}\mathbf{1}^{T}+\mathbf{1}\mathbf{r_{1}}^{T})$ terms: Lemma 3.5.

Proofs of these technical results appear in Appendix B.

Lemma 3.3**.**

Suppose that (1.14) holds. We have

[TABLE]

Lemma 3.4**.**

Suppose that (1.14) holds. We have

[TABLE]

Lemma 3.5**.**

Suppose that (1.14) holds. For any $v\in\mathbb{R}^{n}$, we have

[TABLE]

so that

[TABLE]

When $v=\mathbf{1}$, we have

[TABLE]

The following result summarizes Lemmas 3.3, 3.4, and 3.5, using the assumption $n\epsilon\leq 1$ from (1.14) to simplify some terms.

Theorem 3.6**.**

Suppose that (1.14) holds. We have

[TABLE]

Proof.

The first result (3.23a) follows immediately from Lemma 3.3 when we note that $\epsilon^{2}=n^{-1/2}(\epsilon n^{-1/2})(\epsilon n)\leq\epsilon n^{-1/2}$. The bound (3.23b) is immediate from (3.22) in Lemma 3.5. Lemma 3.4 immediately yields (3.23c). For (3.23d), we use $\epsilon n\leq 1$ and $\epsilon^{2}n^{1/2}=\epsilon(n\epsilon)n^{-1/2}\leq\epsilon n^{-1/2}$ to simplify the coefficients of $I$ and $\mathbf{1}\mathbf{1}^{T}$ in (3.21). ∎

3.4. The Four-Term Recurrence and Convergence Bound for RPCD

In this section we discuss the sequence of $n\times n$ symmetric matrices $\hat{A}_{\epsilon}^{(t)}$ that dominates the sequence $\bar{A}_{\epsilon}^{(t)}$ defined in Section 2.2. Using the results of the previous subsection, together with the four-term parametrization of $\hat{A}_{\epsilon}^{(t)}$ defined in (3.1), we derive a recurrence relationship for the sequence of quadruplets $\{(\hat{\eta}_{t},\hat{\nu}_{t},\hat{\epsilon}_{t},\hat{\tau}_{t})\}_{t=0,1,2,\dotsc}$. By finding the rate at which this sequence decreases to zero, we derive a bound on the expected values of $f$ after each epoch of RPCD.

We now show the main result for recurrence of the representation (3.1).

Theorem 3.7**.**

Suppose that (1.14) holds. Consider a nonnegative sequence of quadruplets $\{(\hat{\eta}_{t},\hat{\nu}_{t},\hat{\epsilon}_{t},\hat{\tau}_{t})\}_{t=0,1,2,\dotsc}$ satisfying

[TABLE]

along with the recurrence

[TABLE]

where

[TABLE]

where each $\rho_{1}$ represents a positive quantity not much greater than $1$ and independent of $n$, $\epsilon$, and $\delta$. Then we have for $\hat{A}_{\epsilon}^{(t)}$ defined by (3.1) that $\hat{A}_{\epsilon}^{(t)}\succeq\bar{A}_{\epsilon}^{(t)}$ for all $t$.

Proof.

By definition, we have $\hat{A}_{\epsilon}^{(0)}\succeq\bar{A}_{\epsilon}^{(0)}$. Supposing that $\hat{A}_{\epsilon}^{(t)}\succeq\bar{A}_{\epsilon}^{(t)}$ for some $t\geq 0$, we have from (2.9) that

[TABLE]

Analogous to (3.1), we define the following matrix, parametrized by the coefficients $(\tilde{\eta}_{t+1},\tilde{\nu}_{t+1},\tilde{\epsilon}_{t+1},\tilde{\tau}_{t+1})$ defined in (3.25):

[TABLE]

Since $(\hat{\eta}_{t},\hat{\nu}_{t},\hat{\epsilon}_{t},\hat{\tau}_{t})\geq 0$, we can use Theorem 3.6 to ensure that

[TABLE]

A little more explanation is needed here. Because the matrices $I$, $\mathbf{1}\mathbf{1}^{T}$, and $D$ (the coefficients of $\hat{\eta}_{t}$, $\hat{\nu}_{t}$, and $\hat{\epsilon}_{t}$, respectively) are positive semidefinite, we can use the upper bounds in (3.23a), (3.23b), and (3.23c) to derive the $\preceq$ relationship. The coefficient of $\hat{\tau}_{t}$ may not be positive definite, but since $\hat{\tau}_{t}\geq 0$, we can still use the bound (3.23d) to establish the $\preceq$ relationship. Moreover, we can assume that $\tilde{\tau}_{t+1}\geq 0$, by replacing $\mathbf{r_{1}}$ by $-\mathbf{r_{1}}$ in the representation (3.28) if necessary. Thus, from (3.25), we have

[TABLE]

By combining this expression with (3.27) and (3.29), we obtain $\bar{A}_{\epsilon}^{(t+1)}\preceq\hat{A}_{\epsilon}^{(t+1)}$, as required. ∎

We now analyze the decay of the sequence of quadruplets $(\hat{\eta}_{t},\hat{\nu}_{t},\hat{\epsilon}_{t},\hat{\tau}_{t})$ generated by this recursion. For purposes of this analysis, we assume that all quantities $\rho_{1}$ that appear in the matrix $\hat{M}$ defined by (3.26) are bounded in magnitude by constant $\bar{\rho}$. From this constant, we define

[TABLE]

We place further restrictions on the allowable regime for values of $n$, $\epsilon$, and $\delta$, in addition to those in (1.14). Specifically, we require

[TABLE]

As immediate consequences of these bounds, in combination with (1.14) and (3.30), we have

[TABLE]

Other useful consequences of (1.14) and (3.31), used repeatedly below, are as follows

[TABLE]

We now define two sequences that can be used to bound in norm the quadruplets $(\hat{\eta}_{t},\hat{\nu}_{t},\hat{\epsilon}_{t},\hat{\tau}_{t})$. These are

[TABLE]

We note immediately by combining with (3.33a) and (3.33b) that

[TABLE]

The following lemma details how the sequence of quadruplets $(\hat{\eta}_{t},\hat{\nu}_{t},\hat{\epsilon}_{t},\hat{\tau}_{t})$ is bounded in terms of the quantities in (3.35). Its proof appears in Appendix C.

Lemma 3.8**.**

Assume that the conditions (1.14) and (3.31) hold, and let $(\hat{\eta}_{t},\hat{\nu}_{t},\hat{\epsilon}_{t},\hat{\tau}_{t})$ be defined as in Theorem 3.7, and $\bar{\eta}_{t}$ and $\bar{\epsilon}_{t}$ be defined as in (3.34). Then the following bounds hold for all $t=1,2,\dotsc$:

[TABLE]

We are now ready to prove the main convergence result.

Theorem 3.9**.**

Suppose that the RPCD version of Algorithm 1 is applied to function $f$ defined by (1.3) with coefficient matrix satisfying (1.13). Suppose that the quantities $\rho_{1}$ in each recurrence matrix $\hat{M}$ in (3.25) are all bounded in magnitude by $\bar{\rho}$, and that conditions (1.14) and (3.31) hold. Then there is a constant $C$ such that for all $t=1,2,\dotsc$, we have

[TABLE]

indicating an asymptotic per-epoch convergence rate approaching $1-1.4\delta$.

Proof.

The proof follows from (2.10) and (3.34) when we use Lemma 3.8, the bound $\|\hat{A}_{\epsilon}^{(t)}\|\leq\bar{C}\|(\hat{\eta}_{t},\hat{\nu}_{t},\hat{\epsilon}_{t},\hat{\tau}_{t})\|$ for some $\bar{C}>0$, and the bound

[TABLE]

for some $\hat{C}>0$, where we used $\delta\leq\epsilon$ in the last step. The final claim follows by taking the ratio of the bound after $t+1$ and $t$ epochs, which approaches $(1-1.4\delta)$ as $t\to\infty$. ∎

3.5. Decrease in the First Iteration

A behavior of all CD variants that we observe in Figures 1-2 is that the objective value decreases dramatically in the very first iteration of the algorithm. The theorem below shows that this phenomenon can be explained for both RCD and RPCD by using an extension of the analysis in [3, Theorem 3.4]. (Similar reasoning also applies for CCD in most cases, but there is no guarantee, since adversarial examples consisting of particular choices of $x^{0}$ can be constructed.) Geometrically, the phenomenon is due to the function (1.3), (1.13) increasing rapidly along just one direction — the all-one direction $\mathbf{1}$ — and more gently in other directions. Thus an exact line search along any coordinate search direction will identify a point near the bottom of this multidimensional “trench.”

Our result for first-iteration decrease is as follows.

Theorem 3.10**.**

Consider solving (1.3) with the matrix $A=A_{\epsilon}$ defined in (1.13), and $\epsilon\in(0,1)$, using CCD, RCD, or RPCD with exact line search. Then after a single iteration, we have

[TABLE]

where $i=i(0,0)$ is the coordinate chosen for updating in the first iteration. When RCD or RPCD is used, we further have that

[TABLE]

Proof.

Suppose that $i\in\{1,2,\dotsc,n\}$ is the component chosen for updating in the first iteration, which is chosen uniformly at random from $\{1,2,\dotsc,n\}$ for RPCD and RCD. After a single step of CD, we have

[TABLE]

Thus, from (1.13), we have

[TABLE]

Since $d_{i}\in[0,1]$, $\delta\in(0,1)$, and $\epsilon\in(0,1)$, it can be shown that

[TABLE]

Thus by substitution into (3.39), we obtain

[TABLE]

Further, by noting that

[TABLE]

the desired result (3.37) is obtained.

The result (3.38) is then obtained by noting that $d_{i}\leq 1$ and that

[TABLE]

whose derivation can be found in the proof of [3, Theorem 3.4]. ∎

We can compare $f(x^{1})$ from this theorem with $f(x^{0})$ obtained by substituting into (1.13), which is

[TABLE]

Note that the second term, which involves $(\mathbf{1}^{T}x^{0})^{2}$, decreases by a factor of approximately $(\delta+\epsilon)$ in the first iteration, whereas the first term, which involves $\|x^{0}\|^{2}$, does not change much from its original value, which is typically already small. For most starting points, the decrease is dramatic.

4. Analysis for CCD and RCD

The analysis for the RCD variant of coordinate descent for (1.3), (1.13) follows from the standard analysis [5]. The modulus of convexity $\mu$ is $\delta$, while the maximum coordinate-wise Lipschitz constant for the gradient $L_{\max}$ is $1+\epsilon$. The per-epoch linear rate of expected improvement in $f$ for RCD on $A_{\epsilon}$ is thus

[TABLE]

yielding a complexity of $O(|\log\hat{\epsilon}|/(\delta(1-\epsilon)))$ iterations for reaching an $\hat{\epsilon}$-accurate objective. The (slightly tighter) complexity of RCD from [5, Section 4] improves this epoch bound by approximately a factor of $2$, to

[TABLE]

That is, the per-epoch convergence rate of a bound on $\rho_{\text{RCD}}$ is approximately $1-2\delta/(1+\epsilon+\delta)$.

It is also shown in [5] that one can get an improved rate by non-uniform sampling of the coordinates when the coordinate-wise Lipschitz constants are not identical. In particular, for (1.3), (1.13), if the probability that the $i$th coordinate is sampled is proportional to the value of $(A_{\epsilon})_{ii}$, the result in [5] improves $\rho_{\text{RCD}}$ to

[TABLE]

Since $d_{\textrm{av}}\in(0,1)$, this improvement is rather insignificant, and the rate is still worse than that of (4.2). (Whether nonuniform sampling can improve the complexity expression (4.2) is unknown.)

For CCD, we note that the iterates have the form

[TABLE]

where $C=-(L+\Delta)^{-1}L^{T}$, where $A=L+\Delta+L^{T}$ is the triangular-diagonal splitting of $A$ (that is, $C=C_{I}$ from (1.4), (1.5)). Thus

[TABLE]

and the asymptotic behavior of the sequence of function values is governed by $\|C^{\ell}\|^{2}$. By Gelfand’s formula [2], the asymptotic per-epoch decrease factor is thus approximately $\rho(C)^{2}$. Proposition 3.1 of [8] yields an upper bound on the per-epoch decrease factor. Noting that the largest eigenvalue of $A_{\epsilon}$ is bounded above by $n(1-\delta)+\delta+\epsilon$, their bound is as follows:

[TABLE]

which is approximately $1-\delta/n^{2}$ for the ranges of values of $\delta$ and $\epsilon$ of interest in this paper. The implied iteration complexity guarantee is about a factor of $n^{2}$ worse than that for RCD. In our computational experiments, we compare empirical observations of CCD convergence rate with $\rho(C)^{2}$ rather than with (4.3).

Note that the upper bounds for convergence rates of RCD and CCD are worst-case guarantees. On the problem class (1.13), we show that the convergence rate of RPCD is similar to the bound for RCD, and we see in the next section that both bounds are quite tight in practice. The worst-case bounds on CCD are looser, in the sense that the computational behavior is not quite as poor as these bounds suggest. Nevertheless, comparison of the worst-case bounds correctly foreshadows that relative behavior of the different variants on these problems, seen in Figure 2: CCD is much slower than RCD or RPCD on this class of problems.

5. Computational Results

We report here on some experiments with variants of CD on problems of the form (1.3), (1.13). Fixing $n=100$, we tried different settings of $\epsilon$ and $\delta$, and ran the three variants CCD, RCD, and RPCD for many epochs. Results are reported in Tables 1 and 2. We obtain empirical estimates of the per-epoch asymptotic convergence rate by geometrically averaging the rate over the last 10 epochs, tabulating these observations as $\rho_{\text{CCD}}(\delta,\text{observed})$, $\rho_{\text{RCD}}(\delta,\text{observed})$, and $\rho_{\text{RPCD}}(\delta,\text{observed})$. Since we report the difference between these quantities and $1$ in the tables, larger numbers correspond to faster rates. (The numbers in the table are reported in scientific notation, with $a(b)$ representing $a\times 10^{b}$.) As noted in Section 4, we use $\rho(C)^{2}$ as the theoretical bound on the convergence rate for CCD, while we use $\rho_{\text{RCD}}$ from [5, Section 4] (which corresponds to the complexity (4.2)) as the theoretical estimate of the convergence rate for RCD. For RPCD, we used $2\delta$ as a “benchmark” value, corresponding to a per-epoch rate of $1-2\delta$, slightly faster than the $1-1.4\delta$ rate proved in Section 3.

Not all the settings of parameters $n$, $\epsilon$, and $\delta$ in these tables satisfy the conditions (1.14), (3.31) that were assumed in our analysis. We mark with an asterisk those entries for which these conditions are not satisfied. We note that the benchmark rate of $1-2\delta$ continues to hold in regimes beyond the reach of our theory. This accords with the observation that the matrix $\hat{M}$ defined in (3.26) indeed has norm very close to $1$, so that behavior of the sequence is governed chiefly by the $(1-\delta)^{2}$ factor in the recurrence (3.25).

These tables confirm that the empirical performance of RCD and RPCD is quite similar, across a wide range of parameter values, and markedly faster than CCD.

Acknowledgments

We are grateful for the careful reading and penetrating comments of two referees, which resulted in significant improvement of the results.

Appendix A Invariance of Coordinate Descent under Diagonal Scaling

Coordinate descent applied to quadratics (1.3) with exact line search at each iterate is invariant under symmetric diagonal scalings of $A$. For any symmetric positive definite $A$ and nonzero diagonal $F$, define

[TABLE]

Note that $\tilde{A}$ is symmetric positive definite. Consider the objective functions (1.3) defined with Hessians $A$ and $\tilde{A}$. For a given $x^{0}$, define $\tilde{x}^{0}=Fx^{0}$. The function values match at these points, that is,

[TABLE]

Considering the iterates generated by Algorithm 1 for the two functions, with $\alpha_{k}$ defined by exact line searches, and the same choices of coordinates $i(\ell,j)$ at each iteration. Assume that $Fx^{t}=\tilde{x}^{t}$ for $t=1,2,\dotsc,k$. Suppose that coordinate $i$ is chosen at iteration $k$, the updates are

[TABLE]

By noting that

[TABLE]

and using the inductive hypothesis, it is easy to verify that $\tilde{x}^{k+1}=Fx^{k+1}$, as required.

Appendix B Proofs of Lemmas from Section 3.3

B.1. Proof of Lemma 3.3

Proof.

From Lemma 3.2, we have

[TABLE]

(We give further details on the $\epsilon^{2}$ term below.) For the $O(1)$ term in (B.1b), we have from (3.5) and $e_{1}^{T}e_{1}=1$ that

[TABLE]

For the first part of the $O(\delta)$ term, we have from (3.5), (3.8), (3.9), and $e_{1}^{T}e_{2}=0$ that

[TABLE]

(By symmetry, the second part of the $O(\delta)$ term will also have expectation $-\frac{1}{n}I$.)

For the $O(\epsilon)$ term in (B.1b), we have from (2.6), (3.8), (3.9), Lemma 3.1, and the fact that $\mathbb{E}_{P}e_{1}^{T}D_{P}e_{1}=d_{\textrm{av}}$ that

[TABLE]

The lower-order terms in the main result follows by substituting this estimate along with (B.2) and (B.3) into (B.1b).

We now address the $\epsilon^{2}$ term in (B.1b). Gathering together all terms with coefficients $\epsilon^{2}$, $\epsilon\delta$, and $\delta^{2}$ from (B.1a), we have

[TABLE]

The first term in this expression (and its transpose) is clearly accounted for by the $\epsilon^{2}$ term in (B.1b). For the other $\epsilon^{2}$ term, and also the $\epsilon\delta$ terms, we use the facts that $\|F\|=1$ and $\|D_{P}(I-F)\|\leq\|D_{P}\|(\|I\|+\|F\|)=\rho_{1}$, along with $\delta\leq\epsilon$, to deduce that these terms too are accounted for by the $\epsilon^{2}$ term in (B.1b). From $\|F\|=1$ and $e_{2}^{T}e_{2}=1$, we can say the same too for the coefficient of $\delta^{2}$.

For the final “$\preceq$” claim in the lemma, we use the facts (1.16), $d\mathbf{1}^{T}=\rho_{1}n^{1/2}\mathbf{r_{1}}\mathbf{1}^{T}$ (from (3.11)), $d_{\textrm{av}}\in(0,1]$, and $D\succeq 0$, ∎

B.2. Proof of Lemma 3.4

Proof.

From Lemma 3.2, we have

[TABLE]

where we used (3.10a) from Lemma 3.1 along with $e_{2}^{T}D_{P}e_{1}=0$ to simplify the coefficient of $\delta$. (Further justification for the form of the $\epsilon^{2}$ term appears below.)

For the $O(1)$ term, we have that

[TABLE]

Thus from (3.5), we have by taking expectations over $P$ that

[TABLE]

as required.

For the coefficient of $\delta$, we have

[TABLE]

Taking expectations with respect to $D$, we have from (3.5) and (3.9) that

[TABLE]

For the coefficient of $\epsilon$, we have

[TABLE]

where we used (3.10a) in Lemma 3.1 to eliminate terms that are multiples of $Fe_{1}=0$.

For the first term in (B.5), we have from (3.9) that

[TABLE]

For the second term in (B.5), we have

[TABLE]

where we used (3.10b) from Lemma 3.1 and the definition of $d_{\textrm{av}}$ in (3.3). The third term in (B.5) is the transpose of this second term. For the final term in (B.5), we have

[TABLE]

By substituting (B.6), (B.7), and (B.8) into (B.5), we obtain the required coefficient of $\epsilon$.

We return to verifying the form of the $\epsilon^{2}$ term in (B.4b). The coefficients of $\epsilon^{2}$, $\epsilon\delta$, and $\delta^{2}$ terms from (B.4a) are are follows:

[TABLE]

By making use of the bounds $\|I\|=\|F\|=1$, $\|D_{P}\|\leq 1$, $\|e_{1}\|=\|e_{2}\|=1$, and $\delta\leq\epsilon$, we see that this expression is accounted for by the coefficient of $\epsilon^{2}$ in (B.4b).

For the final “$\preceq$” relationship, we use $dd^{T}\preceq nI$ to obtain $(2/n)\epsilon dd^{T}\preceq 2\epsilon I$, $\mathbf{1}\mathbf{1}^{T}D^{2}=n^{1/2}\mathbf{1}\mathbf{r_{1}}^{T}$ to bound the terms with $\mathbf{1}\mathbf{1}^{T}D^{2}$ (and similarly for $D^{2}\mathbf{1}\mathbf{1}^{T}$), $D^{2}\succeq 0$, $\mathbf{1}d^{T}=n^{1/2}\mathbf{1}\mathbf{r_{1}}^{T}$ to obtain $-(1/n)\mathbf{1}d^{T}\preceq n^{-1/2}\mathbf{1}\mathbf{r_{1}}^{T}$, and $\mathbf{R_{1}}\preceq I$. ∎

B.3. Proof of Lemma 3.5

Proof.

We have

[TABLE]

where we use $\operatorname{(transpose)}$ to denote the transpose of the explicitly stated terms. From Lemma 3.2 and $P^{T}\mathbf{1}=\mathbf{1}$, we have

[TABLE]

To derive the remainder term (the coefficient of $\epsilon^{2}$ in (B.9b)), we need to consider the coefficients of $\epsilon^{2}$, $\delta\epsilon$, and $\delta^{2}$ from (B.9a). The coefficient of the $\epsilon^{2}$ term is

[TABLE]

Using $(I-\mathbf{1}e_{1}^{T})\mathbf{1}=\mathbf{1}-\mathbf{1}=0$, we see that the first term in this expression vanishes. From (3.7) and (3.8), we have several other identities:

[TABLE]

In the third term, we thus have that $(-D_{P}+D_{P}F)\mathbf{1}=-D_{P}(I-F)\mathbf{1}=-D_{P}e_{n}=\rho_{1}\mathbf{r_{1}}$. We also have that $e_{1}^{T}D_{P}e_{n}=0$ and $v^{T}(-D_{P}+F^{T}D_{P})=\rho_{1}\|v\|\mathbf{r_{1}}^{T}$. Thus (B.10) becomes

[TABLE]

which is accounted for by the $\epsilon^{2}$ term in (B.9b). We turn next to the coefficient of $\epsilon\delta$ in (B.9a). This term consists of the following expression plus its transpose:

[TABLE]

Because $\delta\leq\epsilon$, this term (plus its transpose) can also be accounted for by the $\epsilon^{2}$ term in (B.9b). For the coefficient of $\delta^{2}$ in (B.9a), we have, using (B.11) again,

[TABLE]

which can also be absorbed into the $\epsilon^{2}$ term in (B.9b).

Returning to the lower-order terms in (B.9b), we use again the fact that $(I-\mathbf{1}e_{1}^{T})\mathbf{1}=0$ to eliminate the $O(1)$ term, and also one of the two terms in the coefficients of both $\epsilon$ and $\delta$. We thus obtain

[TABLE]

Additionally, we have from (2.6) that

[TABLE]

By substituting these identities into (B.12), we obtain

[TABLE]

where the second equality follows from

[TABLE]

since $D_{P}$ is diagonal for all $P$.

Taking expectations, we have for the coefficient of $(-\delta)$ in (B.13) that

[TABLE]

where the second equality is from a conditional expectation over permutation matrices $P$ such that $Pe_{1}=j$ and $Pe_{n}=i$, for all $i,j=1,2,\dotsc,n$ with $i\neq j$. By combining (B.14) with its transpose, we obtain the full coefficient of $(-\delta)$ in (3.20), which is

[TABLE]

This verifies the $O(\delta)$ term in (3.20).

We note that the coefficient of $(-\epsilon)$ in (B.13) is the same as the coefficient of $(-\delta)$, except for being multiplied from the left by $D$, which is independent of $P$. Thus the expectation of this term is simply (B.14), multiplied from the left by $D$, that is,

[TABLE]

We obtain the full coefficient of $(-\epsilon)$ in (3.20) by adding this quantity to its transpose, to obtain

[TABLE]

as required.

For the “$\preceq$” result (3.21), we use

[TABLE]

together with $\delta\leq\epsilon$ from (1.14). We also use $\mathbf{r_{1}}\mathbf{1}^{T}+\mathbf{1}\mathbf{r_{1}}^{T}=n^{1/2}\mathbf{R_{1}}$ and $-I\preceq\mathbf{R_{1}}\preceq I$ for symmetric $\mathbf{R_{1}}$ to absorb the $\epsilon^{2}(\mathbf{r_{1}}\mathbf{1}^{T}+\mathbf{1}\mathbf{r_{1}}^{T})$ term in (3.20) into the $I$ term in (3.21).

For (3.22), we have

[TABLE]

where we used the following identities for the third equality:

[TABLE]

and $e_{1}^{T}D_{P}e_{n}=0$ for the fourth equality. By substituting $D_{P}=P^{T}DP$ into (B.15), we obtain

[TABLE]

By taking the outer product of this vector with itself, and using $\delta\leq\epsilon$, we obtain

[TABLE]

where in the remainder term we used $\|Pe_{n}\|=1$ and $\|D\|\leq 1$. By taking expectations over $P$, and using $\mathbb{E}_{P}(Pe_{n})(Pe_{n})^{T}=n^{-1}I$, we obtain

[TABLE]

where we used (1.14) in the last expression to deduce that $\epsilon^{3}n\leq\epsilon^{2}$.

∎

Appendix C Proof of Lemma 3.8

Proof.

Note first that $\hat{\eta}_{t}$, $\hat{\nu}_{t}$, $\hat{\epsilon}_{t}$, and $\hat{\tau}_{t}$ are all nonnegative by definition. In this proof, we use repeatedly that they can be bounded by $|\tilde{\eta}_{t}|$, $|\tilde{\nu}_{t}|$, $|\tilde{\epsilon}_{t}|$, and $|\tilde{\tau}_{t}|$, respectively, though the $|\cdot|$ are unnecessary in the case of $\tilde{\epsilon}_{t}$ (since its exact value can be determined trivially from (3.25)) and in the case of $\tilde{\tau}_{t}$ (which can be assumed without loss of generality to be nonnegative, as mentioned in the proof of Theorem 3.7).

The proof is by induction on $t$. We show first that the bounds (3.36) hold for $t=1$.

We have from (3.33b) and the obvious property $\hat{\epsilon}_{t}=(1-\delta)^{2t}\epsilon$ from (3.25) that

[TABLE]

verifying (3.36b) for $t=1$. For (3.36a), we have from (3.25) with $t=0$, using the initial values (3.24) and the bounds in (3.32) that that

[TABLE]

It follows from $\hat{\rho}\geq 3$ and $\delta\leq.2$ (see (3.32b)) that

[TABLE]

verifying (3.36a) for $t=1$. For (3.36c) with $t=1$, we have

[TABLE]

where for the second inequality we used $n^{-1/2}\leq.5$ and $\epsilon^{2}\leq.04$. Continuing, we use $\bar{\eta}_{1}\geq 4(1-1.4\delta)\delta$ and $\bar{\epsilon}_{1}=(1-1.8\delta)\epsilon$ to write

[TABLE]

which suffices to prove (3.36c) for $t=1$. For (3.36d), we simply use $\epsilon\leq.2$ from (3.32).

For (3.36e) with $t=1$, we have from (3.25) and (3.24), using again $\epsilon^{2}\leq.04$ from (3.32), as well as $d_{\textrm{av}}\leq 1$ from (3.3) that

[TABLE]

which suffices to demonstrate (3.36e) for $t=1$.

Assuming now that (3.36) holds for some $t\geq 1$, we prove that the bounds holds for $t+1$ as well. We start with (3.36c). It follows from (3.25) that

[TABLE]

It then follows from (3.35) that

[TABLE]

as required. As earlier, (3.36d) follows immediately when we note that $\epsilon\leq.2$, from (3.32).

For (3.36e), we have

[TABLE]

where we used (3.36d) for the second inequality. Using now the bounds $\bar{\rho}\epsilon^{2}\leq.05$ and $n^{-1/2}\leq.5$ (from (3.32)), $d_{\textrm{av}}\leq 1$, and $\epsilon n^{-1/2}=(\epsilon n)n^{-3/2}\leq n^{-3/2}\leq.1$, we have

[TABLE]

so that

[TABLE]

as required.

The proof for (3.36b) is trivial, since

[TABLE]

We now prove (3.36a) for $t$ replaced by $t+1$. We have, substituting from the other formulas in (3.36), and using the bounds in (3.32), that

[TABLE]

where we used the definition (3.30) of $\hat{\rho}$ for the final inequality. Thus from (3.33), substituting from (3.34), and using $\epsilon^{2}<\delta$ from (1.14), we have

[TABLE]

as required. This completes the inductive step and hence the proof. ∎