On the Efficiency of Random Permutation for ADMM and Coordinate Descent

Introduction

A simple yet powerful idea for solving large-scale computational problems is to iteratively solve smaller subproblems. The applications of this idea include coordinate descent (CD), POCS (Projection onto Convex Sets), SGD (Stochastic Gradient Descent). They are well suited for large-scale unconstrained optimization problem (see, e.g. Wright , for a recent survey of CD) since it decomposes a large problem into small subproblems. The decomposition idea is crucial for huge problems due to both the cheap per-iteration cost and small memory requirement. Moreover, this idea is “orthogonal” to other large-scale optimization ideas such as first-order methods (using only gradient information) and random projection, and thus can be easily combined with other ideas.

This paper is motivated by a natural question: how should we extend the decomposition idea to solve problems with constraints? We consider a constrained minimization problem with a convex objective function and linear constraints (this is for motivation; our analysis is for a much simpler version):

To apply the decomposition idea to a constrained problem, one possible way is to form the augmented Lagrangian function and perform coordinate descent for the primal problem and a gradient step for the dual problem, i.e. combining BCD with augmented Lagrangian method, to obtain the so-called alternating direction method of multipliers (ADMM). ADMM was originally proposed in Glowinski and Marroco (see also Chan and Glowinski , Gabay and Mercier ) to solve problem (1) when there are only two blocks (i.e. $n=2$) and the objective function is separable. It is natural and computationally beneficial to extend the original ADMM directly to solve the general $n$-block problem (1) via the following procedure:

The convergence of the direct extension of ADMM to multi-block case had been an open question, until a counter-example was recently given in Chen et al. . More specifically, Chen et al. showed that even for the simplest scenario where the objective function is and the number of blocks is $3$, ADMM can be divergent for a certain choice of $A=[A_{1},A_{2},A_{3}]$. There are several proposals to overcome the drawback (see, e.g., ), but they either need to restrict the range of original problems being solved, add additional cost in each step of computation, or limit the stepsize in updating the Lagrange multipliers. These solutions typically slow down the performance of ADMM for solving most practical problems. Moreover, it is not clear how to compare the convergence speed of these algorithms as they typically contain different parameters. One may ask whether a “minimal” modification of cyclic multi-block ADMM (2) can lead to convergence, and whether we can provide some convergence speed analysis that is easy to interpret.

One of the simplest modifications of (2) is to add randomness to the update order. Randomness has been very useful in the analysis of block coordinate descent (BCD) methods and stochastic gradient descent (SGD) methods. In particular, a recent work Sun and Ye showed that randomized CD (R-CD) can be up to $O(n^{2})$ times faster than cyclic CD (C-CD) for quadratic minimization in the worst case, where $n$ is the number of variables Rigorously speaking, these two bounds are not directly comparable since the result for the randomized version only holds with high probability, while the result for the cyclic version always holds; anyhow, this $O(n^{2})$ gap is still meaningful if ignoring this difference between deterministic and randomized algorithm.. Another example is the comparison of IAG (Incremental Aggregated Gradient) in Blatt et al. and its randomized version SAG (Stochastic Average Gradient) : it turns out that the introduction of randomness leads to better iteration complexity bounds. There is also some study on randomly permuted version of pure SGD . These examples show that randomization may improve the algorithm in theory and in practice.

It is important to note that the iteration complexity bounds for randomized algorithms are usually established for independent randomization (sampling with replacement), while in practice, random permutation (sampling without replacement) has been reported to exhibit faster convergence (e.g. Shalev et al. , Recht and Re , Sun ). Interestingly, our simulation shows that for solving linear system of equations, randomly permuted ADMM (RP-ADMM) always converges, but independently randomized versions of ADMM can be divergent even for Gaussian data. Therefore, we focus on the analysis of RP-ADMM in this paper.

Random permutation is known to be notoriously difficult to analyze. Even for unconstrained quadratic minimization, the convergence rate of RP-BCD is poorly understood. Many existing works treated cyclic BCD and RP-BCD together , and thus the best known convergence rate of RP-BCD for general convex problems are in fact the same as that of C-BCD . However, in light of a recent study which established an up to $O(n^{2})$ gap between cyclic CD and R-CD , it is unlikely that RP-CD has the same convergence rate as C-CD since that would imply RP-CD could be $O(n^{2})$-times slower than R-CD. For the special example that demonstrates the gap between C-CD and R-CD, it was shown recently that RP-CD is faster than R-CD This paper appeared after the first version of the current paper. . However, the general quadratic case seems to be quite difficult, probably due to its close connection to a matrix AM-GM (algebraic mean-geometric mean) inequality , the difficulty of which is essentially to prove an inequality in non-commutative algebra.

We consider two extremes of a general RP-ADMM: i) the objective is zero, i.e., RP-ADMM for solving a linear system; ii) the constraint is zero and the objective is a quadratic function, i.e., RP-BCD for solving quadratic minimization. Due to the lack of understanding of random permutation for quadratic minimization as discussed previously, we restrict to the two cases in this paper.

The first result of this paper is the expected convergence of RP-ADMM for solving linear systems. More specifically, when the objective function is zero and the constraint is a non-singular square linear system of equations, the expected output of randomly permuted ADMM converges to the unique primal-dual optimal solution. A major technical result in this proof is that the eigenvalues of the expected iteration matrix of RP-BCD for quadratic problems lie in $(-1/3,1)$, instead of the typical range $(-1,1)$.

The second result is about the expected convergence rate of RP-ADMM for solving linear systems and RP-BCD for solving quadratic problems. We show that RP-BCD for a convex quadratic minimization problem with equal diagonal entries has expected iteration complexity $O(n\frac{\lambda_{\text{avg}}}{\lambda_{\min}}\log(1/\epsilon))$, where $\lambda_{\text{avg}}$ and $\lambda_{\min}$ are the average eigenvalue and the minimum eigenvalue of the coefficient matrix, and one “iteration” here means a cycle of updating all blocks. This improves an existing bound of $O(n^{2}\frac{\lambda_{\text{avg}}}{\lambda_{\min}}\log(1/\epsilon))$ for RP-BCD by a factor of $n$. Built on this result, we further show that RP-ADMM for solving linear systems achieves the same expected iteration complexity bound $O(n\frac{\lambda_{\text{avg}}}{\lambda_{\min}}\log(1/\epsilon))$.

Technically, we provide a simple and clean proof of the expected convergence, by applying a classical result on the eigenvalues of Jordan product. For proving the expected convergence rate, we propose a new variant of the matrix AM-GM inequality conjecture, and prove a weaker version of this conjecture.

Our result shows that random permutation may be a good answer to the question “how to apply the decomposition idea to solve constrained problems”. As multi-block BCD is widely used for large-scale unconstrained problems, we expect multi-block RP-ADMM to be a good candidate for large-scale linearly constrained problems. Our result provides one of the few direct analyzes of random permutation in optimization algorithms, and offers an explanation of the mysterious gap between RP-ADMM and cyclic ADMM. As reflected by the proof, the intuition is that random permutation provides “3-level symmetrization” that adjusts the spectrum of the update matrix. Based on the analysis for RP-ADMM, we are able to improve the best known complexity of RP-BCD for equally-diagonal quadratic problems by a factor of $n$, when expressing the complexity only in terms of the quantity $\frac{\lambda_{\text{avg}}}{\lambda_{\min}}$.

2 Related Works

This paper is a stronger version of a previous technical report Sun et al. which was not published. Another related work is the paper Chen et al. , which modifies the proof of to make it work with a quadratic objective function.

We highlight a few novel contributions of the current paper (neither in the original technical report nor in the paper ).

(i) The current paper provides a much simpler proof for the result of expected convergence.

(ii) The current paper provides the first convergence rate analysis of RP-ADMM. See Theorem 4 and the proof in Section 4.5, Section 7.2 and Section 7.1.

(iii) The current paper provides an improved convergence rate analysis of RP-BCD, See Theorem 3 and the proof in Section 4.3 and Section 7.3.

(iv) The current paper introduces a theory-motivated algorithm Bernoulli-ADMM, which reduces the sampling time yet still achieves the expected convergence. This update order has not appeared before even in other algorithm setups to our knowledge. See Section 2.5 and Proposition 1.

Besides the technical contributions, we want to emphasize that the current paper is not just adding new result to our previous technical report , but actually completes a missing step of the story. From a mathematical point of view, the most striking consequence of our original proof is that the spectral radius of RP-BCD lies in a smaller region $(-1/3,1)$. It is natural to think that this fundamental fact should have an impact on the analysis of original RP-BCD. Our current paper fills this gap by showing that this result can help build an $O(n)$ gap between the (expected) onvergence rate of RP-BCD and cyclic BCD. A general message is that on one hand, to understand constrained optimization we have to understand unconstrained optimization (analyzing ADMM reduces to analyzing BCD); on the other hand, analyzing constrained optimization helps improve the understanding of unconstrained optimization (the analysis of ADMM leads to progress in BCD). We find this interaction between unconstrained optimization (BCD) and constrained optimization (ADMM) fascinating. The whole story is only revealed in the current paper, but not in the previous technical report or Chen et al. .

Besides the above unique aspects, the current paper inherits some interesting numerical findings from the technical report Sun et al. which do not appear in Chen et al. . We find that cyclic ADMM diverges with probability 1 for many random distributions of data, thus showing that the seemingly surprising divergence behavior reported in is quite common. However, it is easy to miss this finding if one uses the Gaussian distribution to generate data. Another interesting finding is that the independently randomized version of ADMM diverges with probability 1 for Gaussian data but not for the counter-example in , preventing us from analyzing the independently randomized version. Without these findings, the motivation of studying RP-ADMM would be less clear. See Section 2.4 and Section 8.

3 Notation and Organization

Organization. In Section 2, we present three versions of randomized ADMM, with an emphasis on RP-ADMM. In Section 3, we present our main results Theorem 1, Theorem 2 and their proofs. The subsequent sections are devoted to the proofs of the two technical results Lemma 1 and Lemma 2, which are used in the proof of Theorem 2. In particular, the proof of Lemma 1 is given in Section 5, and the proof of Lemma 2 is given in Section 6.

Algorithms

In this section, we will present both randomly permuted and independently randomized versions of ADMM for solving (1), and specialize RP-ADMM for solving a square system of equations. We also present a rather novel algorithm Bernoulli-randomized ADMM (motivated by our proof).

In this subsection, we first propose RP-ADMM for solving the general optimization problem (1), then we present the update equation of RP-ADMM for solving a linear system of equations.

At each round, we draw a permutation $\sigma$ of $\{1,\dots,n\}$ uniformly at random from $\Gamma$, and update the primal variables in the order of the permutation, followed by updating the dual variables in a usual way. Obviously, all primal and dual variables are updated exactly once at each round. See Algorithm 1 for the details of RP-ADMM. Note that with a little abuse of notation, the function $\mathcal{L}(x_{\sigma(1)},x_{\sigma(2)},\dots,x_{\sigma(n)};\mu)$ in this algorithm should be understood as $\mathcal{L}(x_{1},x_{2},\dots,x_{n};\mu)$. For example, when $n=3$ and $\sigma=(231)$, $\mathcal{L}(x_{\sigma(1)},x_{\sigma(2)},x_{\sigma(3)};\mu)=\mathcal{L}(x_{2},x_{3},x_{1};\mu)$ should be understood as $\mathcal{L}(x_{1},x_{2},x_{3};\mu)$.

Throughout this paper, we assume $A$ is non-singular. Then the unique solution to (8) is $x=A^{-1}b$, and problem (7) has a unique primal-dual optimal solution $(x,\mu)=(A^{-1}b,0)$. The augmented Lagrangian function (3) for the optimization problem (7) becomes

Throughout this paper, we assume $\beta=1$; note that our algorithms and results can be extended to any $\beta>0$ by simply scaling $\mu$.

1.2 Example of 333-block ADMM

then the update equation (10) becomes $\bar{L}y^{k+1}=\bar{R}y^{k}+\bar{b}$, i.e.

1.3 General Update Equation of RP-ADMM

In general, for the optimization problem (7), the primal update (5) becomes

Replacing $\sigma(i),\sigma(j),\sigma(l)$ by $i,j,l$, we can rewrite the above equation as

where $\sigma^{-1}$ denotes the inverse mapping of a permutation $\sigma$, i.e. $\sigma(i)=t\Leftrightarrow i=\sigma^{-1}(t)$. Denote the output of Algorithm 1 after round $(k-1)$ as

The update equations of Algorithm 1 for solving (7), i.e. (15) and (6), can be written in the matrix form as (when the permutation is $\sigma$ and $\beta=1$)

where $\bar{L}_{\sigma},\bar{R}_{\sigma},L_{\sigma},R_{\sigma},\bar{b}$ are defined by

Another expression of $L_{\sigma}$, equivalent to (19), is the following:

A user-friendly rule for writing $L_{\sigma}$ is described as follows (use $\sigma=(231)$ as an example). Start from a zero matrix. First, find all reverse pairs of $\sigma$; here, we say $(i,j)$ is a reverse pair if $i$ appears after $j$ in $\sigma$. For the permutation $(231)$, all the reverse pairs are $(1,3),(3,2)$ and $(1,2)$. Second, in the positions corresponding to the reverse pairs, write down the corresponding entries of $A^{T}A$, i.e. $a_{1}^{T}a_{3},a_{3}^{T}a_{2}$ and $a_{1}^{T}a_{2}$, respectively. At last, write $a_{i}^{T}a_{i}$ in the diagonal positions. Using this rule, we can write down the expression of $L_{(231)}$ as

A user-friendly rule to quickly check the correctness of an expression of $L_{\sigma}$ is the following (still take $\sigma=(231)$ as an example). According to the order of the permutation $(231)$, the $2$nd row, the $3$rd row and the $1$st row should have a strictly decreasing number of zeros ($2$ zeros, $1$ zero and no zero). In contrast, the $2$nd column, the $3$rd column and the $1$st column should have a strictly increasing number of zeros.

For the general case that $d_{i}\geq 1,\forall i$, we can write down the block partitioned $L_{\sigma}$ in a similar way. For example, when $n=3$ and $\sigma=(231)$, we have

2 Randomly Permuted BCD

RP-ADMM is a generalization of RP-BCD. In fact, when the constraint does not exist, RP-ADMM reduces to RP-BCD. In this subsection, we present RP-BCD for solving convex quadratic problems. Note that RP-ADMM for solving linear systems and RP-BCD fo solving quadratic problems are two extremes of general RP-ADMM: in the former case the objective function is zero, and in the latter case the constraint is zero. Interestingly, the two extreme cases are related as the expected iteration matrix of RP-BCD appears as a component of the expected iteration matrix of RP-ADMM. We will show later that their eigenvalues are closely related.

In the augmented Lagrangian function given in (9), if we delete the first term which depends on the dual variable $\mu$, we obtain the quadratic function $\frac{1}{2}\|Ax-b\|^{2}$. Thus if we eliminate the dual variable $\mu$ in the update equations of RP-ADMM, we will obtain the update equations for RP-BCD. Suppose $x^{k}$ is the iterate after the k-th epoch (i.e. go through all coordinates once), and $\sigma$ is the order used in the $k$-th iteration, then, as a simpler version of (17), we have

where $L_{\sigma}$ and $R_{\sigma}$ are defined as in (21) and (20), and $\sigma$ is a random permutation.

3 Residual Trick for Efficient Implementation of ADMM and BCD

We note here that when $d_{i}=1,\forall i$ (in this case BCD becomes CD), per-epoch computation time of ADMM and CD (no matter what order) is $O(n^{2})$; or in other words, per-coordinate-update time is $O(n)$. For instance, updating $x^{k+1}$ by (24) in RP-BCD or updating $y^{k+1}$ by (17) in RP-ADMM only takes time $O(n^{2})$. As mentioned in Section 3.1 of , the trick is to keep track of the residual. For both efficient practical implementation and calculation of computation complexity, one should use this residual trick, but for the ease of theoretical analysis we use the matrix update forms (17) and (24) in this paper; there is no contradiction as our theory only depends on the value of $x^{k}$ but not the specific procedure to compute $x^{k}$.

For completeness, we briefly explain how this trick works in our settings. Suppose $d_{i}=1,\forall i$, and we use CD methods to solve (23) with a certain update order (could be any order, such as cyclic, randomized or randomly permuted). Suppose the coordinate $i$ is picked, then $x_{i}$ is updated by by

where $A_{-i}$ contains all columns of $A$ except $A_{i}$, $x_{-i}$ contains all elements of $x$ except $x_{i}$ and represents the current values, and $x_{i}^{+}$ represents the new value. A straightforward implementation of (25) requires multiplying $x_{-i}$ by $A_{-i}$ which takes $O(n^{2})$ operations. With the residual trick (e.g. ), we introduce the residual $r=Ax-b$, and replace (25) by

Now the calculation of $x_{i}^{T}$ and $r^{+}$ takes time $O(n)$, and thus one epoch of BCD takes time $O(n^{2})$. The same trick can be applied to the primal update of ADMM; with this trick, the dual update (6) can be rewritten as $\mu^{+}=\mu-\beta r$ which takes time $O(n)$, and thus one epoch of ADMM takes time $O(n^{2})$.

Finally, when $d_{i}>1$, similar update equations can still be used except a minor difference that $\frac{1}{A_{i}^{T}A_{i}}$ should be replaced by $(A_{i}^{T}A_{i})^{-1}$. In a special case that $d_{i}=d,\forall i$ and $N=dn$, each iteration of BCD takes time $O(Nd+d^{3})$ and each epoch takes time $O(N^{2}+Nd^{2})$. This cost can be reduced if we use BCGD (i.e. not solving the subproblem exactly but updating each block of variables by a gradient step). In order not to make the paper more complicated, we will not discuss the inexact versions of BCD and ADMM in this paper.

4 Two Versions of Independently Randomized ADMM

In this subsection, we present two other versions of randomized ADMM which can be divergent according to simulations. The failure of these versions makes us focus on analyzing RP-ADMM in this paper. These versions can be viewed as natural extensions of R-BCD (randomized BCD) and .

In the first algorithm, called primal-dual randomized ADMM (PD-RADMM), the whole dual variable is viewed as the $(n+1)$-th block. In particular, at each iteration, the algorithm draws one index $i$ from $\{1,\dots,n,n+1\}$, then performs the following update: if $i\leq n$, update the $i$-th block of the primal variable; if $i=n+1$, update the whole dual variable. The details are given in Algorithm 2. We have tested PD-RADMM for the counter-example given in Chen et al. , and found that PD-RADMM always diverges (for random initial points).

A variant of PD-RADMM has been proposed in Hong et al. with two differences: first, instead of minimizing the augmented Lagrangian $\mathcal{L}$, that algorithm minimizes a strongly convex upper bound of $\mathcal{L}$; second, that algorithm uses a diminishing dual stepsize. With these two modifications, shows that each limit point of the sequence generated by their algorithm is a primal-dual optimum with probability 1. Note that also proves the same convergence result for the cyclic version of multi-block ADMM with these two modifications, thus it does not show the benefit of randomization.

In the second algorithm, called primal randomized ADMM (P-RADMM), we only perform randomization for the primal variables. In particular, at each round, we first draw $n$ independent random variables $j_{1},\dots,j_{n}$ from the uniform distribution of $\{1,\dots,n\}$ and update $x_{j_{1}},\dots,x_{j_{n}}$ sequentially, then update the dual variable in the usual way. The details are given in Algorithm 3. This algorithm looks quite similar to RP-ADMM as they both update $n$ primal blocks at each round; the difference is that RP-ADMM samples without replacement while this algorithm P-RADMM samples with replacement. In other words, RP-ADMM updates each block exactly once at each round, while P-RADMM may update one block more than one times or does not update one block at each round.

We have tested P-RADMM in various settings. For the counter-example given in Chen et al. , we found that P-RADMM does converge. However, if $n\geq 30$ and $A$ is a Gaussian random matrix (each entry is drawn i.i.d. from $\mathcal{N}(0,1)$), then P-RADMM diverges in almost all cases we have tested. This phenomenon is rather strange since for random Gaussian matrices $A$ the cyclic ADMM actually converges (according to simulations). An implication is that randomized versions do not always outperform their deterministic counterparts in terms of convergence.

Since both Algorithm 2 and Algorithm 3 can diverge in certain cases, we will not further study them in this paper. In the rest of the paper, we will focus on RP-ADMM (i.e. Algorithm 1).

5 Bernoulli-Randomized ADMM

To implement randomly permuted ADMM, one needs to sample from all blocks without replacement. To save the sampling time, we propose another algorithm which we call Bernoulli-randomized ADMM. This algorithm is motivated by the proof of Theorem 1. This updating scheme can be applied to other algorithms such as SGD and coordinate descent methods.

The new update order combines the well-known double-sweep order and Bernoulli-randomization. The original double-sweep order is $(1,2,...,n-1,n,n-1,n-2,...,1)$, meaning that $x_{1},x_{2},\dots,x_{n-1},x_{n},x_{n-1},x_{n-2},\dots,x_{1}$ are updated sequentially in each “cycle”. It combines the normal cyclic order $(1,2,\dots,n)$ and a reverse order $(n,n-1,\dots,1)$. We propose the following updating scheme: add a check box to each block, and in each cycle we perform the following operations.

Phase I: go through the blocks $x_{1},x_{2},\dots,x_{n}$ one by one sequentially as follows: for each block $x_{i}$, flip a fair coin and:

if the outcome is “head”, update the block $x_{i}$ and check the check box;

if the outcome is “tail”, do nothing about $x_{i}$ and uncheck the check box.

Phase II: go through the blocks $x_{n},x_{n_{1}},\dots,x_{1}$ in the reverse order, and update $x_{i}$ if the box is unchecked.

Note that in each cycle we go through each block twice but update each block exactly once so that the number of totally updated blocks remains $n$. For example, when $n=5$, $(35421)$ is a possible update order, as shown in the following diagram.

Similarly, $(13542)$ is also a possible update order. But $(13524)$ and $(35412)$ are not possible. The set of all possible update orders is given by

where $\Gamma$ is the set of permutations of $\{1,2,\dots,n\}$ as defined in (4). In other words, a sequence from $\Gamma_{\text{BR}}$ is a concatenation of an increasing sequence and a decreasing sequence. Note that the permutation $(1,2,...,n)$ is in $\Gamma_{\text{BR}}$ since it can be viewed as the concatenation of an increasing sequence $(1,2,...,n-1)$ and a “decreasing sequence” $(n)$, and we can let $i=n-1$ in the above definition to cover this case. Similarly, the permutation $(n,n-1,\dots,1)$ is also in $\Gamma_{\text{BR}}$ as $i=1$ will cover this case.

The algorithm Bernoulli-randomized ADMM (BR-ADMM) is formally described below. We skip the epoch index $k$ since otherwise the notation would be cumbersome.

For solving linear systems of equations, the update formula is the same as (17), the update formula of RP-ADMM. The difference is that for RP-ADMM $\sigma$ can be an arbitrary permuation, while for BR-ADMM there is some restriction on $\sigma$: it has to be a permuation in $\Gamma_{\text{BR}}$.

Main Results

Let $\sigma_{i}$ denote the permutation used in round $i$ of Algorithm 1, which is a uniform random variable drawn from the set of permutations $\Gamma$. After round $k$, Algorithm 1 generates a random output $y^{k+1}$, which depends on the observed draw of the random variable

We will show that the expected iterate (the iterate $y^{k}$ is defined in (16))

Assume the coefficient matrix $A=[A_{1},\dots,A_{n}]$ of the constraint in (7) is a non-singular square matrix. Suppose Algorithm 1 is used to solve problem (7), then the expected output converges to the unique primal-dual optimal solution to (7), i.e.

Since the update matrix does not depend on previous iterates, we claim (and prove in Section 4.1) that Theorem 1 holds if the expected update matrix has a spectral radius less than 1, i.e. if the following Theorem 2 holds.

where the expectation is taken over the uniform random distribution over $\Gamma$, the set of permutations of $\{1,2,\dots,n\}$. Then the spectral radius of $M$ is smaller than $1$, i.e.

For the counterexample in Chen et al. where $A=[1,1,1;1,1,2;1,2,2]$, it is easy to verify that $\rho(M_{\sigma})>1.02$ for any permutation $\sigma$ of $(1,2,3)$. Interestingly, Theorem 2 shows that even if each $M_{\sigma}$ is “bad” (with spectral radius larger than $1$), the average of them is always “good” (with spectral radius smaller than $1$).

Theorem 2 is just a linear algebra result, and can be understood even without knowing the details of the algorithm. However, the proof of Theorem 2 is rather non-trivial. This proof will be provided in Section 4.2, and the technical results used in this proof will be proved in Section 5 and Section 6.

The convergence rate of RP-ADMM for solving linear systems of equations is closely related to the convergence rate of RP-BCD (randomly permuted BCD) for solving quadratic problems. We will discuss their relation and how our results in this paper improve our understanding for RP-BCD.

A similar convergence result holds for BR-ADMM proposed in Section 2.5, as presented below. The proof is a simple modification of the proof of Theorem 1, and can be found in Section 6.4.

Assume the coefficient matrix $A=[A_{1},\dots,A_{n}]$ of the constraint in (7) is a non-singular square matrix. Suppose Algorithm 4 is used to solve problem (7), then the expected output converges to the unique primal-dual optimal solution to (7).

2 Expected Convergence Rate of RP-ADMM and RP-BCD

There is a close relation between RP-ADMM for solving linear systems and RP-CD for solving quadratic problems (see Lemma 2). Thus it is not surprising that we need to understand RP-BCD before understanding RP-ADMM. We will first present an expected convergence rate of RP-BCD (in terms of the expected iterates) for solving quadratic problems, which improves the best existing convergence rate (one type of rates, to be precise) by a factor of $n$ Rigorously speaking, this is not a fair comparison as the complexity of C-CD is deterministic complexity.. The result is proved via establishing a weak version of matrix AM-GM inequality. This result also establishes a large gap of $O(n)$ between RP-BCD and C-BCD (cyclic BCD). Second, built upon the result for RP-BCD, we establish a convergence rate of RP-ADMM which is similar to RP-BCD and also $n$ times better than that of C-BCD.

The first result is about the expected convergence rate of RP-BCD for the case $A_{i}^{T}A_{i}=I$. This assumption is made so that the expression is simple, and the case for general $A_{i}$ is given in the next result.

(rate of RP-BCD for quadratic functions with identity diagonal blocks) Assume the coefficient matrix $A=[A_{1},\dots,A_{n}]$ is a non-singular square matrix, and $A_{i}^{T}A_{i}=I,\forall\;i.$ Suppose RP-BCD is used to solve problem (23), where $x^{k}$ denotes the variable after $k$ epochs (each epoch represents one cycle of updating all coordinates). Denote the unique optimal solution as $x^{*}=A^{-1}b$. Then

To put this convergence rate result in the context, we consider the simple case that each $d_{i}=1$, i.e., each block consists of a single coordinate. In this case, every diagonal entry of $A^{T}A$ is $1$, thus the average eigenvalue of $A^{T}A$ is $1$. Throughout the paper, we consider the total computation complexity The computation complexity equals the iteration complexity times the per-iteration cost. We do not present iteration complexity since there may be confusion about whether “one iteration” means $n$ coordinate updates or $1$ coordinate update. Presenting iteration complexity is better if one considers a general convex problem, but then one needs to discuss the per-iteration cost. We are considering quadratic problems throughout the paper, so we feel it is more clear to stick to computation complexity.; note that we assume the residual trick as described in 2.3 is always used for all methods.

Our Theorem 3 provides an expected computational complexity upper bound $O(n^{3}\kappa_{\text{CD}}\log\frac{1}{\epsilon})$ for RP-CD, since each epoch takes $O(n^{2})$ time and it requires $O(n^{2}\log\frac{1}{\epsilon})$ epochs to achieve error $\epsilon$ according to (33). It is known that the computational complexity of R-CD (randomized coordinate descent) to achieve relative accuracy $\epsilon$ Here, the relative accuracy $\epsilon$ means $\|E(x^{k})-x^{*}\|/\|x^{0}-x^{*}\|$ or $\|E(f(x^{k}))-f^{*}\|/\|f(x^{0})-f^{*}\|$. is $O(n^{2}\kappa_{\text{CD}}\log\frac{1}{\epsilon})$, where $\kappa_{\text{CD}}=\lambda_{\text{avg}}(A^{T}A)/\lambda_{\min}(A^{T}A)=1/\lambda_{\min}(A^{T}A)$ is the ratio of the average eigenvalue over the minimum eigenvalue. It was recently shown that in terms of $\kappa_{\text{CD}}$ and $n$ only, the worst-case complexity of C-CD (cyclic CD) is $O(n^{4}\kappa_{\text{CD}}\log\frac{1}{\epsilon})$, which is $n^{2}$ times worse than R-CD and $n$ times worse than GD. This shows a large gap between C-CD and R-CD in the worst case.

It was widely conjectured that RP-CD is at least as fast as R-CD, but this conjecture is considered to be rather difficult to prove. For a special class of matrices, recent works validated the conjecture. However, to our knowledge, even for a general quadratic function with equal diagonal entries $1$, the previously best known convergence rate of RP-CD is almost the same as C-CD (see ), which can be $n^{2}$ times worse than that of R-CD. Our Theorem 3 provides an expected computational complexity upper bound $O(n^{3}\kappa_{\text{CD}}\log\frac{1}{\epsilon})$ for RP-CD, which is $n$ times faster than C-CD and $n$ times slower than R-CD. This improves the best existing rate by a factor of $n$ Note that this “improvement” is valid when the convergence rate is characterized by only $\kappa_{\text{CD}}$ and $n$. It is common to use other parameters such as the maximum eigenvalue to characterize the convergence rate (see for a detailed discussion), and our result here does not provide improvement for other kinds of convergence rate.. We summarize the comparison of the complexity for C-CD, R-CD and RP-CD in Table 1.

The following proposition generalizes Theorem 3 to the non-identity-diagonal case, i.e., $A_{i}^{T}A_{i}$ does not need to be an identity matrix.

(rate of RP-BCD for quadratic functions, with non-identity blocks) Assume the coefficient matrix $A=[A_{1},\dots,A_{n}]$ is a non-singular square matrix. Suppose RP-BCD is used to solve problem (23). Denote $D=\text{diag}(A_{1}^{T}A_{1},\dots,A_{n}^{T}A_{n})$ as a block-diagonal matrix, and the norm $\|z\|_{D}=\sqrt{z^{T}Dz}$. Then

The proof of Proposition 2 is given in Section 4.4. One can easily transform the quantity $\lambda_{\min}(D^{1/2}A^{T}AD^{-1/2})$ to certain quantity that only depends on the eigenvalues of $A_{i}^{T}A_{i}$ and $A^{T}A$. However, as noted in , it is far from clear how tight the transformation is, thus we skip the transformation here. In fact, it is related to some open question on the so-called Jacobi-preconditioning. We refer the interested readers to for a detailed discussion of the subtle issues in the non-identity-diagonal case.

At last, we present a result on the expected convergence rate of RP-ADMM for solving linear systems, under the assumption that $A_{i}^{T}A_{i}=I,\;\forall i$. Very similar to Proposition 2, we can also generalize this result to non-identity-diagonal case, i.e., $A_{i}^{T}A_{i}\neq I$, but to save space we skip the generalization here. The proof of Theorem 4 is given in Section 4.5.

(Expected convergence rate of RP-ADMM for linear systems) Assume the coefficient matrix $A=[A_{1},\dots,A_{n}]$ of the constraint in (7) is a non-singular square matrix and $A_{i}^{T}A_{i}=I_{d_{i}}$. Suppose Algorithm 1 is used to solve problem (7). Denote $y^{*}=\begin{bmatrix}A^{-1}b\\ 0\\ \end{bmatrix}$ as the unique primal-dual optimal solution to the problem (7), then

This result implies that similar to RP-CD for solving quadratic problems, the complexity of RP-ADMM in terms of the expected iterates for solving linear systems is also at most

In light of the fact that C-CD has been shown to only achieve a rate $O(n^{4}\kappa_{\text{CD}}\log(1/\epsilon))$ , the rate of RP-ADMM we obtain is already quite good. Nevertheless, we conjecture that this complexity upper bound can be improved to $O(n^{2}\kappa_{\text{CD}}\log(1/\epsilon))$, the same as the conjectured complexity for RP-CD. But an improved rate of RP-ADMM leads to an improved rate of RP-BCD (this should be clear via the comparison of (50) and (57)), thus proving this conjecture is an even more difficult problem than the long-standing open question on RP-CD.

3 Matrix AM-GM Inequality

The original version is more general: the number of matrices does not need to be the same as the dimension of the matrix. For simplicity, we just present a simpler version here.

The matrix AM-GM inequality is a generalization of the well-known AM-GM inequality: for non-negative numbers $a_{1},\dots,a_{n}$, the geometric mean $(a_{1}a_{2}\dots a_{n})^{1/n}$ is no more than the algebraic mean $\frac{1}{n}\sum_{i=1}^{n}a_{i}$. When extending this inequality to matrix domain, the non-commutative nature of matrix multiplication makes the problem rather difficult to prove.

We observe that we only need to prove a matrix AM-GM inequality for projection matrices. We conjecture that the following matrix AM-GM inequality holds.

Compared with (36), our conjecture makes a stronger claim on the relation, but it only applies to projection matrices. We have found examples to show that (37) does not hold for general positive semi-definite matrices, but it holds for projection matrices in all of our experiments.

We are not able to prove the new conjecture – that would solve the open question of the best convergence rate of RP-CD for quadratic problem. Nevertheless, inspired by the new conjecture, we prove a weaker version (see Lemma 3), which can lead to an improved convergence rate estimate for RP-CD.

Proof of Main Results

Denote $\sigma_{k}$ as the permutation used in round $k$, and define $\xi_{k}$ as in (28). Rewrite the update equation (17) below (replacing $\sigma$ by $\sigma_{k}$):

We first prove (30) for the case $b=0$. By (18) we have $\bar{b}=0$, then (38) is simplified to $y^{k+1}=\bar{L}_{\sigma_{k}}^{-1}\bar{R}_{\sigma_{k}}y^{k}$. Taking the expectation of both sides of this equation in $\xi_{k}$ (see its definition in (28)), and note that $y^{k}$ is independent of $\sigma_{k}$, we get

Since the spectral radius of $M$ is less than 1 by Theorem 2, we have that $\{\phi^{k}\}\rightarrow 0$, i.e. (30).

We then prove (30) for general $b$. Let $y^{*}=[A^{-1}b;0]$ denote the optimal solution. Then it is easy to verify that

for all $\sigma_{k}\in\Gamma$ (i.e. the optimal solution is the fixed point of the update equation for any order). Compute the difference between this equation and (38) and letting $\hat{y}^{k}=y^{k}-y^{*}$ , we get $\hat{y}^{k+1}=\bar{L}_{\sigma_{k}}^{-1}\bar{R}_{\sigma_{k}}\hat{y}^{k}$. According to the proof for the case $b=0$, we have $E(\hat{y}^{k})\longrightarrow 0$, which implies $E(y^{k})\longrightarrow y^{*}$.

2 Proof of Theorem 2

The difficulty of proving Theorem 2 (bounding the spectral radius of $M$ defined in (31)) is two-fold. First, $M$ is a non-symmetric matrix, and there are very few tools to bound the spectral radius of a non-symmetric matrix. In fact, spectral radius is neither subadditive nor submultiplicative (see, e.g. Kittaneh ). Note that the spectral norm of $M$ can be much larger than $1$ (there are examples that $\|M\|>2$), thus we cannot bound the spectral radius simply by the spectral norm. Second, although it is possible to explicitly write each entry of $M$ as a function of the entries of $A^{T}A$, these functions are very complicated ($n$-th order polynomials) and it is not clear how to utilize this explicit expression.

The proof outline of Theorem 2 and the main techniques are described below. In Step 0, we provide an expression of the expected update matrix $M$. In Step 1, we establish the relationship between the eigenvalues of $M$ and the eigenvalues of a simple symmetric matrix $AQA^{T}$, where $Q$ is defined in (39). As a consequence, the spectral radius of $M$ is smaller than one iff the eigenvalues of $AQA^{T}$ lie in the region $(0,4/3)$. This step partially resolves the first difficulty, i.e. how to deal with the spectral radius of a non-symmetric matrix. In Step 2, we show that the eigenvalues of $AQA^{T}$ do lie in $(0,4/3)$ using mathematical induction. The induction analysis circumvents the second difficulty, i.e. how to utilize the relation between $M$ and $A$.

Step 0: compute the expression of the expected update matrix $M$. Define

It is easy to prove that $Q$ defined by (39) is symmetric. In fact, note that $L_{\sigma}^{T}=L_{\bar{\sigma}},\forall\sigma\in\Gamma$, where $\bar{\sigma}$ is a reverse permutation of $\sigma$ satisfying $\bar{\sigma}(i)=\sigma(n+1-i),\forall\ i$, thus $Q=\frac{1}{n!}\sum_{\sigma}Q_{\sigma}=(\frac{1}{n!}\sum_{\sigma}Q_{\bar{\sigma}})^{T}=Q^{T},$ where the last step is because the sum of all $Q_{\bar{\sigma}}$ is the same as the sum of all $Q_{\sigma}$.

Substituting the expression of $\bar{L}_{\sigma}^{-1}$ into the above relation, and replacing $R_{\sigma}$ by $L_{\sigma}-A^{T}A$, we obtain

Since $M_{\sigma}$ is linear in $L_{\sigma}^{-1}$, we have

Step 1: relate $M$ to a simple symmetric matrix. The main result of Step 1 is given below, and the proof of this result is relegated to Section 5.

Furthermore, when $Q$ is symmetric, we have

Remark: For our problem, the matrix $Q$ as defined by (39) is symmetric (see the argument after equation (39)), thus the relation (45) indeed holds according to Lemma 1. For a general non-symmetric $Q$, (45) does not need to hold, but the first conclusion (44) still holds.

Step 2: Bound the eigenvalues of $QA^{T}A$. The main result of Step 2 is summarized in the following Lemma 2. The proof of Lemma 2 is given in Section 6.

in which $L_{\sigma}$ is defined by (21) and $\Gamma$ is defined by (4). Then all eigenvalues of $QA^{T}A$ lie in $(0,4/3)$, i.e.

Remark: The upper bound $\frac{4}{3}$ in (47) is probably tight, since we have found numerical examples with $\text{eig}(QA^{T}A)>1.3333$. Now the expected convergence of RP-ADMM seems to be a pleasant coincidence: Lemma 1 shows that to prove the expected convergence we need to prove $\sup_{A}\text{eig}(QA^{T}A)$, a quantity that can be defined without knowing ADMM, is bounded by $4/3$; Lemma 2 and numerical experiments show that this quantity happens to be exactly $4/3$ so that RP-ADMM can converge (in expectation).

Theorem 2 follows immediately from Lemma 1 and Lemma 2.

3 Proof of Theorem 3

We first describe the outline of the proof. The expected update matrix of RP-BCD is $I-QA^{T}A$, and the eigenvalues of this matrix lie in $(-1,1)$. The expected convergence speed of RP-BCD depends on the distance between the eigenvalues and the two extremes $-1$ and $1$. Lemma 2 shows that the distance to $-1$ is at least $1/3$, which is a constant. We will show that the distance to $1$ is at least $\lambda_{\min}(A^{T}A)/n$, by proving a weaker version of matrix AM-GM inequality. Combining the two results, we obtain the expected convergence speed of RP-BCD.

According to (24), we have $x^{k+1}-x^{*}=(I-L_{\sigma}^{-1}A^{T}A)(x^{k}-x^{*})$, where $\sigma$ is the randomly picked permutation at the $k$-th epoch. Therefore, the expected update formula of RP-BCD for solving the least squares problem is

Suppose the eigenvalues of $QA^{T}A$ are $\eta_{1}\geq\eta_{2}\geq\dots\geq\eta_{n}$, then according to Lemma 2,

thus the spectral radius of $I-QA^{T}A$ is

An interesting phenomenon occurs here. The spectral radius is either $1-\eta_{n}$ or $|1-\eta_{1}|$. In the latter case, $\rho(I-QA^{T}A)=|1-\eta_{1}|\leq 1/3$, implying that $\|E(x^{k})-x^{*}\|\leq\frac{1}{3^{k}}\|E(x^{0})-x^{*}\|$, or equivalently, the relative error $|E(x^{k})-x^{*}\|/|E(x^{0})-x^{*}\|$ achieves $\epsilon$ in $\log 3\log(1/\epsilon)$ epochs. We do not even need to compute $\eta_{1}$ since it will only affect the convergence speed when the speed is already very fast. From a theoretical perspective, the improvement from $\log 3$ to $\log(1/(1-|1-\eta_{1}|))$ is just an improvment in the constant. Therefore, it is reasonable to ignore $\eta_{1}$ and focus on the estimate of $1-\eta_{n}$.

To estimate the maximum eigenvalue of $I-QA^{T}A$ (or equivalently, that of $I-AQA^{T}$), we first provide a useful identity that connects $I-AQA^{T}$ and projection matrices $P_{i}=I-A_{i}A_{i}^{T}$.

Suppose $A=[A_{1},\dots,A_{n}]$ is a non-singular square matrix, and $A_{i}^{T}A_{i}=I,\forall\;i.$ For a permutation $\sigma=(\sigma_{1},\dots,\sigma_{n})\in\Gamma$, $L_{\sigma}$ is defined as in (21), and $Q_{\sigma}=L_{\sigma}^{-1}$. Denote $P_{i}=I-A_{i}A_{i}^{T}$, $i=1,\dots,n$. Then we have

The proof of Claim 4.1 is given at the end of this subsection. Claim 4.1 states that $I-AQA^{T}$ is exactly equal to $\frac{1}{n!}\sum_{\sigma=(\sigma_{1},\dots,\sigma_{n})\in\Gamma}P_{\sigma_{n}}P_{\sigma_{n-1}}\dots P_{\sigma_{1}}$, thus we only need to estimate the maximal eigenvalue of the latter expression. This is achieved by the following lemma (the proof is given in Section 7.3).

The above Lemma 3 and Claim 4.1 immediately lead to the following corollary.

Suppose $A=[A_{1},\dots,A_{n}]$ is a non-singular square matrix, and $A_{i}^{T}A_{i}=I,\forall\;i.$ Suppose $P_{i}=I-A_{i}A_{i}^{T},\;\forall\;i.$ $L_{\sigma}$ is defined as in (21), and $Q=E_{\sigma}(L_{\sigma}^{-1})$. Then

Note that $\frac{1}{n}\sum_{i}P_{i}=\frac{1}{n}(nI-\sum_{i}A_{i}A_{i}^{T})=I-\frac{1}{n}AA^{T}$, thus (53) implies

Substituting this relation into (49), we obatain

Remark: There is a coefficient $1/n$ in front of $\lambda_{\min}(AA^{T})$ in (54), and this is why the complexity of RP-CD we establish is $n$ times worse than the conjectured one in Table 1. If Conjecture 3.1 holds, then this factor of $1/n$ would be removed and the conjectured (expected) complexity of RP-CD in Table 1 would hold.

We prove (51a) by induction on $n$. Without loss of generality, we can assume $\sigma=(1,2,\dots,n)$, then $L_{\sigma}=\begin{bmatrix}A_{1}^{T}A_{1}&0&\dots&0\\ A_{2}^{T}A_{1}&A_{2}^{T}A_{2}&\dots&0\\ \vdots&\vdots&\ddots&\vdots\\ A_{n}^{T}A_{1}&A_{n}^{T}A_{2}&\dots&A_{n}^{T}A_{n}\\ \end{bmatrix}.$ In this case, (51a) becomes

The expression obviously holds for $n=1$. Suppose the expression holds for $n-1$, i.e., for $\hat{A}=[A_{1},\dots,A_{n-1}]$, we have

where $\hat{\sigma}=(1,2,\dots,n-1)$ is a permutation of $n-1$ elements and $\hat{L}_{\hat{\sigma}}$ is the counterpart of $L_{\sigma}$ for $n-1$ blocks defined as

The two matrices $L_{\sigma}$ and $\hat{L}_{\sigma^{\prime}}$ are related by

where in the last step we use the induction hypothesis (55). Thus we have proved (51a). Summing up (51a) for all possible permutations $\sigma$ and divide by $n!$, we obtain (51b). $\Box$

4 Proof of Proposition 2

According to (48), the (expected) update equation of RP-BCD is given by $E(x^{k+1})-x^{*}=(I-QA^{T}A)(E(x^{k})-x^{*})=Z(E(x^{k})-x^{*})$, where $Z=I-QA^{T}A=I-E(L_{\sigma}^{-1}A^{T}A)$.

5 Proof of Theorem 4

Now we consider the expected convergence rate of RP-ADMM. The difference with the analysis for RP-BCD is that here we need to consider the distance between the eigenvalues of $I-AQA^{T}$ with $-1/3$ while for RP-BCD what matters is the distance between the eigenvalues of $I-AQA^{T}$ and $-1$ which is at least $2/3$ and thus can be ignored.

Suppose the minimum and maximum eigenvalues of $QA^{T}A$ are $0<\tau_{\min}\leq\tau_{\max}<4/3$. Then

where $z_{+}=max\{z,0\}$. Furthermore, we have

The proof of Claim 4.2 is given in Section 7.1. The next lemma provides a universal estimate of the maximum eigenvalules of $QA^{T}A$.

The maximum eigenvalues of $QA^{T}A$ is at most $\frac{4}{3}-\frac{4}{9}\frac{1}{n+1}$, i.e.,

The proof of Lemma 4 is given in Section 7.2

According to (54), which is established in the proof of the expected convergence rate of RP-BCD, we have

Substituting the bounds (58) and (59) into (57), we obtain

Since $\lambda_{\min}(A^{T}A)\leq 1$, $\frac{1}{2n}\leq\frac{1}{n+1}$, this bound can be simplified to

Remark: The eigenvalues of $QA^{T}A$ lie in the region $(0,4/3)$, which guarantees the expected convergence of RP-ADMM. To obtain the expected convergence rate, we need to know the distance of the spectrum to the two extremes and $4/3$. We conjecture that the bound can be improved to $\rho(M)\leq 1-\frac{1}{2}\lambda_{\min}(A^{T}A)$. This requires more effort than the conjecture of RP-CD: besides showing $\tau_{\min}\geq O(\lambda_{\min}(A^{T}A)),$ we also need to show $\tau_{\max}\leq\frac{4}{3}-O(\lambda_{\min}(A^{T}A))$. This is left as future work.

Proof of Lemma 1

The proof of Lemma 1 relies on two simple techniques. The first technique, as elaborated in the Step 1 below, is to factorize $M$ and rearrange the factors. The second technique, as elaborated in the Step 2 below, is to reduce the dimension by eliminating a variable from the eigenvalue equation.

Step 1: Factorizing $M$ and rearranging the order of multiplication. The following observation is crucial: the matrix $M$ defined by (43) can be factorized as

Switching the order of the products by moving the first component to the last, we get a new matrix

Note that $\text{eig}(XY)=\text{eig}(YX)$ for any two square matrices, thus

Step 2: Relate the eigenvalues of $M^{\prime}$ to the eigenvalues of $QA^{T}A$, i.e. prove (62). This step is simple as we only use the definition of eigenvalues. However, note that, without Step 1, just applying the definition of eigenvalues of the original matrix $M$ may not lead to a simple relationship as (62).

We claim that (63) holds when $v_{1}=0$. In fact, in this case we must have $v_{0}\neq 0$ (otherwise $v=0$ cannot be an eigenvector). By (64b) we have $\lambda v_{0}=v_{0}$, thus $\lambda=1$. By (64a) we have $0=QA^{T}v_{0}=QA^{T}A(A^{-1}v_{0})$, which implies $\frac{(1-\lambda)^{2}}{1-2\lambda}=0\in\text{eig}(QA^{T}A)$, therefore (63) holds in this case.

The equation (64b) implies $(1-\lambda)v_{0}=Av_{1}$. Multiplying both sides of (64a) by $(1-\lambda)$ and invoking this equation, we get

We must have $\lambda\neq\frac{1}{2}$; otherwise, the above relation implies $v_{1}=0$, which contradicts (65). Then (66) becomes

Therefore, $\frac{(1-\lambda)^{2}}{1-2\lambda}$ is an eigenvalue of $QA^{T}A$, with the corresponding eigenvector $v_{1}\neq 0$, which finishes the proof of (63).

The other direction For the purpose of proving Theorem 2, we do not need to prove this direction. Here we present the proof since it is quite straightforward and makes the result more comprehensive.

is easy to prove. Suppose $\frac{(1-\lambda)^{2}}{1-2\lambda}\in\text{eig}(QA^{T}A)$. We consider two cases.

Case 2: $\frac{(1-\lambda)^{2}}{1-2\lambda}\neq 0$, then $\lambda\neq 1$. Let $v_{1}$ be the eigenvector corresponding to $\frac{(1-\lambda)^{2}}{1-2\lambda}$ (i.e. pick $v_{1}$ that satisfies (67)), and define $v_{0}=v_{1}/(1-\lambda)$. It is easy to verify that $v=\begin{bmatrix}v_{1}\\ v_{0}\end{bmatrix}$ satisfies $Mv=\lambda v$, which implies $\lambda\in\text{eig}(M)$.

Step 3: When $Q$ is symmetric, prove (45) by simple algebraic computation.

Note that when $\tau(\tau-1)<0$, the expression $\sqrt{\tau(\tau-1)}$ denotes a complex number $i\sqrt{\tau(1-\tau)}$, where $i$ is the imaginary unit. To prove (45), we only need to prove

Case 1: $\tau<0$. Then $\tau(\tau-1)=|\tau|(|\tau|+1)>0$. In this case, $\lambda_{1}=1+|\tau|+\sqrt{|\tau|(|\tau|+1)}>1.$

Case 2: $0<\tau<1$. Then $\tau(\tau-1)<0$, and (69) can be rewritten as

which implies $|\lambda_{1}|=|\lambda_{2}|=\sqrt{(1-\tau)^{2}+\tau(1-\tau)}=\sqrt{1-\tau}<1$.

Case 3: $\tau>1$. Then $\tau(\tau-1)>0$. According to (69), it is easy to verify $\lambda_{1}>0>\lambda_{2}$ and

Combining the conclusions of the three cases immediately leads to (70).

Proof of Lemma 2

This section is devoted to the proof of Lemma 2. We first give a proof overview in Section 6.1. The formal proof of Lemma 2 is given in Section 6.2. The proofs of the technical results involved in the proof are given in the subsequent subsections.

Without loss of generality, we can assume

In the proof overview, we discuss a few issues one may encounter when proving the result, and how we resolve these issues.

The simulations show that $\|QA^{T}A\|<\frac{4}{3}\ll\|Q\|\|A^{T}A\|$, thus we cannot relax $\|QA^{T}A\|$ to the product of $\|Q\|$ and $\|A^{T}A\|$, and have to treat $QA^{T}A$ as a single subject. However, each entry of $QA^{T}A$ is a complicated function (in fact, a high order polynomial) of the entries of $A^{T}A$. In other words, $Q$ is like a black box. To open the “black box”, we use a simple expression of $Z=I-AQA^{T}$ proved in Claim 4.1, i.e., $Z=E_{\sigma}(P_{\sigma_{1}}\dots,P_{\sigma_{n}}),$ where $P_{i}=I-A_{i}A_{i}^{T}$ is directly related to $A_{i}$. The problem becomes how to connect the eigenvalues of $E_{\sigma}(P_{\sigma_{1}}\dots,P_{\sigma_{n}})$ with those of $AA^{T}=\sum_{i}A_{i}A_{i}^{T}=n-\sum_{i}P_{i}$.

Although this is a clear linear algebra problem, it is not easy to obtain a lower bound of $E_{\sigma}(P_{\sigma_{1}}\dots,P_{\sigma_{n}})$. In fact, even though we know the eigenvalues of $Z=E_{\sigma}(P_{\sigma_{1}}\dots,P_{\sigma_{n}})$ are lower bounded by $-1$ because RP-CD converges, it is not clear how to prove this lower bound directly from a linear algebra perspective.

In our solution, we apply two tricks. The first trick is to view $E_{\sigma}(P_{\sigma_{1}}\dots,P_{\sigma_{n}})$ as an induction formula that connects it and its lower dimensional analogs. This is based on a simple observation that any permutation $(\sigma_{1}\sigma_{2}\dots\sigma_{n})$ can be written as the concatenation of $(\sigma_{1}\sigma_{2}\dots\sigma_{n-1})$ and $\sigma_{n}$, thus the expression of $Z=E_{\sigma}(P_{\sigma_{1}}\dots,P_{\sigma_{n}})$ can be decomposed accordingly. We then reduce the problem to bounding the eigenvalues of a Jordan product $P_{n}\hat{Z}+\hat{Z}P_{n}$, where $P_{n}$ is a projection matrix and $\hat{Z}$ is the lower dimensional analog of $Z$. The second trick is to apply a formula on the eigenvalues of Jordan product developed by Strang in 1962 . Somewhat surprisingly, his formula exactly leads to the desired lower bound of $-1/3$.

2 Proof of Lemma 2

The proof can be divided into three steps: first provide an alternative expression of $AQA^{T}$, then prove an induction formula, and finally apply Strang’s formula to perform mathematical induction. This subsection contains the major part of the proof, and the intermediate technical results will be proved in later subsections.

Step 0: Expression of $I-AQA^{T}$. As proved in Claim 4.1, we have a simple expression of the update matrix $I-AQA^{T}$

Define $W_{k}$ as the $k$-th block-column of $A^{T}A$ excluding the block $A_{k}^{T}A_{k}$, i.e.

Based on the expression of $I-AQA^{T}$ presented before, we build a connection between the update matrix $I-AQA^{T}$ and its lower dimensional analogs. The proof of Proposition 3 is given in Section 6.3.

where $Q$ is defined as in (39), $\hat{A}_{k}=[A_{1},\dots,A_{k-1},A_{k+1},\dots,A_{n}]$, and $\hat{Q}_{k}$ is defined in (73), and $P_{k}=I-A_{k}A_{k}^{T}$. Then we have

Step 2: Applying Strang’s result on Jordan product to perform mathematical induction.

It is obvious that the product of two symmetric matrices is not necessarily symmetric, so it is common to encounter the symmetrized product $XY+YX$, which is called Jordan product of two matrices $X$ and $Y$. Our induction formula basically states that $Z$ is the average of the Jordan product of the lower dimensional analog and $P_{k}$.

The eigenvalues of the Jordan product of two matrices have been studied before. The following result is proved in Strang .

([44, Theorem 1]; eigenvalues of Jordan product) Suppose two symmetric positive-semidefinite matrices $X$ and $Y$ satisfy

then the maximal (resp. minimal) eigenvalue of the Jordan product $XY+YX$ are the largest (resp. smallest) of the set

Let us come back to the proof of Lemma 2. We use mathematical induction to prove Lemma 2. For the basis of the induction ($n=1$), Lemma 2 holds since $QA^{T}A=I_{d_{1}\times d_{1}}$. Assume Lemma 2 holds for $n-1$, we will prove Lemma 2 for $n$.

Consider one term of (75) $P_{k}\hat{Z}_{k}+\hat{Z}_{k}P_{k}$. Note that $P_{k}=I-A_{k}A_{k}^{T}$ is a projection matrix, since we have assumed $A_{k}^{T}A_{k}=I$. Combining with the induction hypothesis, we have

Let $\alpha_{1}=0,\alpha_{n}=1,\beta_{1}=-1/3,\beta_{n}=1$, then the set (76) becomes (keep the repeated values)

Note that since by the induction hypothesis the eigenvalues of $\hat{Z}_{k}$ cannot achieve the extreme values of region $(-1/3,1)$, the eigenvalues of $\frac{1}{2}(P_{k}\hat{Z}_{k}+\hat{Z}_{k}P_{k})$ also cannot A more detailed argument is as follows. Since $-I/3\preceq\hat{Z}_{k}$, we can let $\beta_{1}=-1/3+\epsilon$ for a sufficiently small positive number $\epsilon$, while keeping $\alpha_{1}=0,\alpha_{n}=1,\beta_{n}=1$. The set (76) now becomes $\{0,0,-2/3+2\epsilon,2,-\frac{(4/3-\epsilon)^{2}}{4(2/3+\epsilon)}\}.$ Both $-2/3+2\epsilon$ and $-\frac{(4/3-\epsilon)^{2}}{4(2/3+\epsilon)}$ are strictly larger than $2/3$, thus the extreme value $-2/3$ cannot be achieved. By a similar argument the other extreme value $2$ also cannot be achieved. . So we have

Remark: Where does the magical number $-1/3$ come from? It is actually the strange and complicated term $\frac{16\alpha_{1}\alpha_{n}\beta_{1}\beta_{n}-(\beta_{1}-\beta_{n})^{2}(\alpha_{1}-\alpha_{n}^{2})}{4(\alpha_{1}+\alpha_{n})(\beta_{1}+\beta_{n})}$ in Strang’s result (76), which occurs due to the special structure of the Jordan product.

3 Proof of Proposition 3 (the induction formula)

It is easy to build an induction formula from the expression (51b). For example, when $n=3$, the matrix $\sum_{\sigma}P_{\sigma_{1}}P_{\sigma_{2}}P_{\sigma_{3}}$ can be decomposed as the sum of $P_{1}(P_{2}P_{3}+P_{3}P_{2})+(P_{2}P_{3}+P_{3}P_{2})P_{1}$ and two other similar terms (changing the outside part $P_{1}$ to $P_{2},P_{3}$ and the inside part $P_{2}P_{3}+P_{3}P_{2}$ correspondingly). The inside part $P_{2}P_{3}+P_{3}P_{2}$ only involves two matrices, thus is a lower-dimensional analog. To make this even easier to see, denote $X=P_{1},Y=P_{2},Z=P_{3},$ then

A rigorous argument based on the above intuition is given as follows. Applying the formula (51b) to the matrix $P_{1},\dots,P_{k-1},P_{k+1},\dots,P_{n}$, and by the definition $\hat{A}_{k}=[A_{1},\dots,A_{k-1},A_{k+1},\dots,A_{n}]$ and the definition of $\hat{Q}_{k}$ in (73), we have

4 Proof of Proposition 1

We provide the proof of the expected convergence of BR-ADMM here, as this proof is a slightly smaller subset of the proof of Theorem 1. We will just describe the necessary modifications.

We only need to prove a similar version of Theorem 2, i.e., the spectral radius of the expected update matrix of BR-ADMM is less than 1. Throughout the proof, we need to change the matrix $Q=\frac{1}{|\Gamma|}\sum_{\sigma\in\Gamma}Q_{\sigma}$ to another one defined as

where $\Gamma^{\text{BR}}$ denotes the set of all possible permutations according to the Bernoulli randomization rule. It is easy to see that $|\Gamma^{\text{BR}}|=2^{n}$. Other matrices such as $M$ should be changed accordingly.

The proof of Theorem 2 mainly consists of Lemma 1 and Lemma 2. Since Lemma 1 has nothing to do with the specific expression of $Q$, so we only need to prove Lemma 2 for BR-ADMM, i.e., the matrix $AQ^{\text{BR}}A^{T}$ has all eigenvalues in the region $(0,4/3)$. Following the proof of Lemma 2, we divide the proof into three steps.

Step 0: Expression of $Z^{\text{BR}}\triangleq I-AQ^{\text{BR}}A^{T}$. In Claim 4.1, we have prove the expression (51a) that $I-AQ_{\sigma}A^{T}=P_{\sigma_{n}}P_{\sigma_{n-1}}\dots P_{\sigma_{1}}$ for any permutation $\sigma$, which implies

Step 1: Induction formula. Notice that a characteristic of the Bernoulli randomization rule is: the first block is either updated first or updated last. For instance, when $n=4$, $(1,3,4,2)$ is a feasible permutation in $\Gamma^{\text{BR}}$ and $(3,4,2,1)$ is also a feasible permutation, but $(3,1,4,2)$ is not feasible. After removing the first block, the rest $n-1$ blocks form a permutation in $\hat{\Gamma}_{\text{BR}},$ where $\hat{\Gamma}_{\text{BR}}$ is the set of all permutation of $2,3,\dots,n$ according to the Bernoulli randomization rule. In other words, we have $\Gamma^{\text{BR}}=\{(1,\hat{\sigma}),(\hat{\sigma},1),\text{ where }\hat{\sigma}\in\hat{\Gamma}^{\text{BR}}\}$. Thus we have an induction formula

where $\hat{Z}^{\text{BR}}$ is the lower dimensional analog of $Z^{\text{BR}}$ for the rest $n-1$ blocks (after removing the first block).

Step 2: Applying mathematical induction. This step is almost the same as Step 2 of the proof of Lemma 2. More specifically, combining the induction hypothesis that $\text{eig}(\hat{Z}^{\text{BR}})\in(-1/3,1),$ Strang’s result Lemma 5 and (78), we obtain the desired result $\text{eig}(Z^{\text{BR}})\in(-1/3,1).$ This finishes the proof.

Proof of Technical Results for Expected Convergence Rates

Suppose all the distinct eigenvalues of $I-QA^{T}A$ are $0<\tau_{N^{\prime}}<\dots<\tau_{1}<4/3$, where $1\leq N^{\prime}\leq N$. Denote $\tau_{\min}=\tau_{N^{\prime}},\tau_{\max}=\tau_{1}.$ According to Lemma 1, the expected update matrix of RP-ADMM $M$ has $2N^{\prime}$ distinct eigenvalues $\lambda_{k,1},\lambda_{k,2}$ given by