Faster proximal algorithms for matrix optimization using Jacobi-based eigenvalue methods
Authors: Hamza Fawzi
Harry Goulbourne
Venue: NeurIPS 2021 Poster (2021)
Keywords: proximal splitting algorithms
matrix optimization
semidefinite programming
eigenvalue computation
Jacobi's method for eigenvalue decomposition

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.

Faster proximal algorithms for matrix optimization using Jacobi-based eigenvalue methods
Hamza Fawzi
Harry Goulbourne
Department of Applied Mathematics and Theoretical Physics University of Cambridge, UK {hf323, hmg42) .c.uk
Abstract
We consider proximal splitting algorithms for convex optimization problems over matrices. A significant computational bottleneck many these algorithms is the need to compute full eigenvalue or singular value decomposition at each iteration for the evaluation of proximal operator.
In this paper we propose to use an old and surprisingly simple method due to Jacobi compute these eigenvalue and singular value decompositions, and we demonstrate that it can lead to substantial gains in terms of computation time compared standard approaches. We rely on three essential properties of this method: (a) its ability to exploit an approximate decomposition which in the case of iterative optimization algorithms can be obtained from the previous iterate; (b) its parallel nature which makes it great fit for hardware accelerators such 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, 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 proximal operators obtained using approximate eigenvalue or singular value decompositions.
an initial point,
GPU, we can obtain to 10x speed-ups in
1
Introduction
Many problems in statistics and learning can be formulated as convex optimization over matrices. This includes kernel learning [LCB+044 low-rank matrix recovery [RFPIO CR09), sparse covariance selection [BGAN06] FHT08), community detection [BBV16], and many others. Proximal splitting algorithms such (accelerated) proximal gradient methods [Nes03] BTT9), the alternating direction method of multipliers (ADMM) [EB92 [PCH]] 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)+G(X) X
via calls to the proximal operators of F and/or G. When the variable X is symmetric or rectangular matrix, is often the case that the functions F' and/or G are symmetric functions the eigenvalues or singular values of X; i.e., so-called spectral or othugonally-invarinn functions. Typical examples are F(X) [(X > 0) = Ei (() 2 0) the convex indicator function for the positive semidefinite cone, or F(X)= = XX| - (X) the nuclear norm.
Computing the proximal operator of such spectral functions F requires a full eigenvalue/singular value decomposition which contributes large computational cost of least n² floating-poin
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

operations per iteration. The most common algorithms for such decompositions are the QR method, and the (LAPACK) that are used in standard numerical computing software such as MATLAB or NumPy.
algorithm [GL[33] which are implemented in linear algebra packages
Contributions In this paper we propose to use 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 eigenvectorsssingular 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 now very common in machine learning.
(iii) It has 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 cigenvaluc/singular value decomposition. When combined with existing results on the convergence of proximal splitting methods with inexact proximal oracles, we able 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 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 standard approaches on CPU and GPU.
Related work Several approaches have been considered to avoid computing a full eigenvalue or singular value decomposition each iteration, when the solution is known be low-rank. Popular approaches include the Lanczos algorithm, see c.g., (CCSIO), randomized SVD [HMTI]], 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 (KnyOT)) 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 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 CCV88 and [V21] 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. á is a strongly-convex function, the o- d-proximal operator of F is defined by
prox (1: a) argming (F(x) + (a, x) p(x) 2))
The paper [CVI8 considers optimization problems over symmetric positive semidefinite Toeplitz matrices and shows that by choosing & to the kernel for the so-called takura-Saito distance one can compute the generalized proximal operator in O(n²) time, instead of (n²)). In [V2T] 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 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 [KKMNOT) e.g., (ZFP+20 MKL15). We note that these contributions are specific to sparse the present paper because they still require an eigenvalue decomposition to be computed at each iteration.
semidefinite programs and are complementary
Finally let us mention that there are first-order methods that do not require full eigenvalue decompo- sition at each iteration but only the computation of the smallest cigennaane(s)/eiieneeetrr at each iteration such as c.g., (HR00 Ren 14 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 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/singuaar value decompositions yield guaranteed approximations of proximal operators for spectral functions. Finally we present numerical experiments in Section
2 Background
In this section we present some notations and background results concerning spectral functions. If ff : E + R is convex function on Euclidean space E' we let 0f(x) be the subdifferential of at z, namely (f()) = {g E E: f(v) 2 f(x) + {g. y - =)}. An element g E 0f(x) is called a subgradient of f at The proximal operator of f at defined as
prox = argmin (f(x) + (1/2)][ - y/(²}
We let S be the Euclidean space of n symmetric matrices with the trace inner product (X,Y) tr(XY). We let S (resp. +) be the cone of positive semidefinite (resp. definite) matrices in S". Also (n) denotes the set of n orthogonal matrices.
Spectral functions A function F S + RR of symmetric matrices is called spectral if there permutation-innarian function f RD I such that F(X) f(X(X) where X(X) (X)( A, (X)) ) is the vector of eigenvalues of X (in any order). It is clear that spectral functions orthogonally-invariin,, i.e., F(U XU) F F(X), because U XU and X have the same cigenvalues; conversely any orthogonally- -nnaariallyyinariiantt function F is the spectral function associated the permutation-invarine about spectral functions is a result due to Davis [Dav57]. which asserts that J' is convex as soon as f is, and gives a variational formula for F certifying this convexity.
function f(x) = (z)). The main remarkable fact that we will
Theorem (Dav557,, see also [BTNOT p.282, claim DJ). Let
R" + : invariant function, and let F S + R. F(X)= f(X()) be the associated spectral function. If is convex on R", then F is convex on S", and we have the following variational formulation for F:
be a permutation-
F(X)=
max [(diag(U* XU)).
UEO(n)
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 E S" have
prox; diag(prox (V))UT
where Y = U diag(y)U" an eigenvalue decomposition of Y.
Functions of singular values Similar results hold for functions F R + R of rectangular matrices that are orthogonally
i.e., FV) F(X) for all U 0(m) and V E Ó()). Any such function can be written FXX = (X)()) where oXX E gmin(m,.) is the vector of singular values and f [gmin(m.m) + R is an absolutely symmetric function, i.e., it invariant under signed permutations: ( (T)) = f (æ1...., .2) for 1), =1 and permutation of {1....,A k} where = min(m, n) [ew95]. this case too, it holds that F is convex as soon as f is, and we have:
max F(X)= UEO(m))
fddiag(U
(VT)) XV)).
(3)
The variational formulas 2 and are consequence of well-known results concerning majorization relations between the diagonal of matrix and its eigenvalues or singular values [MOA79] Chapter 9]. For completeness we have included 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-invarinn functions.
3

---

3 Jacobi-based methods for eigenvalue and singular value decompositions
In this section we present Jacobi-based methods to compute eigenvalue and singular value decompo- sitions, 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
3.1 The method
Jacobi eigenvalue methods compute the eigenvalue decomposition of real symmetric matrix by iteratively climinating off-diagonal elements using Givens rotations, More precisely, let A be real symmetric matrix that we want to diagonalize. The Jacobi method begins by selecting an off-ciagonal index (i,j), with i < j. such that |Aiil > 0, and constructs 2 x 2 rotation matrix Au Ag Agj The rotation matrix can be embedded (i,j) which is otherwise equal to the identity
[-sine 9 sin @
sin cos 6 that diagonalizes the 2 x 2 submatrix
in a larger n X π. Givens rotation matrix at indices matrix:
G. = G5 = cos
and
Gpq = 6p.9 for Torp' & (i.j} or & (i,1),
Gig sin 0, Gj = sin 0
where 8p.9 if p q and 0 otherwise. This Givens rotation matrix satisfies (GAG*),; The Jacobi algorithm then proceeds by updating A GAG 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)
Σ ISitiin
APy
(4)
The algorithm can thus be summarized as follows:
Algorithm 1: Jacobi method for diagonalizing symmetric matrix A
while of((A) A) 2 tolerance do choose an off-diagonal index (i. such that |Ais) > 0
construct a Givens matrix G such that (GAG*) )g = 0 update A +- GAGT
end
In the classical Jacobi method, the off-diagonal entry (i,j) selected at each iteration the one for which [A,y is largest. It can be easily shown that with this choice, the quantity off satisfies after one iteration
(GAGT) < (1 (/n) (A)
where N = n(n 1)/2. As such off converges linearly to O. It has been shown by Schönhage, and van Kempen [Sh664 VK66al 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 n matrix the cost of this step O(n²) flops. Cyclic Jacobi methods instead cycle through all off- 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] [W1622] [V666). In practice, the number of sweeps needed to diagonalize matrix of size n 1000 is usually no more than 15.
The cyclic Jacobi method can be easily parallelized. Indeed, instead of eliminating a single off off-diagonal entry at each iteration, one can eliminate in parallel [n/2] off-diagonal entries that correspond to distinct rows/columns. More precisely by selecting indices (i1,j1), ..., (im Jm)
Parallelization

---

are all distinct one can compute in parallel the m Givens rotations that are needed eliminate the entries can then be performed in n steps where at each step k = 1, ..... 1 one chooses indices (i [ [k])...... (mm Jm (B:) that are pairwise distinct such 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].
where m |n/2] such that
i A single sweep
Other eigenvalue algorithms In general, the Jacobi method is not competitive with more tradi- tional eigenvalue algorithms such as the QR method or the ing floating point operations only, the QR method is roughly equivalent to two sweeps of Jacobi [GVL.13 TB197 One important difference between the Jacobi method and QR/divide-and-conoeee 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 directly. One consequence of this is that the Jacobi method can have faster convergence if off (AA) is small, whereas this is not the case of other algorithms that start by a reduction into tridiagonal form. This ability to exploit approximate diagonalization of A is what makes the Jacobi method particularly interesting for iterative optimization algorithms.
method. Count-
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 SVD problems, instead of 2 x 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 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 [GVL3 Section 8.6.4] for more details and references.
3.2 Jacobi-based methods in proximal algorithms
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 to use the eigenvector matrix from the previous iteration as warm-start for the current iterate. More precisely, assume that iteration one needs compute the proximal operator p (Yk) where F : S" - R is
the spectral function associated to the permutation permutationiivvarian invariant function f R" + R. We proceed as follows:
Let Uk-1 be the orthogonal eigenvector matrix computed at iterate
Apply the Jacobi method U Y=Un-1 with tolerance eg. This returns an orthogonal matrix W such that
off ( U YUU-UW S c
Set U4 = W (eigenvector matrix of Y2) and yı = diag(U? (U Y&Ux) (eigenvalues of Y2) Output UR diag(prox ((x)) as the proximal operator prox (Yx).
As demonstrate in the numerical experiments section below, very small number of sweeps (often just one) are in general needed to converge to the eigenvalue decomposition of Y& in the second step above, because the matrix U. YAU is usually very close 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 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 Us not to be normalized after a large number of iterations. To resolve this issue one can orthonormalize the columns of Un periodically using a QR decomposition.
4 Approximate proximal evaluation
Most proximal algorithms (e-g., ADMM, proximal gradient method, ...) are still convergent even only approximation of the proximal operator is computed at each iteration, provided the approxi-
5

---

mation quality suitably decays to zero, see e.g., [B92 SRBIT Y BTB20] A natural question to understand how the tolerance parameter of the Jacobi method G₁ affects the approximation quality of the proximal operator prox f (V&)? This section addresses this question.
Measures of quality of prox There are multiple ways to quantify the approximation quality of proximal operator, see e.g., [BTB20] which compiles a list of different ways that were considered in the literature. this section we will focus on two simple notions that have been used convergence guarantees for proximal splitting methods. Let F: E R be convex function for which we want to compute the prox at a point Y E. Let
(x; F(X) 1 YI
be the prox objective so that prox xY) = argmin ((;;)). & E' is a candidate value for pY) the two measures of approximations we consider are: (i) Euclidean distance to the true prox, i.e., 8 || proxp (Y))| and (ii) suboptimality in prox objective value, i.e.,
? ({; (prrox x (Y)); Y). The former notion was used e.g., to guarantee convergence of ADMM algorithms, assuming that the sequence of errors (6x) is summable [EB92 Yáoló]; and the latter was used in [SRBII] guarantee convergence of (fast) proximal gradient methods, under decaying assumption on (77²)
Remark 2. Note that since & is 1-strongly convex we always have the inequality:
+ (x:; Y) ((prox 2 5 - - proxx
This means that any bound the approximation quality in terms of prox objective value, implies bound on the approximation quality using the Euclidean norm.
Approximate prox using approximate cigenvalue/singular decomposition Our main result gives bounds on the quality of the prox of any convex spectral function computed using an ap- proximate eigenvalue decomopsition. More precisely, assume that F S" + R is convex spectral function
function of symmetric matrices associated to the permutation-invarie Y S" for which want to compute pYo) We know that prox, = U prex, (y)UT where y (yi,..., Vo) is the vector of eigenvalues of Y, and U is the orthogonal matrix that diagonalizes Y. Assume that Ü is orthogonal matrix such that ofT(Ü=YÚ) S e, i.e.
Let
YTYū = diag(() + where
{
|| F S F4 = for i = 1,...,n
(6)
Our main theoretical result, stated below, quantifies the error in the approximation X proxp(()) where
x U diag(prox,
Theorem 2. Under the assumptions above, we have:
(i) || proxp(Y))l S and
(ii) (x;; Y) (prox, (Y); < 2/2.
A similar result can be stated for Theorem 3. Let F mRxx - RR. be convex orthogonall-invariin function asxociated with the function fx) F(diag(=:) defined Te m.. Let Y Rmxn and assume that Ü E (m)) and E (n) are orthogonal matrices such that
functions of rectangular matrices.
[||E] S
ÚTYV = =diag(0) ( + E where
Eu = 0 f min(m, n).
Then, with X = t diag(prox, (i)) VT, the following holds: (i) 1 ( S and
(ii) + 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 [GN22 The main result in [GNR20] shows that when F the indicator function the positive semidefinite cone (so that pr is
the Euclidean projection on 5 then || - prox, (Y))l & /2e. 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 [].
We believe Theorems
and are remarkable in number of ways:
They bound the distance to the true prox (a quantity that hard to compute) using simple quantity which is off YŪ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 Ú 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/eigenneeoor computation (e.g.. the Davis- Kahan theorem) suggest that the error would depend 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 c.
Implications If we combine our theorems with [EB92] and [RBI1] 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 F(X) G(X)
X
where F(X) is convex spectral function and G(X) is convex and smooth with L-Lipschitz gradient. Consider the (accelerated) proximal gradient method
Xx+1 = (2 - (1/L)GG()))
where Z Xx for the standard proximal gradient method, and Z Xx
(X X
E+2
for the accelerated version. Assume that the proximal operator in (9) is evaluated approximately according to ( with accuracy Cg- Then, using [RBII] Props. and 2] and Theorem Lii). we immediately get: (a) if (cx) is summable (e.g., C₂ = ((1/k++)) then the (basic) proximal gradient method converges at a rate O(1/k) to the optimal value of (8) and (b) if (kea) is summable (e.g.,
< = 01/k22+) then the accelerated proximal gradient method converges at a rate (114k²) to the optimal value of For the ADMM algorithm, we can combine the result of EB92 Theorem 8] with Theorem i) to guarantee convergence when (fx) is summable.
Numerical experiments
In this section present numerical experiments with facobi-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 (NV120) 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 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 i used is NVIDIA TITAN Xp with CUDA 11.0, and the CPU is 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) use 32-bit floats, instead of 64- doubles, which are much faster GPU. On the examples we have tried this had little or no effect on the quality of the solutions, especially because are looking problems where moderate accuracy is satisfactory.
+ 17] functions linalg eigh and linalg svd. The GPU we
We note that the proof of inequality (i) does actually use the assumption that En 0 for 1,...
7

---

Performance of Jacobi's method on GPU Before showing the results on optimization problems, we first compare in Table [ the time needed for a single Jacobi sweep vs. the time needed for full eigenvalue (resp. singular value) decomposition using NumPy/CuPy's linalg. eigh (resp. linalg. svd), on 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 Z 2000.
Symmetric eigenvalue decomposition
Singular value decomposition
default (CPU)
default (GPU)
jacobi-Isw (GPU)
default (CPU)
default (GPU)
jacobi-Isw (GPU)
n
128
2.18
5.91
0.83
3.75
12.23
1.14
256
6.86
12.72
1.43
12.93
39.93
1.85
512
26.24
31.78
3.29
62.46
107.47
3.99
1024
152.44
63.28
11.45
385.59
476.22
15.17
2048
1180.78
163.62
59.94
3,453.76
3,223.60
77.18
4096
8491.13
532.92
682.72
24,869.11
19,649.37
500.99
Table I: Times, in milliseconds, for one sweep of the Jacobi method on a GPU vs. full eigenvalue/Singullr-vvlle decomposition using the default method CPU/GPU.
5.1 Sparse inverse covariance selection
We consider here the well-known statistical problem of estimating unknown covariance matrix I E S from samples ",,..., IN E IR" of Gaussian distribution (0, M), under the assumption that Σ is sparse. The approach proposed in [BGdN06 FHT08 consists of solving the following convex optimization problem:
min XEST
tr(šX) - llog det X + X|X|||
(10)
where A > 0 is a regularization parameter, S E ST. is the given empirical covariance matrix, and ||||| = Eig X,yl is the entrywise & -norm of X. The ADMM algorithm one of the most popular methods for this problem, see e.g., (SMGIO BPCIT] and references therein. The problem (T0) can
be naturally put in the splitting form (I) with F(X) tr($X) log det X and G(X) x| l| The ADMM iterates can be written as:
X = prox plog dot (Zp= - Un - ps)
Zk+1 = -Sx/o(X+t] + Us)
(11)
Uk+1 = UK + Xx+1 Zx+1
where the proximal operator of -p det is prox-plog der (Y) = V diag (I / + 44)/22) vT where Y V VT is the eigenvalue decomposition of Y; and St is the entrywise soft- thresholding operation defined by (x = 0 if |x|| S and x - sign(x)t otherwise. The main computational cost in the iterations (11) is the proximal update rule for X which requires the eigenvalue decomposition of Zμ UR S 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.22
We have implemented this approach with different values of n in the range 200-1000, where the covariance matrix Σ S was chosen be banded matrix with constant bandwidth A 5. We have used the parameters and stopping criterion of the Matlab implementation of (TT) from [BCTT] available htttpp:://eb.statttetee ord. werdddddduuddppdd papers//p////////p///e admm/
We compare the Jacobi method
It is not clear to us why the performance of the Jacobi- Jacobi-symmetric symmetric eigenvalue routine syevj significantly degrades for sizes > 2000. Note that, surprisingly, the Jacobi-SV 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/k², more precisely we set Eş | || r1(1+ k²), where Yú, is the matrix to diagonalize at iteration k. We see that the average number of sweeps per iteration is ze confirming the fact that the eigenvector matrix from iterate k 1 a good initial point for the k'th iterate.
iterations
X X-upate time per iter (ms)
default (CPU)
default (GPU)
jacobi (GPU)
default (CPU)
default (GPU)
jacobi (GPU)
sweeps/iter
200
308
308
309
6.20
7.43
1.81
0.97
400
269
269
270
25.83
19.56
3.06
0.97
600
248
248
249
62.64
33.41
5.62
0.98
800
236
236
237
100.06
46.55
10.7
0.99
1000
227
227
227
234.32
67.42
16.4
0.99
Table 2: Sparse covariance selection problem. Time per iteration for the computation of the X-
proximal update in [II) using the Jacobi method compared the standard approaches on CPU/GPU.
5.2 Semidefinite programming
We consider here semidefinite programs in primal-dual form,
min zER" * s.t.
x
max
-tr(BY)
YES'
+ + S
S.L.
Yz (),Y) -q (i ,.....n) n)
(12)
where c E R", in (CCPB16) 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 S The projection step be computed using the Jacobi- 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,
E and B E Sh primal-dual ADMM-based algorithm was proposed
||&al= f k2),, where Ys is the matrix to diagonalize iteration k. We see that the time spent 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 1.
iterations
proj time/iter (ms)
linsolve/iter (ms)
default (CPU)
default (GPU)
jacobi (GPU)
default (CPU)
default (GPU)
jacobi (GPU)
sweeps
k:
per iter
200
60
60
60
5.28
8.01
11.63
1.00
2.46
400
80
80
80
21.75
20.57
3.70
1.01
10.07
600
90
90
66.05
32.62
7.97
1.01
25.76
800
90
90
90
95.20
47.9
11.65
1.01
39.10
1000
150
150
150
247.76
74.24
18.55
1.01
58.76
Table 3: Semidefinite programming using primal-dual ADMM approach of [OCPB16] on random Ig. and the A; are sparse random matrices with density k. Time per iteration (in milliseconds) to compute projection of & x l symmetric matrix on the positive semidefinite cone using the Jacobi method vs. standard approaches on CPU/GPU. For comparison, we also show the time needed for the linear system solve at each iteration (in milliseconds). We used the same termination criterion as the one in (OCPB16), with a tolerance of 10-3 on the normalized primal and dual residuals, and duality gap.
instances of (12) where n = 3k:

---

5.3 Nuclear norm minimization using the proximal gradient method
consider the low-rank matrix recovery problem
min XER*Xm
(1/2)]lA(X) - -b1 + |||X||~
(13)
where A: R -> IRm E IR", is the regularization parameter, and = (x is the nuclear norm. A popular method for solving (3) is the (fast) proximal gradient method, with iterates:
X = SVT -A b))
(14)
where >> 0 is the step size, and SVT,x 1 is the proximal operator associated to the nuclear norm which given by the soft-thresholding of singular values. the standard proximal gradient method, we take Z X2 in (14), whereas for the fast proximal gradient method we take Z Xx + / (Xx - X
Clearly, the main computational cost in the implementation of [14 decomposition of n X n matrix.
compute the singular value
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 the proximal operator, on random matrix-completion instances of [3). Like for the previous examples, used relative tolerance for the Jacobi method decaying like 1/k², except for the largest instance where we fixed the number of sweeps to be at each iteration. The number of iterations (14 was fixed and set to 1000, with time
limit of hour on the total running time. For the largest instance the two default methods exceeded
the time limit and the objective value report is the one obtained after 1000 iterations.
objective
SVD time per iter (ms)
default CPU
default GPU
jacobi GPU
default CPU
default GPU
jacobi GPU
jacobi sweeps per iter
n
200
7.64c+01
7.64c+01
7.64c+01
7.1
16.8
2.1
1.4
400
3.15e+02
3.15e+02
3.15e+02
29.4
53.1
4.9
1.9
600
7.24++02
7.24e+02
7.24e+02
71.1
105.4
10.1
2.3
800
1.29e+03
1.29e+03
1.29e+03
149.5
177.7
17.4
2.5
1000
1.99e+03
1,99c+03
1.99c+03
306.2
327.7
29.3
2.6
4000
4.65e+04*
3.94e+04"
3.41c+04
22,641.5
10,676.7
597.5
1.0
Table Nuclear norm minimization using fast proximal gradient method (14). Time to compute the SVD at each iteration using the Jacobi's method vs. the default methods on CPU/GPU. Objective value is after 1000 iterations (except for 4000, default CPU/GPU).
'Ojeectvve value after the time limit of 1h was reached.
"The number of sweeps of the Jacobi method was fixed to for the instance n 4000.
Acknowledgments
H.G. supported by GlaxoSmithKline PhD studentship. H.F. would like thank Oisín Faust and Alhussein Fawzi for comments on an earlier draft of this paper.
References
[BBV16]
Afonso S Bandeira, Nicolas Boumal, and Vladislav Voroninski. On the low-rank approach for semidefinite programs arising synchronization and community detection. Conference on learning theory, pages 361-382, PMLR, 2016.
[BGdN06]
Onureena Banerjee, Laurent El Ghaoui, Alexandre A'spremont, and Georges Natsoulis. Convex optimization techniques for fitting sparse gaussian graphical models. In Proceedings of the 23rd international conference on Machine learning, pages 89 96, 2006.
[BM03]
Samuel Burer and Renato DC Monteiro. A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization. Mathematical Programming, 95(2):329-33 2003.
10

---

[BPC11]
Stephen Boyd. Neal Parikh, and Eric Chu. Distributed optimization and statistical learning via the alternating direction method multipliers. Now Publishers Inc, 2011.
[BT09]
Amir Beck and Marc Tehoulle. A fast iterative shrinkage-threshollinng algorithm for linear inverse problems. SIAM Journal on Imaging Sciences, 2(1):183-202, 2009.
[BTB20]
Mathieu Barré, Adrien Taylor, and Francis Bach. Principled analyses design of first-order methods with inexact proximal operators. arXiv preprint arXiv:2006.06041 2020.
[BTNOI]
Aharon Ben-Tal Tal and Arkadi Nemirovski. Lectures on modern conver optimization: analysis,
algorithms, and engineering applications. SIAM, 2001.
[CCS10]
Jian-Feng Cai, Emmanuel J Candès, and Zuowei Shen. A singular value thresholding algorithm for matrix completion. SIAM Journal on optimization,
2010.
[CR09]
Emmanuel J Candès and Benjamin Recht. Exact matrix completion via convex optimization. Foundations of Computational mathematicy, 96):717-772, 2009.
[CV18]
Hsiao-Has Chao and Lieven Vandenberghe. Entropic proximal operators for nonnegative trigonometric polynomials. JEEE Transactions on Signal Processing, 66 66(18):4826-4838 2018.
[Dav57]
Chandler Davis. All convex invariant functions of hermitian matrices. Archiv der Mathematik, 8(4):276-278, 1957.
[FB92]
Jonathan Bekstein and Dimitri P Bertsekas. On the Rachlos-Rachfrrd splitting method and the proximal point algorithm for maximal monotone operators. Mathematical Programming, 55()):2933318 1992.
[FH60]
George Elmer Forsythe and Peter Henrici. The cyclic Jacobi method for computing the principal values of a complex matrix. Transactions of the American Mathematical Society, 94(1):1-23, 1960.
[FHTO8]
Jerome Friedman, Trevor Hastie, and Robert Tibshirani. Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3):432-44) 2412,4411101 2008.
[FKMNOI]
Mituhiro Fukuda, Masakazu Kojima, Kazuo Murota, and Kazuhide Nakata. Exploiting sparsity in semidefinite programming via matrix completion i: General framework. SIAM Journal on Optimization, 113):647-674 2001.
[GNR20]
Paul Goulart, Yuji Nakatsukasa, and Nikitas Rontsis Accuracy of approximate projection to the semidefinite cone. Linear Algebra and its Applications, 594:177-192, 2020.
[GVI.13]
Gene H Golub and Charles Van Loan, Matrix computations, volume 3. JHU press, 2013.
[HMT11]
Nathan Halko, Per-Gunnar Martinsson, and Joel A Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions, SIAM review, 53(2):217-288, 2011.
[HMvdW+ 20] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel Smith, Robert Kern, Matti Picus. Stephan Hoyer. Marten H. van Kerkwijk. Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre
Kevin Sheppard, Tyler Reddy. Warren Weckesser, Hameer Abbasi. Christoph Gohlke, and Travis E. Oliphant. Array programming with NumPy. Nature, 585(7825):357-362 September 2020.
[HR00]
Christoph Helmberg and Franz Rendl. A spectral bundle method for semidefinite programming. SIAM Journal on Optimization, 10(3):633-60 696, 2000.
[Jac46]
Carl Gustav Jacob Jacobi. Über leichtes Verfahren die in der Theorie der Säcularstörungen vorkommenden Gleichungen numerisch aufzulösen. J. Reine Angew. Math, 30:51-95, 1846.
[JV21]
Xin Jiang and Lieven Vandenberghe. Bregman primal-dual first-order method and application to sparse semidefinite programming. 2021.
[Kny01]
Andrew V Knyazev. Toward the optimal preconditioned eigensolver: Locally optimal block preconditioned conjugate gradient method. SIAM journal scientific computing, 23(2)5517 541, 2001.
[LCB+04]
Gert RG Lanckriet, Nello Cristianini, Peter Bartlett, Laurent El Ghaoui, and Michael I Jordan. Learning the kernel matrix with semidefinite programming. Journal of Machine learning research, 57727): 2004.
[Lew95]
Adrian S Lewis. The convex analysis of unitarily invariant matrix functions. Journal of Convex Analysis, 2(1):173-18)) 1995.
[MKL15]
Ramtin Madani, Abdulrahman Kalbat, and Javad Lavaei. ADMM for sparse semidefinite programming with applications to optimal power flow problem. In 2015 54th IEEE Conference Decision and Control (CDC), pages 5932-5933 5339. IEEE, 2015.
11

---

[MOA79]
Albert W Marshall, Ingram Olkin, and Barry C Arold. Inequalities: theory of mujorization and its applications. volume 143. Springer, 1979.
[Nes03]
Yurii Nesterov. Introductory lectures on convex optimization: basic course, volume 87. Springer Science & Business Media, 2003.
[NV120]
NVIDIA. cusolver library. 2020. httpe::/ddocs.nvidda.. wrw.oocss.iiiiii.....ooooo wwwpaiii..o/////////SS/......... pdf
[OCPB16]
B. O'Donoguue, E. Chu, N. Parikh, and S. Boyd. Conic optimization via operator splitting and homogeneous self-dual embedding. Journal of Optimization Theory and Applications,
169(3):1042-1000 10680220888, June 2016.
[OCPB19]
B. O'Donoghue, E. Chu, N. Parikh, and S. Boyd. SCS: Splitting conic solver, version 2.1.3.
sttpe:///gttaa...coo//ccrrssssss November 2019.
[OUN+17] 17]
Ryosuke Okuta, Yuya Unno, Daisuke Nishino, Shohci Hido, and Crissman Loomis. CuPy: A HumPyy-Coppatibe Library for NVIDIA GPU Calculations. In Proceedings of Workshop Machine Learning Systems (LearningSys) in The Thirty-first Annual Conference on Neural Information Processing Systems (NIPS), 2017.
[PB14]
Neal Parikh and Stephen Boyd. Proximal algorithms. Foundations and Trends in optimization, 1(3):127 239, 2014.
[Rea14]
James Renegar. Efficient first-order methods for linear programming and semidefinite program- ming. arXiv preprint 1.0Xvv140998822 2014.
[RFP10]
Benjamin Recht, Maryam Fazel, and Pablo A Parrilo. Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization. SIAM review, 52(3):471-555 1501, 2010.
[RGN19]
Nikitas Rontsis. Paul J Goulart, and Yuji Nakatsukasa. Efficient semidefinite programming with approximate admm. arXiv preprint arXiv:1912.02760 2019.
[Sch64]
A Schönhage. Zur quadratischen Konvergenz des 1964.
Numerische Mathematik,
[SMG10]
Katya Scheinberg. Shiqian Ma, and Donald Goldfarb. Sparse inverse covariance selection via alternating linearization methods. arXiv preprint 1.0097 2010.
[]RB11]
Mark Schmidt, Nicolas Ie Roux, and Francis Bach. Convergence rates of inexact proximal- gradient methods for convex optimization. Neural Information Processing Systems, 2011.
[TB197]
Lloyd N Trefethen and David Bau III. Numerical linear algebra, volume 50. SIAM, 1997.
[VK66a]
HPM Van Kempen. On the convergence matrices with non-distinct eigenvalues. Numerische Mathematik, 9(1):11-18 1966.
classical jacobi method for real symmetric
[VK66h]
HPM Van Kempen. On the quadratic convergence of the special cyclic Jacohi method. Nu- merische Mathematik, 9(1::19-22, 1966.
[Wil62]
J. H. Wilkinson. Note the Quadratic Convergence of the Cyclic Jacobi Process. Numerische Mathematik 296- 300, 1962.
[Yao16]
Wang Yao. Approximate versions of the alternating direction method of multipliers. PhD thesis, Rutgers University-Graduate Graduate School-New Brunswick, 2016.
[YTF*21]
Alp Yurtsever, Joel Tropp. Olivier Fercoq, Madeleine Udell, and Volkan Cevher. Sealable semidefinite programming. SIAM Journal on Mathematics of Data Science. 3(1):171-200, 2021.
[ZFP*20] 20)
Yang Zheng, Giovanni Fantuzzi, Antonis Papachristodoulou, Paul Goulart, and Andrew Wynn. Chordal decomposition in operator-pllitingg methods for sparse semidefinite programs, Mathe- matical Programming, 180(1):489-533 2020.
12


Reviews:
[{"review_id": "xgM9t74WYLp", "reviewer": "Reviewer_gv2k", "summary": "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.  \nNumerical experiments show the significant improvement over the state-of-art algorithms.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This paper is well-organized and well-written in general.  I have the following main comments:\n1. The property of the spectral function in Theorem 1 will play a key role in proving Theorem 2. I suggest clarifying it in section  2.\n\n2. The assumptions on $F(X)$ should be put in  the assumption environment so that the statement of Theorem 2 becomes more clear. \n\n3. It says 'One consequence of this is that the Jacobi method can have faster convergence if off(A) is small, whereas\n  this is not the case of other algorithms that start by a reduction into tridiagonal form.' I'm not sure about this. Could you give more details in why the warm-start does not work in other eigenvalue algorithms? For example, can this be verified by experiments? The current implementation in NumPy solves the eigen decomposition exactly.\n\nMinor comment:\nThe stopping criterion of in Section 5.1 and 5.3 should be given. \n\n\n ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Not apply.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "dPxG6q05Ia0", "reviewer": "Reviewer_KM7u", "summary": "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).\n\nProximal 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.\n\nThe 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.\n\nThe 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.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Originality:\nI've not seen the Jacobi algorithm applied in this way, so as far as I can tell, the result is original. Though most results in this paper are not particularly new, there is merit in combining the ideas of using a few Jacobi sweeps (with guaranteed and easy to compute convergence criteria).\n\nQuality:\nOverall, the results of the paper appear sound, though I have some comments about certain results and experiments.\n  - Beginning of section 5: I do not see the significance of comparing a single sweep of the Jacobi algorithm against the computation of the full eigenvalue/singular value decomposition. A single sweep will not compute anything meaningful and obviously does less work than any full decomposition, so Table 1 is not helpful and in general misleading\n  - Numerical experiments in general: these problems are quite small by today's standards, and the iteration time per problem is quite fast (on the order of ms) for all methods involved. Especially given the message of the paper, which focuses on exploiting parallelism for large scale problems, it is necessary to increase the size of the problems in the experiments section.\n  - Remark 1: what floating point precision was used in your experiments? It is clear that continuously applying givens rotations would result in loss of orthogonality over time, but this can be addressed by re-orthogonalizing U periodically (which s equivalent to 1 sweep of Jacobi; in fact, it should be faster since it can be implemented using BLAS 3 operations on the GPU).\n\nSome other comments:\n  - It would be helpful to include convergence results for the cyclic Jacobi methods. The authors provide convergence results for the classical method that picks the largest off-diagonal entry, but this algorithm isn't used in this paper.\n  - Out of curiosity, are there convergence results for a randomized Jacobi algorithm.\n  - For each experiment, what was the fraction of time spent running the Jacobi sweep vs. computing the Frobenius norm of the error? In fact, how often was the Frobenius error computed to determine whether sufficient accuracy was reached?\n\nClarity:\nThe paper is overall well-written, and I think communicates its ideas well.\n\nSignificance:\nIf this paper can be shown to work well on truly large scale problems, then it's significance would be improved. However, with the current experiments, it's difficult to say that it's a truly significant result.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "As discussed in the previous section, my main concern is over the numerical experiments, and seeing this applied on truly large scale problems, rather than those that can be solved with conventional approaches quickly already.\n\nI do not see any issues regarding this work and potential negative impacts on society.", "review_text": "", "extra_scores": {"ethical_concerns": "N/A.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "c6_mf9qMgZ", "reviewer": "Reviewer_6wmm", "summary": "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).\n\nConvergence results are provided and experiments validate the approach and its efficacy. ", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is very well written and a pleasure to read. It provides a very nice introduction to the Jacobi eigenvalue algorithm and uses it within the framework of (approximate) proximal methods. I did not check the proofs but the theorems seem correct. The core contribution of the paper is to bring Jacobi's eigenvalue algorithm \n\nSome concerns though: \n1. The paper uses the resulting algorithm for solving three problems in the experiments: (i) sparse inverse covariance selection, (ii) semidefinite programming, and (iii) low-rank matrix recovery. All experiments are run on synthetic data. There are better algorithms to solve these problems. For (i), use [1] or [2]. For (ii), make use of interior point methods for problems of this scale since they considerably need fewer iterations. For (iii), use a Frank-Wolfe type algorithm. These are the standard methods to solve each of these problems and they are considerably faster than the ADMM approach presented here. While the presented approach is nice and general, one would not use it in order to solve these problems on real-world machine learning instances. Hence,  I do not know, whether NeurIPS is the right venue for such a paper. Maybe some optimization journal/conference would be better suited.\n\n2. The improvement in speed-up for each eigenvalue / SVD decomposition is provided. It would be nice to also report the total solving time, since this speedup does not translate to an overall speedup, e.g., in the case of the SDP experiment, it seems a factor of two improvement in the end.\n\n3. Comparison to state-of-the-art approaches on real-world problems/data would be nice for an ML conference.\n\n4. It is stated in the introduction that there are first-order methods that do not require solving a full eigenvalue decomposition / SVD but instead only need one eigenvector (Frank-Wolfe type algos for instance). It would be nice to also compare to such algorithms in the experimental section, since they are usually much more efficient and to provide a better global picture. \n\n5. What about interior point algorithms? They need much fewer iterations and since a full eigenvalue decomposition is computed in each iteration, the per-iteration complexity is comparable. Especially the size of the problems (up to 1000 x 1000)in the evaluation suggests using IP methods instead.\n\nSo all in all, the paper is really nice, but I am not sure if it is a suitable contribution for the ML community. \n\nReferences:\n[1] Cho-Jui Hsieh, Matyas A. Sustik, Inderjit S. Dhillon, Pradeep Ravikumar. Sparse Inverse Covariance Matrix Estimation Using Quadratic Approximation. In NIPS 2011.\n[2] Cho-Jui Hsieh, Matyas A. Sustik, Inderjit S. Dhillon, Pradeep K. Ravikumar, Russell Poldrack. BIG & QUIC: Sparse Inverse Covariance Estimation for a Million Variables. In NIPS 2013.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "N/A.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3 hours"}, "extra_text": {}}, {"review_id": "7jxwtb5Kx2t", "reviewer": "Reviewer_BV6g", "summary": "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.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The Jacobi-based approximate eigenvalue decomposition method is not novel, and its cyclic and parallelized implementation is also not new. These are introduced in section 3.1. There is certain novelty in applying the method to compute proximal operators in proximal algorithms, which is introduced in section 3.2. In order to better clarify the contribution, it is suggested that section 3.1 and 3.2 should be separated into two different sections, since the former is a review of existing work, and the latter is the main contribution.\n\nThe significance and novelty can be improved by bringing more theoretical insights. For example, is the number of sweeps bounded by the increment between the previous and current iterates? How many sweeps would be required for a certain class of proximal algorithms such as ADMM?\n\nRandomized approximate eigenvalue decomposition (or singular value decomposition) has been an active research topic. How is the adopted Jacobi-based method compared with these randomized methods? It is suggested that such discussion should be added in the related work, and in the experiments the proposed approach should be compared with these methods. \n\nFrom the experimental results in Table 1, one sweep of the Jacobi method is about 10 times faster when n is moderate (256 or 512), but the improvement is only about 3 times when n is large (2048). Does it indicate that the advantage would become less and less when n scales up?\n\nIn the experiments, there are some parameters in the accuracy $\\epsilon_k$, namely, the scale taking values 0.1 and 1 and the power taking value 2.5. Are these parameters carefully tuned? Are their choices crucial to the performance?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "It would be better to provide certain bound on the complexity of the method and to compare it with existing randomized approximate eigenvalue decomposition methods, in order to better demonstrate its advantage and limitation.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "8 hours"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The submission studies the application of the Jacobi method to compute approximate eigendecompositions or SVD in the context of solving composite optimization problems over matrices, where the objective functions are spectral. The main theoretical contribution of the paper shows how the error in the approximate decomposition (in terms of the $\ell_2$-norm of the off-diagonal entries) propagates into the analysis of proximal gradient descent and proximal accelerated gradient descent. The AC agrees with the reviewers that this is an interesting result, albeit fairly simple given previous work of inexact oracles for proximal methods. Moreover, the reviewers and the AC all find that the paper is well-written, making the technical contribution easy to grasp.

The main practical contribution is an experimental evaluation of the resulting Jacobi-based method against the  analogue QR-based methods, which generally shows the superiority of the Jacobi method. This is partly explained by the fact that the Jacobi method can effectively use previous decompositions as warm starts for the next. The AC agrees with reviewers BV6g and KM7u that the significance of such experiments would be increased by including larger matrices and real-world datasets. 

In conclusion, the AC believes that the interesting and well-explained theoretical contribution, together with the positive results on the experimental part, outweigh the lack of larger experiments and recommends that the paper be accepted.
Decision: Accept (Poster)