Faster proximal algorithms for matrix optimization using Jacobi-based eigenvalue methods
Year: 2022

Abstract:
We consider proximal splitting algorithms for convex optimization problems over matrices. A significant computational bottleneck in many of these algorithms is the need to compute a full eigenvalue or singular value decomposition at each iteration for the evaluation of a proximal operator. In this paper we propose to use an old and surprisingly simple method due to Jacobi to compute these eigenvalue and singular value decompositions, and we demonstrate that it can lead to substantial gains in terms of computation time compared to standard approaches. We rely on three essential properties of this method: (a) its ability to exploit an approximate decomposition as an initial point, which in the case of iterative optimization algorithms can be obtained from the previous iterate; (b) its parallel nature which makes it a great fit for hardware accelerators such as GPUs, now common in machine learning, and (c) its simple termination criterion which allows us to trade-off accuracy with computation time. We demonstrate the efficacy of this approach on a variety of algorithms and problems, and show that, on a GPU, we can obtain 5 to 10x speed-ups in the evaluation of proximal operators compared to standard CPU or GPU linear algebra routines. Our findings are supported by new theoretical results providing guarantees on the approximation quality of proximal operators obtained using approximate eigenvalue or singular value decompositions.

heading: 1 Introduction
text: Many problems in statistics and learning can be formulated as convex optimization over matrices. This includes kernel learning [LCB+04], low-rank matrix recovery [RFP10, CR09], sparse covariance selection [BGdN06, FHT08], community detection [BBV16], and many others. Proximal splitting algorithms such as (accelerated) proximal gradient methods [Nes03, BT09], or the alternating direction method of multipliers (ADMM) [EB92, BPC11], are among the most popular methods to deal with large-scale instances of such problems. These methods solve generic convex optimization problems of the form
min X∈E F (X) +G(X) (1)
via calls to the proximal operators of F and/or G. When the variable X is a symmetric or rectangular matrix, it is often the case that the functions F and/or G are symmetric functions of the eigenvalues or singular values of X; i.e., so-called spectral or othogonally-invariant functions. Typical examples are F (X) = I(X 0) = ∑n i=1 I(λi(X) ≥ 0) the convex indicator function for the positive
semidefinite cone, or F (X) = ‖X‖∗ = ∑n i=1 σi(X) the nuclear norm.
Computing the proximal operator of such spectral functions F requires a full eigenvalue/singular value decomposition which contributes a large computational cost of at least n3 floating-point
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
operations per iteration. The most common algorithms for such decompositions are the QR method, and the divide-and-conquer algorithm [GVL13] which are implemented in linear algebra packages (LAPACK) that are used in standard numerical computing software such as MATLAB or NumPy.
Contributions In this paper we propose to use a simple and lesser known method due to Jacobi from 1846 [Jac46] for the computation of these decompositions. Several features of this method make it suitable for its use within iterative proximal methods:
(i) It can effectively exploit an initial approximate decomposition of the target matrix. In iterative proximal methods, the eigenvectors/singular vectors from the previous iterate are often a very good initial point that almost diagonalize the matrix of the current iterate.
(ii) The method can be very easily parallelized, and is a very good fit for parallel hardware accelerators, notably GPUs, which are now very common in machine learning.
(iii) It has a simple termination criterion that allows us to trade-off computation time with the accuracy of the desired decomposition.
The present paper has two main contributions. First, we establish new theoretical results providing guarantees on the approximation quality of the proximal operator of a convex spectral function when computed with an approximate eigenvalue/singular value decomposition. When combined with existing results on the convergence of proximal splitting methods with inexact proximal oracles, we are able to appropriately set the tolerance of Jacobi-based methods and gain computation time, while guaranteeing global convergence of the optimization method. Second, we demonstrate using the observations above, and on a variety of problems and splitting algorithms, that the method of Jacobi on a GPU can yield speed-ups of 5-10x in the evaluation of proximal operators of spectral functions, compared to standard approaches on CPU and GPU.
Related work Several approaches have been considered to avoid computing a full eigenvalue or singular value decomposition at each iteration, when the solution is known to be low-rank. Popular approaches include the Lanczos algorithm, see e.g., [CCS10], randomized SVD [HMT11], or nonconvex low-rank formulations [BM03]. However these approaches usually rely on an a priori knowledge of the rank of the solution.
The recent paper [RGN19] focusses on the ADMM method for SDPs and proposes the use of iterative eigensolvers (LOBPCG [Kny01]) to approximately project a symmetric matrix on the positive semidefinite cone. The proposed method can achieve 10x speedups, especially when the primal or dual solution has low-rank. From the theoretical point of view, the approach relies on error bounds derived in [GNR20] concerning approximate projections on the positive semidefinite cone. Our theoretical results in Section 4 constitute a significant generalization of these error bounds to proximal operators of convex spectral functions, and to different notions of approximations.
The papers [CV18] and [JV21] consider proximal splitting algorithms with a generalized Bregman divergence (i.e., different than the squared Euclidean distance), and show that these proximal operators can sometimes be evaluated without explicitly computing a full eigenvalue decomposition. If φ is a strongly-convex function, the φ-proximal operator of F is defined by
proxφF (y; a) = argminx {F (x) + 〈a, x〉+ φ(x)− 〈∇φ(y), x〉} .
The paper [CV18] considers optimization problems over symmetric positive semidefinite Toeplitz matrices and shows that by choosing φ to be the kernel for the so-called Itakura-Saito distance one can compute the generalized proximal operator in O(n2) time, instead of O(n3). In [JV21] it is shown that for sparse semidefinite programs, an adequate choice of φ leads to a proximal operator that can be computed at the cost of a Cholesky decomposition of a sparse matrix which in general is much smaller than the cost of a full eigenvalue decomposition. Concerning sparse semidefinite programs, there has been a lot of work on devising fast ADMM algorithms using the clique decomposition method [FKMN01], e.g., [ZFP+20, MKL15]. We note that these contributions are specific to sparse semidefinite programs and are complementary to the present paper because they still require an eigenvalue decomposition to be computed at each iteration.
Finally let us mention that there are first-order methods that do not require a full eigenvalue decomposition at each iteration but only the computation of the smallest eigenvalue(s)/eigenvector(s) at each iteration such as e.g., [HR00, Ren14, YTF+21] and for which Lanczos methods are more suited.
Organization In Section 2 we recall some necessary background material about convex spectral functions and their proximal operators. In Section 3 we present the Jacobi method to diagonalize symmetric matrices, and explain how it can be used with proximal algorithms. In Section 4 we present the main theoretical result of the paper and show that approximate eigenvalue/singular value decompositions yield guaranteed approximations of proximal operators for spectral functions. Finally we present numerical experiments in Section 5.

heading: 2 Background
text: In this section we present some notations and background results concerning spectral functions. If f : E → R is a convex function on Euclidean space E we let ∂f(x) be the subdifferential of f at x, namely ∂f(x) = {g ∈ E : f(y) ≥ f(x) + 〈g, y − x〉}. An element g ∈ ∂f(x) is called a subgradient of f at x. The proximal operator of f at y is defined as
proxf (y) = argminx∈E { f(x) + (1/2)‖x− y‖2 } .
We let Sn be the Euclidean space of n×n symmetric matrices with the trace inner product 〈X,Y 〉 = tr(XY ). We let Sn+ (resp. S n ++) be the cone of positive semidefinite (resp. definite) matrices in S
n. Also O(n) denotes the set of n× n orthogonal matrices.
Spectral functions A function F : Sn → R of symmetric matrices is called spectral if there is a permutation-invariant function f : Rn → R such that F (X) = f(λ(X)) where λ(X) = (λ1(X), . . . , λn(X)) is the vector of eigenvalues of X (in any order). It is clear that spectral functions are orthogonally-invariant, i.e., F (UTXU) = F (X), because UTXU and X have the same eigenvalues; conversely any orthogonally-invariant function F is the spectral function associated to the permutation-invariant function f(x) = F (diag(x)). The main remarkable fact that we will use about spectral functions is a result due to Davis [Dav57], which asserts that F is convex as soon as f is, and gives a variational formula for F certifying this convexity.
Theorem 1 ([Dav57], see also [BTN01, p.282, claim D]). Let f : Rn → R be a permutationinvariant function, and let F : Sn → R, F (X) = f(λ(X)) be the associated spectral function. If f is convex on Rn, then F is convex on Sn, and we have the following variational formulation for F :
F (X) = max U∈O(n)
f(diag(UTXU)). (2)
An important, and well-known consequence of the theorem above, is that the proximal operator of F can be obtained from that of f , namely for any Y ∈ Sn we have
proxF (Y ) = U diag(proxf (y))U T
where Y = U diag(y)UT is an eigenvalue decomposition of Y .
Functions of singular values Similar results hold for functions F : Rm×n → R of rectangular matrices that are orthogonally-invariant, i.e., F (UXV ) = F (X) for all U ∈ O(m) and V ∈ O(n). Any such function can be written as F (X) = f(σ(X)) where σ(X) ∈ Rmin(m,n) is the vector of singular values and f : Rmin(m,n) → R is an absolutely symmetric function, i.e., it is invariant under signed permutations: f(η1xπ(1), . . . , ηkxπ(k)) = f(x1, . . . , xk) for ηi = ±1 and π a permutation of {1, . . . , k} where k = min(m,n) [Lew95]. In this case too, it holds that F is convex as soon as f is, and we have:
F (X) = max U∈O(m),V ∈O(n)
f(diag(UTXV )). (3)
The variational formulas (2) and (3) are a consequence of well-known results concerning majorization relations between the diagonal of a matrix and its eigenvalues or singular values [MOA79, Chapter 9]. For completeness we have included a proof of these formulas in the supplementary material. We refer the reader to [PB14, Section 6.7] and references therein for more details on the proximal operator of spectral and orthogonally-invariant functions.

heading: 3 Jacobi-based methods for eigenvalue and singular value decompositions
text: In this section we present Jacobi-based methods to compute eigenvalue and singular value decompositions, and explain how they can be used in proximal optimization algorithms. For concreteness we will focus on the eigenvalue decomposition of symmetric matrices, but as will be clear, the extension to the SVD of rectangular matrices will be straightforward.

heading: 3.1 The method
text: Jacobi eigenvalue methods compute the eigenvalue decomposition of a real symmetric matrix by iteratively eliminating off-diagonal elements using Givens rotations. More precisely, let A be a n × n real symmetric matrix that we want to diagonalize. The Jacobi method begins by selecting an off-diagonal index (i, j), with i < j, such that |Aij | > 0, and constructs a 2× 2 rotation matrix[
cos θ sin θ − sin θ cos θ
] that diagonalizes the 2× 2 submatrix [ Aii Aij Aij Ajj ] . The rotation matrix can be embedded
in a larger n × n Givens rotation matrix at indices (i, j) which is otherwise equal to the identity matrix: {
Gii = Gjj = cos θ Gij = sin θ,Gji = − sin θ and Gpq = δp,q for p /∈ {i, j} or q /∈ {i, j},
where δp,q = 1 if p = q and 0 otherwise. This Givens rotation matrix satisfies (GAGT )ij = 0. The Jacobi algorithm then proceeds by updating A← GAGT and performs the same process until a convergence criterion is satisfied. Usually the convergence criterion is given in terms of off(A), defined to be the square root of the sum of the off-diagonal elements of A:
off(A) := √ ∑ 1≤i 6=j≤n A2ij . (4)
The algorithm can thus be summarized as follows:
Algorithm 1: Jacobi method for diagonalizing symmetric matrix A while off(A) ≥ tolerance do
choose an off-diagonal index (i, j) such that |Aij | > 0 construct a Givens matrix G such that (GAGT )ij = 0 update A← GAGT
end
In the classical Jacobi method, the off-diagonal entry (i, j) selected at each iteration is the one for which |Aij | is largest. It can be easily shown that with this choice, the quantity off(A) satisfies after one iteration
off(GAGT )2 ≤ (1− 1/N) off(A)2
where N = n(n− 1)/2. As such off(A) converges linearly to 0. It has been shown by Schönhage, and van Kempen [Sch64, VK66a] that the method actually attains quadratic convergence locally.
Cyclic method The classical Jacobi method requires searching for a largest element among the off-diagonal elements at each iteration. For an n×nmatrix the cost of this step isO(n2) flops. Cyclic Jacobi methods instead cycle through all off-diagonal elements of the matrix in some fixed order, computing the associated Givens matrices and updating accordingly. The algorithm proceeds by sweeps, where in each sweep each off-diagonal entry is processed exactly once. Various convergence results have been proved for cyclic methods most notably local quadratic convergence, see e.g., [FH60, Wil62, VK66b]. In practice, the number of sweeps needed to diagonalize a matrix of size n . 1000 is usually no more than 15.
Parallelization The cyclic Jacobi method can be easily parallelized. Indeed, instead of eliminating a single off-diagonal entry at each iteration, one can eliminate in parallel bn/2c off-diagonal entries that correspond to distinct rows/columns. More precisely by selecting indices (i1, j1), . . . , (im, jm)
where m = bn/2c such that i1, j1, i2, j2, . . . , im, jm are all distinct one can compute in parallel the m Givens rotations that are needed to eliminate the entries Ai1j1 , . . . , Aimjm . A single sweep can then be performed in n − 1 steps where at each step k = 1, . . . , n − 1 one chooses indices (i1[k], j1[k]), . . . , (im[k], jm[k]) that are pairwise distinct in such a way that all the off-diagonal entries are processed at least once. For more details on the choice of these indices we refer the reader to [GVL13, Section 8.5].
Other eigenvalue algorithms In general, the Jacobi method is not competitive with more traditional eigenvalue algorithms such as the QR method or the divide-and-conquer method. Counting floating point operations only, the QR method is roughly equivalent to two sweeps of Jacobi [GVL13, TBI97]. One important difference between the Jacobi method and QR/divide-and-conquer algorithms is that the latter start by putting the matrix A into tridiagonal form by a sequence of orthogonal (Householder) transformations, whereas the Jacobi method operates on A directly. One consequence of this is that the Jacobi method can have faster convergence if off(A) is small, whereas this is not the case of other algorithms that start by a reduction into tridiagonal form. This ability to exploit an approximate diagonalization of A is what makes the Jacobi method particularly interesting for iterative optimization algorithms.
Singular value decomposition The Jacobi algorithm can be very easily adapted to compute the singular value decomposition of rectangular matrices. The idea is to solve a sequence of 2 × 2 SVD problems, instead of 2× 2 symmetric eigenvalue problems. This leads to the two-sided Jacobi SVD algorithm. Another approach, called the one-sided Jacobi SVD algorithm, performs orthogonal transformations on the right to make the columns of A orthogonal, by iteratively orthogonalizing pairs of columns. From there, one can extract the singular values as the Euclidean norms of the columns of A, and the left singular vector matrix easily. We refer the reader to [GVL13, Section 8.6.4] for more details and references.

heading: 3.2 Jacobi-based methods in proximal algorithms
text: We now explain how the Jacobi method can be used in iterative optimization algorithms for the computation of proximal operators of spectral functions. The main idea is to use the eigenvector matrix from the previous iteration as a warm-start for the current iterate. More precisely, assume that at iteration k one needs to compute the proximal operator proxF (Yk) where F : S
n → R is the spectral function associated to the permutation-invariant function f : Rn → R. We proceed as follows:
• Let Uk−1 be the orthogonal eigenvector matrix computed at iterate k − 1 • Apply the Jacobi method to UTk−1YkUk−1 with tolerance k. This returns an orthogonal
matrix W such that off(WTUTk−1YkUk−1W ) ≤ k
• Set Uk = Uk−1W (eigenvector matrix of Yk) and yk = diag(UTk YkUk) (eigenvalues of Yk) • Output Uk diag(proxf (yk))UTk as the proximal operator proxF (Yk).
As we demonstrate in the numerical experiments section below, a very small number of sweeps (often just one) are in general needed to converge to the eigenvalue decomposition of Yk in the second step above, because the matrix UTk−1YkUk−1 is usually very close to being diagonal, especially after the first few iterations of the optimization algorithm.
It is straightforward to extend the warm-start approach described here to the situation where we want to compute singular value decompositions, and we omit the details for brevity. Remark 1. When numerically implementing the procedure above, the accumulation of errors can cause the columns of Uk not to be normalized after a large number of iterations. To resolve this issue one can orthonormalize the columns of Uk periodically using a QR decomposition.

heading: 4 Approximate proximal evaluation
text: Most proximal algorithms (e.g., ADMM, proximal gradient method, ...) are still convergent even if only an approximation of the proximal operator is computed at each iteration, provided the approxi-
mation quality suitably decays to zero, see e.g., [EB92, SRB11, Yao16, BTB20]. A natural question is to understand how the tolerance parameter of the Jacobi method k affects the approximation quality of the proximal operator proxF (Yk)? This section addresses this question.
Measures of quality of prox There are multiple ways to quantify the approximation quality of a proximal operator, see e.g., [BTB20] which compiles a list of different ways that were considered in the literature. In this section we will focus on two simple notions that have been used in convergence guarantees for proximal splitting methods. Let F : E → R be a convex function for which we want to compute the prox at a point Y ∈ E. Let
Φ(X;Y ) = F (X) + (1/2)‖X − Y ‖2F
be the prox objective so that proxF (Y ) = argminX∈E Φ(X;Y ). If X̂ ∈ E is a candidate value for proxF (Y ) the two measures of approximations we consider are: (i) Euclidean distance to the true prox, i.e., δ = ‖X̂ − proxF (Y )‖F and (ii) suboptimality in prox objective value, i.e., η = Φ(X̂;Y ) − Φ(proxF (Y );Y ). The former notion was used e.g., to guarantee convergence of ADMM algorithms, assuming that the sequence of errors (δk) is summable [EB92, Yao16]; and the latter was used in [SRB11] to guarantee convergence of (fast) proximal gradient methods, under a decaying assumption on (ηk). Remark 2. Note that since Φ is 1-strongly convex we always have the inequality:
Φ(X̂;Y )− Φ(proxF (Y );Y ) ≥ 1
2 ‖X̂ − proxF (Y )‖2. (5)
This means that any bound on the approximation quality in terms of prox objective value, implies a bound on the approximation quality using the Euclidean norm.
Approximate prox using approximate eigenvalue/singular decomposition Our main result gives bounds on the quality of the prox of any convex spectral function computed using an approximate eigenvalue decomopsition. More precisely, assume that F : Sn → R is a convex spectral function of symmetric matrices associated to the permutation-invariant function f(x1, . . . , xn). Let Y ∈ Sn for which we want to compute proxF (Y ). We know that proxF (Y ) = U proxf (y)UT where y = (y1, . . . , yn) is the vector of eigenvalues of Y , and U is the orthogonal matrix that diagonalizes Y . Assume that Û is an orthogonal matrix such that off(ÛTY Û) ≤ , i.e.,
ÛTY Û = diag(ŷ) + E where { ‖E‖F ≤ Eii = 0 for i = 1, . . . , n.
(6)
Our main theoretical result, stated below, quantifies the error in the approximation X̂ ≈ proxF (Y ) where
X̂ = Û diag(proxf (ŷ))Û T . (7)
Theorem 2. Under the assumptions above, we have: (i) ‖X̂ − proxF (Y )‖F ≤ ; and (ii) Φ(X̂;Y )− Φ(proxF (Y );Y ) ≤ 2/2.
A similar result can be stated for orthogonally-invariant functions of rectangular matrices. Theorem 3. Let F : Rm×n → R be a convex orthogonally-invariant function associated with the absolutely-symmetric function f(x) = F (diag(x)) defined on Rmin(m,n). Let Y ∈ Rm×n, and assume that Û ∈ O(m) and V ∈ O(n) are orthogonal matrices such that
ÛTY V̂ = diag(ŷ) + E where { ‖E‖F ≤ Eii = 0 for i = 1, . . . ,min(m,n).
Then, with X̂ = Û diag(proxf (ŷ))V̂ T , the following holds: (i) ‖X̂ − proxF (Y )‖F ≤ , and (ii) Φ(X̂;Y )− Φ(proxF (Y );Y ) ≤ 2/2.
The proofs of both theorems are in the supplementary material. Some remarks are in order:
• Inequality (i) generalizes the main result in [GNR20].1 The main result in [GNR20] shows that when F is the indicator function of the positive semidefinite cone (so that proxF is the Euclidean projection on Sn+) then ‖X̂ − proxF (Y )‖F ≤ √ 2 . Our result significantly
generalizes the latter to the proximal operator of arbitrary convex spectral functions F . • Note that inequality (i) actually follows from (ii) by (5).
We believe Theorems 2 and 3 are remarkable in a number of ways:
• They bound the distance to the true prox (a quantity that is hard to compute) using a simple quantity which is = off(ÛTY Û).
• the bounds apply to all convex spectral functions F and do not depend on any Lipschitz assumption on F ;
• the bounds do not depend on how the orthogonal matrix Û is obtained; in this paper we consider the Jacobi method for obtaining approximate diagonalizations, however the bounds can be used with any other method;
• well-known perturbation bounds for eigenvalue/eigenvector computation (e.g., the DavisKahan theorem) suggest that the error would depend on additional factors related to the spectrum of Y . Such a dependence would make the theorem very hard to apply. Our bounds are free of such dependence and only depend on .
Implications If we combine our theorems with [EB92] and [SRB11] we get guarantees on the convergence rates of ADMM and (fast) proximal gradient methods, when the proximal operators are evaluated using an approximate eigenvector matrix. To illustrate this, consider for instance the following optimization over symmetric matrices
min X∈Sn F (X) +G(X) (8)
where F (X) is a convex spectral function and G(X) is convex and smooth with L-Lipschitz gradient. Consider the (accelerated) proximal gradient method
Xk+1 = proxL−1F (Z − (1/L)∇G(Z)) (9)
where Z = Xk for the standard proximal gradient method, and Z = Xk + k−1k+2 (Xk − Xk−1) for the accelerated version. Assume that the proximal operator in (9) is evaluated approximately according to (6)-(7) with accuracy k. Then, using [SRB11, Props. 1 and 2] and Theorem 2(ii), we immediately get: (a) if ( k) is summable (e.g., k = O(1/k1+δ)) then the (basic) proximal gradient method converges at a rate O(1/k) to the optimal value of (8); and (b) if (k k) is summable (e.g., k = O(1/k
2+δ)) then the accelerated proximal gradient method converges at a rate O(1/k2) to the optimal value of (8). For the ADMM algorithm, we can combine the result of [EB92, Theorem 8] with Theorem 2(i) to guarantee convergence when ( k) is summable.

heading: 5 Numerical experiments
text: In this section we present numerical experiments with Jacobi-based methods on different problems and splitting algorithms.
Setup We use the parallel version of the Jacobi method implemented on GPU, and available through the NVIDIA CUSOLVER library [NVI20] in the functions syevj and gesvdj for symmetric eigenvalue decomposition and singular value decomposition respectively. We compare it with the standard functions on both CPU/GPU implemented in LAPACK/CUSOLVER’s functions syevd and gesvd. All our code is in Python and the latter functions are invoked respectively by NumPy’s [HMvdW+20] and CuPy’s [OUN+17] functions linalg.eigh and linalg.svd. The GPU we used is a NVIDIA TITAN Xp with CUDA 11.0, and the CPU is an Intel Core i5-6500 with 4 cores at 3.20 GHz and 32GB RAM, running Ubuntu 18.04. For GPU computations (default eigenvalue solver and Jacobi’s method) we use 32-bit floats, instead of 64-bit doubles, which are much faster on GPU. On the examples we have tried this had little or no effect on the quality of the solutions, especially because we are looking at problems where moderate accuracy is satisfactory.
1We note that the proof of inequality (i) does not actually use the assumption that Eii = 0 for i = 1, . . . , n.
Performance of Jacobi’s method on GPU Before showing the results on optimization problems, we first compare in Table 1 the time needed for a single Jacobi sweep vs. the time needed for a full eigenvalue (resp. singular value) decomposition using NumPy/CuPy’s linalg.eigh (resp. linalg.svd), on a random Gaussian n × n matrix. For symmetric eigenvalue decompositions, we note that for matrices of size ∼ 500, one sweep of Jacobi is about 10x faster than the default methods. However the default method becomes faster for matrices of size 2000.2 On the other hand, for singular value decomposition the ratio is more pronounced and attains 40x for matrices of size & 2000.

heading: 5.1 Sparse inverse covariance selection
text: We consider here the well-known statistical problem of estimating an unknown covariance matrix Σ ∈ Sn++ from samples x1, . . . , xN ∈ Rn of a Gaussian distribution N (0,Σ), under the assumption that Σ−1 is sparse. The approach proposed in [BGdN06, FHT08] consists of solving the following convex optimization problem:
min X∈Sn++
tr(ŜX)− log detX + λ‖X‖1 (10)
where λ > 0 is a regularization parameter, Ŝ ∈ Sn+ is the given empirical covariance matrix, and ‖X‖1 = ∑ ij |Xij | is the entrywise `1-norm of X . The ADMM algorithm is one of the most popular methods for this problem, see e.g., [SMG10, BPC11] and references therein. The problem (10) can be naturally put in the splitting form (1) with F (X) = tr(ŜX) − log detX and G(X) = λ‖X‖1. The ADMM iterates can be written as: Xk+1 = prox−ρ log det(Zk − Uk − ρŜ) Zk+1 = Sλ/ρ(Xk+1 + Uk)
Uk+1 = Uk +Xk+1 − Zk+1
(11)
where the proximal operator of −ρ log det is prox−ρ log det(Y ) = V diag ( (y + √ y2 + 4ρ)/2 ) V T
where Y = V diag(y)V T is the eigenvalue decomposition of Y ; and St is the entrywise softthresholding operation defined by St(x) = 0 if |x| ≤ t and x − sign(x)t otherwise. The main computational cost in the iterations (11) is the proximal update rule for Xk which requires the eigenvalue decomposition of Zk − Uk − ρŜ at each iteration. We can use the Jacobi method to compute this eigenvalue decomposition, using the previous iterate’s orthogonal transformation as a warm start, as suggested in Section 3.2.
We have implemented this approach with different values of n in the range 200-1000, where the covariance matrix Σ ∈ Sn++ was chosen to be a banded matrix with constant bandwidth k = 5. We have used the parameters and stopping criterion of the Matlab implementation of (11) from [BPC11] available at https://web.stanford.edu/~boyd/papers/admm/. We compare the Jacobi method
2It is not clear to us why the performance of the Jacobi-symmetric eigenvalue routine syevj significantly degrades for sizes 2000. Note that, surprisingly, the Jacobi-SVD routine on matrices of size 4096 is faster than the Jacobi-symmetric eigenvalue routine on the same size.
with the standard routines in linear algebra packages for eigenvalue decomposition on CPU and GPU (NumPy/CuPy’s linalg.eigh).
The results of the numerical experiments are shown in Table 2. Depending on the size of the matrix, the Jacobi method can be more than 5x faster than the standard approaches on both CPU and GPU. For the Jacobi method we used a relative tolerance that decays like 1/k2, more precisely we set k = ‖Yk‖F /(1 + k2), where Yk is the matrix to diagonalize at iteration k. We see that the average number of sweeps per iteration is ≈ 1 confirming the fact that the eigenvector matrix from iterate k − 1 is a good initial point for the k’th iterate.

heading: 5.2 Semidefinite programming
text: We consider here semidefinite programs in primal-dual form,
min x∈Rn,S∈Sk
cTx s.t. ∑n i=1 xiAi + S = B, S 0
max Y ∈Sk − tr(BY ) s.t. Y 0, 〈Ai, Y 〉 = −ci (i = 1, . . . , n)
(12) where c ∈ Rn, A1, . . . , An ∈ Sk and B ∈ Sk. A primal-dual ADMM-based algorithm was proposed in [OCPB16] and implemented in the SCS package [OCPB19]. Each iteration of the algorithm consists of two steps, one projection onto a linear subspace, and one projection onto the positive semidefinite cone Sk+. The projection step can be computed using the Jacobi-based method. Table 3 shows the time per iteration required for the projection step using Jacobi’s method compared to the default approaches, on random instances of (12). Like in the previous example, we used a relative tolerance for the Jacobi method that decays like 1/k2, more precisely, k = ‖Yk‖F /(1 + k2), where Yk is the matrix to diagonalize at iteration k. We see that the time spent in the computation of the projection on the positive semidefinite is up to 5x faster than the default approaches, and that the number of sweeps per iteration is close to 1.

heading: 5.3 Nuclear norm minimization using the proximal gradient method
text: We consider the low-rank matrix recovery problem
min X∈Rn×n
(1/2)‖A(X)− b‖22 + λ‖X‖∗ (13)
where A : Rn×n → Rm, b ∈ Rm, λ is the regularization parameter, and ‖X‖∗ = ∑n i=1 σi(X) is the nuclear norm. A popular method for solving (13) is the (fast) proximal gradient method, with iterates:
Xk+1 = SVTtλ(Z − tA∗(A(Z)− b)) (14) where t > 0 is the step size, and SVTtλ is the proximal operator associated to the nuclear norm which is given by the soft-thresholding of singular values. In the standard proximal gradient method, we take Z = Xk in (14), whereas for the fast proximal gradient method we take Z = Xk+ k−1k+2 (Xk−Xk−1). Clearly, the main computational cost in the implementation of (14) is to compute the singular value decomposition of a n× n matrix. Table 4 shows the result of using the fast proximal gradient method with the Jacobi-based method vs. the default SVD routines on CPU/GPU for the computation of the proximal operator, on random matrix-completion instances of (13). Like for the previous examples, we used a relative tolerance for the Jacobi method decaying like 1/k2, except for the largest instance where we fixed the number of sweeps to be 1 at each iteration. The number of iterations (14) was fixed and set to 1000, with a time limit of 1 hour on the total running time. For the largest instance the two default methods exceeded the time limit and the objective value we report is the one obtained after 1000 iterations.

heading: Acknowledgments
text: H.G. is supported by a GlaxoSmithKline PhD studentship. H.F. would like to thank Oisín Faust and Alhussein Fawzi for comments on an earlier draft of this paper.

Review summaries:
This work proposes to apply a Jacobi-based approximate eigenvalue decomposition method in proximal slitting algorithms for minimization of convex spectral function with matrix variable. The adopted Jacobi-based method is used to compute proximal operators in each iteration, and by taking the previous iterate as the initial point, the Jacobi-based method appears to reach required accuracy within fewer sweeps. The method can be parallelized, and its GPU implementation shows improvement in efficiency compared with standard eigenvalue decomposition functions in LAPACK and CUSOLVER in numerical experiments.

This paper proposes using a classical but often unused method for computing the eigenvalue decomposition of a symmetric matrix: the Jacobi method. The method simply iterates by zeroing out off-diagonal entries using Givens rotations one at a time, until the matrix converges to a diagonal matrix. It is further noted that it is possible to parallelize the removal of ~n/2 entries at a time if disjoint indices are picked (for a matrix of size n-by-n).  Proximal methods are a popular choice for solving convex optimization problems over matrices. A commonly required oracle for such problems is the computation of the (full) eigenvalue decomposition of the current iterate, e.g., to compute the projection of the iterate onto the positive semidefinite cone. It is observed that as the iterations progress, the distance between iterates decreases and so the eigenvalue decomposition of consecutive iterates is not perturbed too much. Therefore, it is expected that few sweeps of the Jacobi method are required to compute the (approximate) eigenvalue decomposition of the next iterate.  The paper further provides bounds on the error of the prox operator as a function of the Frobenius norm of the matrix of off-diagonal entries of the iterate of the Jacobi method. This allows one to establish termination criteria in order to run the Jacobi method to obtain an approximate eigenvalue decomposition, while guaranteeing convergence of the (approximate) proximal algorithm.  The authors provide some numerical experiments on several matrix optimization problems such as solving sparse inverse covariance selection, SDPs, nuclear norm minimization, for small to moderately sized problems. Their experiments show that using the Jacobi algorithm for a small number of sweeps obtains sufficient accuracy to proceed to the next iteration of the proximal algorithm, and therefore reduces the overall wall-clock time to convergence, compared to using an eigenvalue/singular value decomposition implemented on CPUs or GPUs.

This paper proposes a new method to accelerate proximal-mapping based algorithms  where the proximal mapping usually needs computing the eigen/singular value decomposition.   The authors discusses how to use the Jacobi method to get the eigen/singular value decomposition approximately and parallelly,  and the convergence can be guaranteed by easily controlling the norm of the off diagonal elements.   Numerical experiments show the significant improvement over the state-of-art algorithms.

The paper considers proximal algorithms for convex optimization problems over matrices. When the function to be optimized involves a function of the eigenvalues or singular values of a matrix, then prox methods need to evaluate the proximal operator in each iteration which usually boils down to a full eigenvalue decomposition or SVD. This can be very costly. The paper proposes to use the Jacobi eigenvalue algorithm for computing such decompositions. The Jacobi eigenvalue decomposition can be warm started with the decomposition from the last iteration. On average, a good enough approximation to the eigenvalue decomposition can thus be computed by less than one sweep ($n/2$ parallel Givens rotation) over the matrix. Hence, the evaluation of the proximal operator in ADMM or prox algorithms can be speeded up significantly (by a factor of 5 to 20).  Convergence results are provided and experiments validate the approach and its efficacy.