An SDP-Based Algorithm for Linear-Sized Spectral Sparsification
Yin Tat Lee He Sun







          






































 
 


 






 








fminipage
  
  

algbox[0]0.2in
6.3in


0.2in













































































































 










A


𝔼
𝐏
𝐃
0

𝐗























argmax
argmin
mod







 

























thmTheorem[section]  
theorem[thm]Theorem
fact[thm]Fact
lem[thm]Lemma
lemma[thm]Lemma
obs[thm]Observation
clm[thm]Claim
cor[thm]Corollary
fac[thm]Fact
rem[thm]Remark
remark[thm]Remark
*rem*Remark

pro[thm]Proposition
con[thm]Conjecture

constructionConstruction
protocolProtocol
defi[thm]Definition
assumption[thm]Assumption
probProblem





equationsection




















alpha
















An SDP-Based Algorithm for

 Linear-Sized Spectral Sparsification
    Yin Tat Lee

Microsoft Research
 Redmond, USA


He Sun

The University of Bristol

Bristol, UK



    
=====================================================================================================






For any undirected and weighted graph G=(V,E,w) with n vertices and m edges, we call a sparse subgraph H of G, with proper reweighting of the edges,  a (1+ε)-spectral sparsifier if

    (1-)x^L_Gx≤ x^ L_H x≤ (1+ε) x^ L_Gx

holds for any x∈^n,
where L_G and L_H are the respective Laplacian matrices of G and H. Noticing that Ω(m) time is needed for any algorithm to construct a spectral sparsifier and a spectral sparsifier of G  requires  Ω(n) edges, a natural question is to investigate, for any constant , if a  (1+)-spectral sparsifier of G with O(n) edges can be constructed in Õ(m) time, where the Õ notation suppresses polylogarithmic factors. 

All previous constructions on spectral sparsification <cit.> require either super-linear number of edges or m^1+Ω(1) time.



  
  
  In this work we answer this question affirmatively by presenting an algorithm that, for 
  any undirected graph G and >0, outputs a (1+)-spectral sparsifier of G with O(n/^2) edges in Õ(m/^O(1)) time. 

Our algorithm is based on three novel techniques: (1) a new potential function which is much easier to compute yet has similar guarantees as the potential functions used in previous references; 
(2) an efficient reduction from a two-sided spectral sparsifier to a one-sided spectral sparsifier; (3) constructing a one-sided spectral sparsifier by a semi-definite program.  
   



 Keywords:  spectral graph theory, spectral sparsification







empty










§ INTRODUCTION






A sparse graph is one whose number of edges is reasonably viewed as being proportional to the number of vertices. 
Since most algorithms run faster on sparse instances of graphs and it is more space-efficient to store sparse graphs, 
it is useful to obtain a sparse representation H of G so that certain properties between G and H are preserved, see Figure <ref> for an illustration.

 Over the past three decades, different notions of graph sparsification have been proposed and widely used to design approximation algorithms.  For instance, a spanner H
 of a graph G is a subgraph of G so that 
the shortest path distance between any pair of vertices is approximately preserved <cit.>.

Benczúr and Karger <cit.> defined a cut sparsifier of a graph G to be a sparse subgraph H such that the  value of any cut between G and H are approximately the same. 
In particular, Spielman and Teng <cit.> introduced a spectral sparsifer, which is a sparse subgraph H of an undirected graph G such that many spectral properties of the Laplacian matrices between G and H are approximately preserved.   
Formally, for any undirected  graph G with n vertices and m edges, we call a subgraph H of G, with proper reweighting of the edges,  a (1+ε)-spectral sparsifier if 

    (1-)x^L_Gx≤ x^ L_H x≤ (1+ε) x^ L_Gx

holds for any x∈^n,
where L_G and L_H are the respective Laplacian matrices of G and H. 
Spectral sparsification has been proven  to be a remarkably useful tool in algorithm design, linear algebra, combinatorial optimisation, machine learning, and network analysis.







In the seminal work on spectral sparsification, Spielman and Teng <cit.> showed that, for any undirected graph G of n vertices, a spectral sparsifier of G with only  O(nlog^cn/ε^2) edges 
exists and can be constructed in nearly-linear time[We say a graph algorithm runs in nearly-linear time if the algorithm runs in O(m·polylog n) time, where m and n are the number of edges and vertices of the input graph.], where c≥ 2 is some constant. Both the runtime of their algorithm  and the number of edges in the output graph involve large poly-logarithmic factors, and this motivates a sequence of  simpler and faster constructions of spectral sparsifiers with fewer edges <cit.>.  In particular, since any constant-degree expander graph of O(n) edges is a spectral sparsifier of an n-vertex complete graph, a natural question is to study, for any n-vertex undirected  graph G and constant >0, if a (1+)-spectral sparsifier of G with O(n) edges can be constructed in nearly-linear time. Being considered as one of the most important open question about spectral sparsification by Batson et al. <cit.>, there has been many efforts  for fast constructions of linear-sized spectral sparsifiers, e.g. <cit.>, however the original problem posed in <cit.> has remained open.

In this work we answer this question affirmatively by 
presenting the first nearly-linear time algorithm for constructing a linear-sized spectral sparsifier. The formal description of our result is as follows:





 Let G be any undirected graph with n vertices and m edges. For any 0<ε<1, there is an algorithm that runs in Õ(m/^O(1)) work, Õ(1/^O(1)) depth, and produces a (1+)-spectral sparsifier of G with O(n/^2) edges[Here, the notation Õ(·) hides a factor of log^cn for some positive constant c.]. 


maingraph  shows that  a linear-sized spectral sparsifier can be constructed in nearly-linear time in a single machine setting, and in polylogarithmic time in a parallel setting. The same algorithm can be applied to the matrix setting, whose result is summarised as follows:

 
Given a set of m   matrices {M_i}_i=1^m, where M_i∈^n× n. Let M = ∑_i=1^m M_i and Z=∑_i=1^mnnz(M_i), where nnz(M_i) is the number of non-zero entries in M_i.
For any 1>ε>0, there is an algorithm that runs in Õ((Z + n^ω)/^O(1)) work, Õ(1/^O(1)) depth and produces a (1+)-spectral sparsifier of M with O(n/^2) components, i.e., there is an non-negative coefficients {c_i}_i=1^m such that
| { c_i | c_i≠ 0}| =O(n/^2), and

    (1-)· M≼∑_i=1^m c_i M_i ≼
    (1+)· M.

Here  ω is the matrix multiplication constant.



 §.§ Related work





In the seminal paper on spectral sparsification, Spielman and Teng <cit.>  showed that a spectral sparsifier of any undirected graph G can be constructed by decomposing G into multiple nearly expander graphs, and sparsifying each subgraph individually. This method leads to the first nearly-linear time algorithm for constructing a spectral sparsifier with O(nlog^cn/ε^2) edges for some c≥ 2. 
However, both the runtime of their algorithm and the number of edges in the output graph involve large poly-logarithmic factors.  Spielman and Srivastava <cit.> showed that a (1+)-spectral sparsifier of G  with O(nlog n/^2) edges can be constructed by sampling the edges of G with probability proportional to their effective resistances, which is conceptually much simpler than the algorithm presented in <cit.>. 


 Noticing that any constant-degree expander graph of O(n) edges is a spectral sparsifier of an n-vertex complete graph, Spielman and Srivastava <cit.> asked 
if any n-vertex graph has a spectral sparsifier with O(n) edges.  To answer this question, Batson, Spielman and Srivastava <cit.> presented a polynomial-time algorithm that, for any undirected graph G of n vertices, produces a spectral sparsifier of G with O(n) edges.  At a high level, their algorithm, a.k.a. the , proceeds for O(n) iterations, and in each iteration one edge is chosen deterministically to “optimise” the change of some potential function. Allen-Zhu et al. <cit.> noticed that a less “optimal" edge, based on a different potential function, can be found in almost-linear time and this leads to an almost-quadratic time algorithm. Generalising their techniques, Lee and Sun <cit.> showed that a linear-sized spectral sparsifier can be constructed in time O(m^1+c) for an arbitrary small constant c.

All  of these algorithms proceed for Ω(n^c) iterations, and
every iteration takes  Ω(m^1+c) time for some constant c>0.  Hence, to break the Ω(m^1+c) runtime barrier faced in all previous constructions multiple new techniques are needed.












 §.§ Organisation


The remaining part of the paper is organised as follows. We  introduce necessary notions about matrices and graphs in  pre. In overview we  overview our algorithm and proof techniques.  For readability, more detailed discussions and technical proofs 
 are deferred to details.


§ PRELIMINARIES







 §.§ Matrices


For any n× n real and symmetric matrix A, let λ_min(A)=λ_1(A)≤⋯≤λ_n(A)=λ_max(A) be the eigenvalues of A,  where λ_min(A) and λ_max(A) represent the minimum and maximum eigenvalues of A. 
We call a matrix A  positive semi-definite () if  x^Ax≥ 0 holds for any x∈ℝ^n, and a matrix A positive definite  if  x^Ax> 0 holds for any x∈ℝ^n.  For any positive definite matrix A, we define the corresponding ellipsoid by
𝖤𝗅𝗅𝗂𝗉(A)≜{ x :  x^A^-1x≤ 1 }. 
















 §.§ Graph Laplacian


Let
  G=(V,E,w) be a connected,   undirected and weighted  graph with n vertices, m edges, and weight function w: E→ℝ_≥ 0. We fix an arbitrary orientation of the edges in G, and let B∈^m× n be the signed edge-vertex incidence matrix defined by 
 
    B_G(e,v)={
             1        
    
             -1        
    
                     0       .

We  define an m× m diagonal matrix W_G by W_G(e,e)=w_e for any edge e∈ E[G].

 
 The Laplacian matrix of G is an n× n matrix L defined by

    L_G(u,v)={
             -w(u,v)        
    (u)        
    
                     0       .

where (v)=∑_u∼ v w(u,v).
It is easy to verify 
 that 
    x^L_G x = x^B_G^W_GB_Gx =∑_u∼ v w_u,v(x_u-x_v)^2 ≥ 0
 holds
for any x∈ℝ^n. Hence, the Laplacian matrix of any undirected graph is a  matrix.  
Notice that, by setting x_u=1 if u∈ S and x_u=0 otherwise, x^L_Gx equals to the value of the cut between S and V∖ S. Hence, a spectral sparsifier is a stronger notion than a cut sparsifer.










 §.§ Other notations

For any sequence {α_i}_i=1^m, we use nnz(α) to denote the number of non-zeros in {α_i}_i=1^m.
For any two matrices A and B, we write A≼ B to represent B-A is , and A≺ B to represent B-A
 is positive definite.  For any two matrices A and B of the same dimension, let A∙ B≜(A^B), and
 
    A⊕ B =
     [ [ A 0; 0 B;   ]].

 





§ OVERVIEW OF OUR ALGORITHM


Without loss of generality we study the problem of sparsifying the sum of  matrices. The one-to-one correspondence between the
construction of a graph sparsifier and the following Problem <ref> was presented in <cit.>.


Given a set S of m  matrices 
 M_1,⋯,M_m with ∑_i=1^m M_i=I and 0<<1, find non-negative coefficients {c_i}_i=1^m such that
| { c_i | c_i≠ 0}| =O(n/^2), and

    (1-)· I≼∑_i=1^m c_i M_i ≼
    (1+)· I.



For intuition, one can think all M_i are rank-1 matrices, i.e., M_i = v_i v^_i for some v_i∈^n.
Given the correspondence between  matrices and ellipsoids, Problem <ref> essentially asks to use O(n/^2) vectors from S to construct an ellipsoid, whose shape is close to be  a sphere.
To construct such an ellipsoid with desired shape, all previous algorithms <cit.> proceed by iterations:  in each iteration j the
algorithm chooses  one or more vectors, denoted by v_j_1,⋯, v_j_k, and adds Δ_j≜∑_t=1^k v_j_t v^_j_t to the 
 currently constructed matrix by setting A_j=A_j-1+Δ_j.  To control the shape of the constructed ellipsoid,  two barrier values, the upper barrier u_j and the lower barrier ℓ_j, are maintained  such that  the constructed ellipsoid (A_j) 
is  sandwiched between  the outer sphere u_j· I and the  inner sphere ℓ_j· I for any iteration j. That is, the following invariant always maintains:

    ℓ_j· I≺ A_j≺  u_j· I.



To ensure  invariant holds, two barrier values ℓ_j and u_j are increased properly after each iteration, i.e., 

    u_j+1= u_j +δ_u,j,     ℓ_j+1 = ℓ_j + δ_ℓ,j

 for some positive values δ_u,j and δ_ℓ,j. The algorithm  continues this process, until  after T iterations (A_T) is close to be a sphere. This implies that A_T is a solution of  
 Problem <ref>, see bsspic for an illustration.
  

However, turning the scheme  described above into an efficient algorithm 
we need to consider the following  issues:


  * Which vectors  should we pick in each iteration?

  * How many vectors can be added in each iteration? 

  * How should we update u_j and ℓ_j properly so that the invariant invariant always holds? 

These three questions closely relate to each other: on  one hand, one can always  pick a single “optimal” vector in each iteration based on some metric, and such conservative approach requires a linear number of iterations T=Ω(n/^2) and
super-quadric time for each iteration. On the other hand, one can choose multiple less “optimal" vectors to construct Δ_j in iteration j, but this makes the update of barrier values more challenging to ensure the invariant invariant holds. 
Indeed, the previous constructions <cit.> speed up their algorithms at the cost of increasing the sparsity, i.e., the number of edges in a sparsifier, by more than a multiplicative constant.















To address these, we introduce three novel techniques for constructing a spectral sparsifier:
First of all, we define a new potential function which is much easier to compute yet has similar guarantee as the potential function introduced in <cit.>.

Secondly we show that solving Problem <ref> with two-sided constraints in sscondition can be reduced to a similar problem  with only one-sided constraints.  Thirdly we prove that the problem with one-sided constraints can be solved by a semi-definite program.



 §.§ A new potential function


To ensure that the constructed ellipsoid A is always inside the outer sphere u· I, we introduce a potential function Φ_u(A) defined by 

    Φ_u(A)≜exp((uI-A)^-1)=∑_i=1^n exp( 1/u-λ_i(A)).

It is easy to see that, when  (A) gets closer to the outer sphere, λ_i(u· I-A) becomes smaller and the value of  Φ_u(A) increases. Hence, a bounded value of Φ_u(A) ensures that (A) is inside the sphere u· I. For the same reason, we introduce a potential function Φ_ℓ(A)
defined by 

    Φ_ℓ(A)≜exp((A-ℓ I)^-1)=∑_i=1^n exp( 1/λ_i(A)-ℓ).

to ensure that the inner sphere is always inside (A). We also define 

    Φ_u,ℓ(A)≜Φ_u(A) + Φ_ℓ(A),

as a bounded value of Φ_u,ℓ(A)
 implies that the two events occur simultaneously. Our goal is to design a proper update rule 
 to construct {A_j} inductively, so that Φ_u_j, ℓ_j(A_j) is monotone non-increasing after each iteration. Assuming this, a bounded value of the initial potential function guarantees that the invariant invariant always holds.
 
 To analyse the change of the potential function, we first notice that
 
    Φ_u,ℓ(+Δ)  ≥Φ_u,ℓ() + ( e^(u-)^-1 (u-)^-2Δ) - (e^(-ℓ)^-1(-ℓ)^-2Δ)

by the convexity of the function Φ_u,ℓ.
We  prove that, as long as the matrix Δ satisfies 0≼Δ≼δ(uI-)^2 and 0≼Δ≼δ(-ℓ I)^2 for some small δ, the first-order approximation gives  a good approximation.
 
 

Let A be a symmetric matrix. Let u,ℓ be the barrier values such that u-ℓ≤1 and ℓ≺≺ u.
Assume
that Δ≻ 0,
 Δ≼δ(uI-)^2 and Δ≼δ(-ℓ I)^2
for δ≤ 1/10. Then, it holds that

    Φ_u,ℓ(+Δ)  ≤   Φ_u,ℓ() + (1+2δ)( e^(u-)^-1 (u-)^-2Δ)
            -(1-2δ)(e^(-ℓ)^-1(-ℓ)^-2Δ).



We remark that this is not the first paper to use a potential function to guide the growth of the ellipsoid. In <cit.>, the potential function

    Λ_u,ℓ,p() = ((uI-A)^-p) + ((A-ℓ I)^-p)

is used
with p= 1. The main drawback is that Λ_u,ℓ,1 does not differentiate the
following two cases: 


  * Multiple eigenvalues of A are close to the boundary (both u
and ℓ).

  * One of the eigenvalues of A is very close to the boundary (either
u or ℓ).

It is known that, when one of the eigenvalues of A is very close to the boundary, it is more difficult to find 
an “optimal” vector. 
<cit.> shows this problem can be alleviated by using p ≫ 1. 
However, this choice of p makes the function Λ_u,ℓ,p less smooth and hence one has to take a smaller step size δ, as shown in the following lemma by <cit.>.



Let A be a symmetric matrix and Δ be a rank-1 matrix. Let u,ℓ be the barrier values such that ℓ≺≺ u.
Assume that Δ≻ 0, Δ≼δ(uI-) and Δ≼δ(-ℓ I)
for δ≤ 1/(10 p) and p ≥ 10. Then, it holds that 

    Λ_u,ℓ,p(+Δ) ≤   Λ_u,ℓ,p() + p (1+p δ)((u-)^-(p+1)Δ)- p(1-p δ)((-ℓ)^-(p+1)Δ).



Notice that, comparing with the potential function oldpotential, our new potential function defpotential blows up much faster when the eigenvalues of A are closer to the boundaries ℓ and u.
This allows the problem of finding an “optimal” vector much easier than using Λ_u,ℓ,p. At the same time, we avoid the problem of taking a small step Δ≼ 1/p· (uI-) by taking a “non-linear" step Δ≼ (uI-)^2. As there cannot be too many eigenvalues close to the boundaries, this “non-linear” step allows us to take a large step except on a few directions.




 §.§ A simple construction based on 


The second technique we introduce is the reduction from 
a  spectral sparsifier with two-sided constraints
 to the one with one-sided constraints.  Geometrically,
it is equivalent to require the constructed ellipsoid inside another ellipsoid, instead of being sandwiched between two spheres as depicted in bsspic.
Ideally, we want to reduce the two-sided problem to the following problem: for a set of  matrices ={M_i }_i=1^m such that ∑_i=1^mM_i≼ I,  find a sparse representation Δ=∑_i=1^m α_i M_i with small (α) such that Δ≼ I. 
However, in the reduction we use such matrix  Δ=∑_i=1^m α_i M_i to update A_j and we need the length of (Δ) is large on the direction that (A_j) is small.

 To encode this information, we define the generalised one-sided problem as follows:












[One-sided Oracle]
Let  ≼ B ≼ I, ≽, ≽ be  symmetric matrices,  
ℳ={M_i}_i=1^m be  a set of matrices such that ∑_i=1^m_i =. 
We call a randomised algorithm (, B,, ) a one-sided oracle with speed S∈ (0,1] and error >0, if (, B,, )  outputs  a matrix Δ=∑_i=1^mα_i_i such that 


  * nnz(α)≤λ_min () ·(^-1).

  * Δ≼ and α_i ≥ 0 for all i.

  * [∙Δ] ≥ S ·λ_min() ·() -  S ·λ_min() ·(), where =-, and =+.


We  show in impleSDP the existence of a one-sided oracle with  speed S=Ω(1) and error =0, in which case the oracle only requires  as input, instead of   and . However, to construct such an oracle efficiently an additional error is introduced, which depends on +.



For the main algorithm 𝚂𝚙𝚊𝚛𝚜𝚒𝚏𝚢( ,ε), we  maintain
the matrix A_j inductively as we 
 discussed at the beginning of overview. By employing the ⊕ operator, we reduce the problem of constructing Δ_j with two-sided constraints to the problem of  constructing Δ_j⊕Δ_j with  one-sided constraints. 
We also use   to  ensure that the length of (Δ) is large on the direction where the length of (A_j) is small. See Algorithm <ref> for formal description.







To analyse Algorithm <ref>, we use the fact that the returned Δ_j satisfies the preconditions of potentialchange and prove in potential_decreasing that


    [Φ_u_j+1,ℓ_j+1(_j+1)]≤Φ_u_j,ℓ_j(_j)

for any iteration j.
Hence, with high probability the bounded ratio between u_T and ℓ_T after the final iteration T implies that the (A_T) is close to be a sphere. In particular,  for any < 1/20, a (1+O())-spectral sparsfier can be constructed by calling    O(log^2n/ε^2· S)
times, which is  described in the lemma below.





Let 0<ε<1/20. Suppose we have one-sided oracle 𝙾𝚛𝚊𝚌𝚕𝚎 with speed S and error . Then,
with constant probability 
 the algorithm
𝚂𝚙𝚊𝚛𝚜𝚒𝚏𝚢( ,ε) outputs a (1+O(ε))-spectral sparsifier with  O(n/ε^2· S) vectors
by calling  O(log^2n/ε^2· S)
times.






 
 
 





 §.§ Solving  via SDP









Now we  show that the required solution of  (, B, ) indeed exists[As the goal here is to prove the existence of  with error =0, the input here is  instead of C_+ and C_-.], and can be further solved in nearly-linear time by  a semi-definite program.
We first prove that the required matrix satisfying the conditions of Definition <ref> exists for some absolute constant S=Ω(1) and =0. To make a parallel discussion between Algorithm <ref> and the algorithm we will present later,
 we  use A to denote the output of the  instead of Δ. We adopt the ideas between the ellipsoid and two spheres discussed before, but only  consider one sphere for the one-sided case. Hence, we introduce a barrier value u_j for each iteration j, where u_0=0. We will use the  potential function 

    Ψ_j=(u_jB-A_j )^-1

in our analysis, where u_j  is increased by
δ_j≜ (Ψ_j·λ_min(B))^-1
after  iteration j. 
Moreover, since we only need to prove the existence of the required matrix A=∑_i=1^m α_i M_i, our process proceeds for T iterations, where only one vector is chosen  in each iteration.  To find  a desired vector, we perform random sampling, where  each matrix M_i is sampled with probability  𝗉𝗋𝗈𝖻(M_i) proportional to M_i∙, i.e.,

    𝗉𝗋𝗈𝖻(M_i)≜(M_i∙)^+/∑_t=1^m (M_t∙)^+,

where x^+≜max{x,0}. Notice that, since 
our goal is to construct A such that 𝔼[∙ A] is lower bounded by some threshold as stated in Definition <ref>, we should not pick any matrix M_i with M_i∙ < 0. This random sampling procedure is described in Algorithm <ref>, and the properties of the output matrix is summarised in one_sided. 

 



 
Let 0≼ B≼ I and  be symmetric matrices, and ={M_i}_i=1^m be a set of  matrices such that 

∑_i=1^m M_i = I. Then 
𝚂𝚘𝚕𝚞𝚝𝚒𝚘𝚗𝙴𝚡𝚒𝚜𝚝𝚎𝚗𝚌𝚎(,,) outputs
a matrix =∑_i=1^mα_i _i such that the following  holds:


  * nnz(α)=⌊λ_min()·(^-1)⌋.

  * ≼, and α_i ≥ 0 for all i.

  * [∙]≥1/32·λ_min()·( ).



one_sided shows that the required matrix A defined in Definition <ref> exists, and can be found by  random sampling  described in Algorithm <ref>. Our key observation is that such matrix A can be constructed  by  a semi-definite program.



Let 0≼ B≼ I,  be symmetric matrices, and ={M_i}_i=1^m be a set of matrices such that ∑_i=1^m M_i = I.

Let S⊆ [m] be a random set of ⌊λ_min()·(^-1)⌋ coordinates, where every index i is picked with probability 𝗉𝗋𝗈𝖻(M_i). Let ^⋆
be the solution of the following semidefinite program

    max_α_i≥0∙(∑_i∈ Sα_i_i) subject to =∑_i∈ Sα_i_i≼

Then, we have [∙^⋆] ≥1/32·λ_min()· ().


Taking the SDP formuation SDP and the specific constraints of the 's input into account, the next lemma  shows that the required matrix used in each iteration of   can be computed efficiently by solving a semidefinite program.


The Oracle used in Algorithm 
can be implemented in 
    Õ((Z+n^ω)·ε^-O(1))  work  and  Õ(ε^-O(1)) depth

where Z=∑_i=1^mnnz(M_i) is the total number of non-zeros in _i. When the
matrix ∑_i=1^m M_i= comes from spectral sparsification of graphs,
each iteration of 
can be implemented in 

    Õ(mε^-O(1))work  and  Õ(ε^-O(1)) depth.

Furthermore, the speed of this one-sided oracle is Ω(1) and the error of this one-sided oracle is . 





Combining reductionlem and oracle_implementation gives us the proof of the main result.



oracle_implementation shows that we can construct  with Ω(1) speed and  error that runs in 

    Õ((Z+n^ω)·ε^-O(1))  work  and  Õ(ε^-O(1)) depth

for the matrix setting and 

    Õ(mε^-O(1))work  and  Õ(ε^-O(1)) depth.

for the graph setting. 
Combining this with   reductionlem, which states that  it suffices to call   Õ(1/^2) times,
the main statements hold.







 §.§ Further discussion
  

Before presenting a more detailed analysis of our algorithm,  we   compare our new approach with the previous ones for constructing a linear-sized spectral sparsifier, and see how we 
address the bottlenecks faced in previous constructions. Notice that all 
 previous algorithms require super poly-logarithmic number of  iterations, and super linear-time for each iteration. For instance, our previous algorithm <cit.> for constructing a sparsifier with O(pn) edges requires Ω(n^1/p) iterations and  Ω(n^1+1/p) time per iteration for the following reasons:


  * Ω(n^1+1/p) time is needed per iteration: Each iteration takes n/g^Ω(1) time to pick the vector(s) when (ℓ + g) I ≼ A ≼ (u-g) I. To avoid eigenvalues of A getting too close to the boundary  u or ℓ, i.e., g being too small, we choose the potential function whose value dramatically increases when the eigenvalues of A get close u or ℓ. As the cost, we need to scale down the added vectors by a n^1/p factor. 
 
  * Ω(n^1/p) iterations are needed:  By random sampling, we choose O(n^1-1/p) vectors each iteration and use the matrix Chernoff bound to show that the “quality” of added  O(n^1-1/p) vectors is just p=Θ(1) times worse than adding a single vector. Hence, this requires Ω(n^1/p) iterations.

In contrast, our new approach breaks these two barriers through the following way:


  * A “non-linear” step: Instead of rescaling down the vectors we add uniformly, we pick much fewer vectors on the direction that blows up, i.e., we impose the condition Δ≼ (uI-A)^2 instead of Δ≼ 1/p· (uI-A). This allows us to use the new potential function defpotential with form exp(x^-1) to control the eigenvalues in a more aggressive way.

  * SDP filtering:  By matrix Chernoof bound, we know that the probability that we sample a few “bad” vectors is small. Informally,
 we apply semi-definite programming to filter out those bad vectors, and this allows us to add Ω(n/log^O(1)(n)) vectors in each iteration.




§ DETAILED ANALYSIS 


In this section we give  detailed analysis for the statements presented in overview.




 §.§ Analysis of the potential function


Now we analyse the properties of the potential function defpotential, and prove potentialchange.
 The following two facts from matrix analysis will be used in our analysis.


 
Let A∈ℝ^n× n, U∈ℝ^n× k, C∈ℝ^k× k and
V∈ℝ^k× n be matrices. Suppose that A, C and C^-1 + VA^-1U are invertible, it holds that

    (A+UCV)^-1= A^-1 - A^-1U( C^-1 + VA^-1U )^-1 VA^-1.




It holds for any symmetric matrices A and B that

    (e^A+B) ≤(e^A·e^B ).





We analyse the change of Φ_u(·) and Φ_ℓ(·) individually.  First of all, notice that

    (uI-A-Δ)^-1=(uI-A)^-1/2( I-(uI-A)^-1/2Δ(uI-A)^-1/2)^-1(uI-A)^-1/2.

We define Π=(uI-A)^-1/2Δ(uI-A)^-1/2. Since 0≼Δ≼δ(uI-A)^2 and u-ℓ≤ 1, it holds that

    Π≼δ (uI-A)≼δ (uI-ℓ I) ≼δ I,

and therefore

    (I-Π)^-1≼ I+ 1/1-δ·Π.

Hence, it holds that

    (uI-A-Δ)^-1   ≼ (uI-A)^-1/2( I+ 1/1-δ·Π)(uI-A)^-1/2
        = (uI-A)^-1 + 1/1-δ· (uI-A)^-1Δ(uI-A)^-1.

By the fact that exp is monotone and the Golden-Thompson Inequality (Golden), we have that

    Φ_u(A+Δ)    = exp( (uI-A-Δ)^-1) 
       ≤exp( (uI-A)^-1 + 1/1-δ· (uI-A)^-1Δ(uI-A)^-1) 
       ≤(exp(uI- A)^-1exp(1/1-δ· (uI- A)^-1Δ(uI-A)^-1)).

Since 0 ≼Δ≼δ (uI-A)^2 and δ≤ 1/10 by assumption, we have that
(uI-A)^-1Δ(uI-A)^-1≼δ I,
and

    exp(1/1-δ· (uI-A)^-1Δ(uI-A)^-1)≼ I+(1+2δ)· (uI-A)^-1Δ(uI-A)^-1.

Hence, it holds that

    Φ_u(+Δ)    ≤(^(uI- A)^-1·( I+(1+2δ)(uI-A)^-1Δ(uI-A)^-1) )
        =  Φ_u()+(1+2δ)·(e^(uI-A )^-1(uI-A)^-2Δ).


By the same analysis, we have that 

    Φ_ℓ(+Δ)    ≤(^(A-ℓ I)^-1·( I-(1-2δ)(A- ℓ I)^-1Δ(A-ℓ I)^-1) )
        =  Φ_ℓ()-(1-2δ)·(e^(A-ℓ I )^-1(A- ℓ I)^-2Δ).

Combining the analysis on Φ_u(A+Δ) and Φ_ℓ(A+Δ) finishes the proof. 




Let A be a symmetric matrix. Let u,ℓ be the barrier values such that u-ℓ≤1 and ℓ≺≺ u. Assume
that 0≤δ_u≤δ·λ_min(uI-A)^2 and 0≤δ_ℓ≤δ·λ_min(A-ℓ I)^2 for δ≤ 1/10. Then, it holds that 

    Φ_u+δ_u,ℓ+δ_ℓ()  ≤   Φ_u,ℓ()  - (1-2δ) δ_u·(^(u-)^-1(u-)^-2)
            +(1+2δ)δ_ℓ·(^(-ℓ)^-1(-ℓ)^-2).




Since 0≤δ_u≤δ·λ_min(uI-A)^2 and 0≤δ_ℓ≤δ·λ_min(A-ℓ I)^2,
 we have that δ_u· I≼δ· (uI-A)^2 and
δ_ℓ· I≼δ· (A-ℓ I)^2.
The statement follows by a similar analysis for proving potentialchange.




 §.§ Analysis of the reduction


Now we present the detailed analysis for the reduction from a spectral sparsifier to a one-sided oracle. We first analyse Algorithm <ref>, and prove that in expectation the value of the potential function is not increasing. Based on this fact, we will give a proof of reductionlem, which shows that a (1+O(ε))-spectral sparsifier can be constructed by calling  O(log^2n/ε^2· S)
times.



 Let A_j and A_j+1 be the matrices constructed by Algorithm <ref> in iteration j and j+1, and 
assume that 0 ≤ε≤ 1/20. Then, it holds that

    [Φ_u_j+1,ℓ_j+1(_j+1)]≤Φ_u_j,ℓ_j(_j).



By the description of Algorithm <ref> and Definition <ref>, it holds that

    Δ_j⊕Δ_j ≼ (u_j-_j)^2⊕(_j-ℓ_j)^2,

which implies that
Δ_j≼  (u_j-_j)^2 and Δ_j≼  (_j-ℓ_j)^2. Since u_j-ℓ_j ≤ 1 by the algorithm description and 0≤≤ 1/20, by   setting Δ=ε·Δ_j in potentialchange, we have 

    Φ_u_j, ℓ_j (_j + ε·Δ_j)
       ≤Φ_u_j,ℓ_j(_j)   + ε (1+2ε)·(e^(u_j-_j)^-1(u_j-_j)^-2Δ_j)
                - ε (1-2ε)·(e^(_j-ℓ_j)^-1(_j-ℓ_j)^-2Δ_j)
        = Φ_u_j,ℓ_j(_j)-ε·∙Δ_j.

Notice that the matrices {M_i⊕ M_i}_i=1^m as the input of 𝙾𝚛𝚊𝚌𝚕𝚎 always satisfy ∑_i=1^mM_i⊕ M_i=I⊕ I.
Using this and the definition of 𝙾𝚛𝚊𝚌𝚕𝚎, we know that

    [∙Δ_j] ≥ 
     S·λ_min(B_j)· ()-· S·λ_min(B_j)·( ),

Let α_j=ε· S·λ_min(_j). Then we have that 

    [Φ_u_j,ℓ_j(_j+1) ]    ≤Φ_u_j,ℓ_j(_j) - ε·[∙Δ_j]
       ≤Φ_u_j,ℓ_j(_j)  +  (1+2ε) (1+)·α_j·(e^(u_j-_j)^-1(u_j-_j)^-2) 
                - (1-2ε)  (1-)·α_j·(e^(_j-ℓ_j)^-1(_j-ℓ_j)^-2).


On the other hand, using that 0≤≤ 1/20, S ≤ 1 and Δ_j≼  (u_j-_j), we have that

    δ_u,j≤ε·(1+2ε)(1+)/1-4ε·λ_min(u_j-_j)^2≤ 2 ·λ_min(u_j-_j+1)^2

and

    δ_ℓ,j≤ε·(1-2ε)(1-)/1+4ε·λ_min(_j-ℓ_j)^2≤ 2 ·λ_min(_j+1-ℓ_j)^2.















Hence, potential2 shows that

    Φ_u_j +δ_u,j,ℓ_j+δ_ℓ,j(_j+1)   ≤Φ_u_j,ℓ_j(_j+1)
     - (1 - 4ε)δ_u,j·(e^(u_j-_j+1)^-1(u_j-_j+1)^-2)
                     + (1 + 4ε)δ_ℓ,j·(e^(_j+1-ℓ_j)^-1(_j+1-ℓ_j)^-2) 
       ≤Φ_u_j,ℓ_j(_j+1) - (1 - 4ε)δ_u,j·(e^(u_j-_j)^-1(u_j-_j)^-2)
                    + (1 + 4ε)δ_ℓ,j·(e^(_j-ℓ_j)^-1(_j-ℓ_j)^-2).

By combining echange1, echange2, and setting (1-4ε)δ_u,j=(1+2ε)(1+)α_j,
 (1+4ε)δ_ℓ,j=(1-2ε)(1-)α_j, we have that
𝔼[Φ_u_j+1,ℓ_j+1(_j+1)]≤Φ_u_j,ℓ_j(_j).




We first bound the number of times the algorithm calls the oracle. 
Notice that
Φ_u_0, ℓ_0 = 2·exp((1/4I)^-1) = 2 ^4· n.
Hence, by potential_decreasing we have  𝔼[ Φ_u_j, ℓ_j (A_j) ]=O(n) for any iteration j. By Markov's inequality, it holds that Φ_u_j, ℓ_j   (A_j)=n^O(1) with high probability in n.
In the remainder of the proof, we assume that this event occurs.

Since _j=(u_j-_j)^2⊕(_j-ℓ_j)^2 by definition, it holds that 

    exp(  (λ_min(_j ))^-1/2) ≤Φ_u_j, ℓ_j(A_j) =   n^O(1),

which implies that 

    λ_min( _j )=Ω(log^-2n).

On the other hand, in iteration j the gap between u_j and ℓ_j is increased by

    δ_u,j-δ_ℓ,j=Ω(ε^2· S·λ_min(_j) ).

Combining this with lblambdab gives us that
 
    δ_u,j-δ_ℓ,j  =  Ω(ε^2· S/log^2n)

for any j.
Since u_0-ℓ_0=1/2 and the algorithm terminates once u_j - ℓ_j>1 for some j, with high probability in n, the algorithm terminates in
O( log^2 n/ε^2· S)
iterations.

Next we prove that the number of M_i's involved in the output is at most  O(n/^2· S).
By the properties of 𝙾𝚛𝚊𝚌𝚕𝚎, the number of matrices in iteration j is at most 
λ_min(B_j)·( B_j^-1). Since x^-2≤exp( x^-1) for all x>0, it holds for any iteration j that

    (_j^-1)     =((u_j-_j)^-2) +((_j-ℓ_j)^-2) ≤Φ_u_j,ℓ_j(_j)

By gapincrease, we know that for added matrix M_i in iteration j, the average increase of the  gap u_j-ℓ_j for each added matrix is 
Ω(ε^2· S/Φ_u_j,ℓ_j(_j)).
Since 𝔼[Φ_u_j, ℓ_j(A_j)]=O(n), for every new added matrix, in expectation the gap between u_j and ℓ_j is increased by
Ω(ε^2· S/ n).
By the ending condition of the algorithm, i.e., u_j-ℓ_j>1, and Markov's inequality, the number of matrices picked in total is at most 
O (n/ε^2· S)
with constant probability.

Finally we prove that the output is a (1+O(ε))-spectral sparsifier. Since the condition number of the output matrix A_j is at most

    u_j/ℓ_j =
    ( 1 - u_j-ℓ_j/u_j)^-1,

it suffices to prove that
(u_j-ℓ_j)/u_j=O(ε) and this  easily follows from the ending condition of the algorithm and

    δ_u,j -δ_ℓ, j/δ_u,j= O().































 §.§ Existence proof for 




The property on nnz(α) follows from the algorithm description.
For the second property, notice that  every chosen matrix Δ_j in iteration j satisfies 
Δ_j≼1/2(u_j-_j),
which implies that 
_j≼ u_j holds for any iteration j. Hence,
=1/u_T_T≼, and α_i ≥ 0 since Ψ_j≥ 0. 


Now we prove the third statement. Let 

    β=∑_i=1^m (_i∙)^+.

Then, for  each matrix _i_j picked in iteration j, ∙_j is increased by

    ∙Δ_j    =1/4Ψ_j·𝗉𝗋𝗈𝖻(M_i_j)·∙_i_j=β/4Ψ_j.

On the other hand, it holds that

    u_T= u_0 + ∑_j=0^T-1δ_j = 1+∑_j=0^T-1( Ψ_j·λ_min(B) )^-1.

Hence,
we have that

    ∙    =1/u_T· C∙(∑_j=0^T-1Δ_j ) 
    = ∑_j=0^T-1β·(4Ψ_j)^-1/1+∑_j=0^T-1 (Ψ_j·λ_min())^-1
        =  βλ_min() /4·∑_j=0^T-1Ψ_j^-1/λ_min() +∑_j=0^T-1Ψ_j^-1
       ≥βλ_min()/4·∑_j=0^T-1(Ψ_j+Ψ_0)^-1/λ_min() +∑_j=0^T-1(Ψ_j+Ψ_0)^-1
       ≥βλ_min() /4·∑_j=0^T-1(Ψ_j+Ψ_0)^-1/λ_min() + T ·Ψ_0^-1≥β/8·∑_j=0^T-1(Ψ_j+Ψ_0)^-1,

where the last inequality follows by the choice of T.
Hence, it suffices to bound Ψ_j. 

Since Δ_j≼1/2(u_j-_j)≼1/2(u_j+1-_j),
we have that 

    Ψ_j+1≤((u_j+1-_j)^-1)+2·Δ_j∙(u_j+1-_j)^-2.

Since (u-_j)^-1 is convex in u, we have that

    ((u_j-_j)^-1) ≥((u_j+1-_j)^-1) + δ_j ((u_j+1-_j)^-2 B)











Combining mideq1 and mideq2, we have that 

    Ψ_j+1   ≤((u_jB-A_j)^-1) -δ_j·tr( ( u_j+1 B -A_j)^-2B ) +2·Δ_j∙(u_j+1-_j)^-2
        = Ψ_j  -δ_j·λ_min(B)·((u_j+1B-A_j)^-2)+2·Δ_j∙(u_j+1-_j)^-2

Let ℰ_j be the event that Δ_j≼1/2(u_j-_j).  Notice that our picked Δ_j in each iteration always satisfies ℰ_j by algorithm description.
 Since 

    [Δ_j∙(u_j-_j)^-1]=∑_i: (M_i∙ C)^+>0𝗉𝗋𝗈𝖻(M_i)·1/4Ψ_j·𝗉𝗋𝗈𝖻(M_i)  ·_i∙(u_j-_j)^-1≤1/4,

by Markov inequality it holds that 

    ℙ[ℰ_j]= ℙ[ Δ_j≼1/2(u_j-_j)] ≥1/2,

and therefore

    [Δ_j∙(u_j+1-_j)^-2 | ℰ_j]    ≤[Δ_j∙(u_j+1-_j)^-2]/ℙ(ℰ_j)
       ≤ 2·[Δ_j∙(u_j+1-_j)^-2]
       = 1/2Ψ_j·∑_i: (M_i∙ C)^+>0_i∙(u_j+1-_j)^-2
       ≤1/2Ψ_j·(u_j+1-_j)^-2.
 
Combining the inequality above,  upphinext, and the fact that every Δ_j picked by the algorithm satisfies ℰ, we have that 

    [ Ψ_j+1]≤Ψ_j+(1/Ψ_j-δ_j·λ_min())·(u_j+1-_j)^-2.

By our choice of δ_j, it holds for any iteration j that [ Ψ_j+1] ≤Ψ_j, and 

    [ (Ψ_j+1 + Ψ_0 )^-1]≥(Ψ_j + Ψ_0 )^-1≥1/2·Ψ_0.

Combining this with CA_bound, it holds that 

    [∙] ≥β/8∑_j=0^T-1[(Ψ_j+Ψ_0)^-1]
     ≥β/16·T/Ψ_0 = β/16·T/( ^-1)≥() /16·T/( ^-1) .

The result follows from the fact that T ≥λ_min() ( ^-1)/2.

Using the lemma above, we can prove that such  can be solved
by a semidefinite program.



Note that the probability we used in the statement is the same as the probability we used in 𝚂𝚘𝚕𝚞𝚝𝚒𝚘𝚗𝙴𝚡𝚒𝚜𝚝𝚎𝚗𝚌𝚎(
ℳ,,). Therefore, Lemma <ref> shows that there is a matrix A of the form ∑_i=1^mα_i_i such that 

    [ ∙] ≥1/32·λ_min()·(),

The statement follows by the fact that A^⋆ is the solution of the semidefinite program SDP that maximises ∙.





 §.§ Implementing the SDP in nearly-linear time


Now, we discuss how to solve the SDP SDP in nearly-linear time.
Since this SDP is a packing SDP, it is known how to solve it in nearly-constant depth <cit.>.
The following result will be used in our analysis.
Given a SDP 

    max_x≥0c^x subject to ∑_i=1^mx_i_i≼

with _i≽, ≽ and c∈^m.
Suppose that we are given a direct access to the vector c∈^m and an
indirect access to _i and  via an oracle 𝒪_L,δ
which inputs a vector x∈^m and outputs a vector v∈^m
such that

    v_i∈(1±δ/2)[_i∙^-1/2exp(L·^-1/2(∑_ix_i_i-)^-1/2)^-1/2]

in 𝒲_L,δ work and 𝒟_L,δ depth for any x such that x_i≥0
and ∑_i=1^mx_i_i≼2. Then, we can output
x such that 

    [c^x ]≥(1-O(δ))𝖮𝖯𝖳  with  ∑_i=1^mx_i_i≼

 in

    O(𝒲_L,δlog m·log(nm/δ)/δ^3) work  and   O(𝒟_L,δlog m·log(nm/δ)/δ) depth,

where L=(4/δ)·log(nm/δ).













Since we are only interested in a fast implementation of the one-sided oracle used 
in ,
it suffices to solve the SDP SDP for this particular situation. 










































































Our basic idea is to use solve_SDP as the one-sided oracle.
Notice that each iteration of 
uses the one-sided oracle with the input 

    =1-2ε/2((-ℓ)^-2exp(-ℓ)^-1⊕(-ℓ)^-2exp(-ℓ)^-1),
        =1+2ε/2((u-)^-2exp(u-)^-1⊕(u-)^-2exp(u-)^-1),
    
    B     =(u-)^2⊕(-ℓ)^2,

where we drop the subscript j indicating iterations here for simplicity.
To apply sdp, we first sample a subset S⊆ [m],
then solve the SDP

    max_β_i≥0,β_i=0 on i∉ S
    (-)∙(∑_i=1^mβ_i_i⊕_i)   subject to ∑_i=1^mβ_i_i⊕_i≼.

By ignoring the matrices _i with i∉S, this SDP
is equivalent to the SDP 

    max_β_i≥0c^⊤β   subject to ∑_i=1^mβ_i_i⊕_i≼

where c_i=(-)∙(_i⊕_i). Now assume that  
 (i) we can approximate c_i with δ(+)∙(_i⊕_i) additive
error, and (ii) for any x such that ∑_i=1^mx_i_i⊕_i≼2 and L =Õ(1/δ),
we can approximate

    (_i⊕_i)∙^-1/2exp(L·^-1/2(∑_ix_i_i⊕_i-)^-1/2)^-1/2

with 1±δ multiplicative error.
Then, by solve_SDP we can find a vector β
such that

    [ c^β]    ≥(1-O(δ))𝖮𝖯𝖳-δ∑_i=1^mβ_i(+)∙(_i⊕_i)
       ≥(1-O(δ))𝖮𝖯𝖳-δ∑_i=1^mβ_i(+)∙
       ≥𝖮𝖯𝖳-O(δ)(+)∙.

where we used that _i⊕_i≼ and 𝖮𝖯𝖳≤(+)∙.
Since u-ℓ≤1, we have that B≼ I⊕ I and hence

    [ c^β]    ≥𝖮𝖯𝖳-O(δ)·(+)
       ≥1/32λ_min()·(()-O(δlog^2n)·(+))

where we apply sdp and (<ref>) at
the last line. Therefore, this gives an oracle with speed 1/32
 and ε error by setting δ=ε/log^2n.

The problem of approximating sample probabilities, {c_i}, as well as implementing the oracle 𝒪_L,δ
is similar with approximating leverage scores <cit.>, and relative leverage scores <cit.>. All these references  use 
the Johnson-Lindenstrauss lemma

to reduce the problem of approximating
matrix dot product or trace to matrix vector multiplication. The only
difference is that, instead of computing (-ℓ)^-(q+1)x and (u-)^-(q+1)x for a given vector x in other references, we compute (-ℓ)^-2exp(-ℓ)^-1x and 
(u-)^-2exp(u-)^-1x. These can be approximated by Taylor expansion and the number of terms required for Taylor expansion depends on how close the eigenvalues
of A are to the boundary (u or ℓ). In particular, we show in taylor that Õ(1/g^2) terms in Taylor expansion suffices,
where the gap g is the largest number such that (ℓ +g)I ≼ A ≼ (u - g)I. Since 1/g^2 = Õ(1) by (<ref>), each iteration can be implemented via solving Õ(1/^O(1))  linear systems and Õ(1/^O(1))  matrix vector multiplication. For the matrices coming from graph sparsification, this can be done in nearly-linear work and nearly-constant depth <cit.>. For general matrices, this can be done in input sparsity time and nearly-constant depth <cit.>.





 §.§ Taylor Expansion of x^-2exp(x^-1) 


Suppose f is holomorphic
on a neighborhood of the ball B≜{z∈ℂ : |z-s|≤ r},
then we have that

    |f^(k)(s)|≤k!/r^ksup_z∈ B|f(z)|.



Let f(x)=x^-2exp(x^-1). For any 0<x≤ 1,
we have that

    |f(x)-∑_k=0^d1/k!f^(k)(1)(x-1)^k|≤8(d+1)e^5/x-xd.

In particular, if d≥c/x^2log(1/xε)
for some large enough universal constant c, we have that

    |f(x)-∑_k=0^d1/k!f^(k)(1)(x-1)^k|≤ε.



By the formula of the remainder term in Taylor series, we have that

    f(x)=∑_k=0^d1/k!f^(k)(1)(x-1)^k+1/d!∫_1^xf^(d+1)(s)(x-s)^dds.

For any s∈[x,1], we define D(s)={z∈ℂ : |z-s|≤ s-x/2}.
Since |f(z)|≤(x/2)^-2exp(2/x) on z∈ D(s),
Cauchy's estimates (Theorem <ref>) shows that

    |f^(d+1)(s)|≤(d+1)!/(s-x/2)^d+1sup_z∈ B(s)|f(z)|≤(d+1)!/(s-x/2)^d+14/x^2exp(2/x).

Hence, we have that

    |f(x)-∑_k=0^d1/k!f^(k)(t)(x-t)^k|≤   1/d!|∫_1^x(d+1)!/(s-x/2)^d+14/x^2exp(2/x)(x-s)^dds|
    
    =    4(d+1)e^2/x/x^2∫_x^1(s-x)^d/(s-x/2)^d+1ds
    ≤   8(d+1)e^2/x/x^3∫_x^1(1-x)^dds
    ≤    8(d+1)·e^5/x-xd.






















































































§ ACKNOWLEDGEMENT

The authors  would like to thank Michael Cohen for helpful discussions and suggesting the ideas to improve the sparsity from O(n/^3) to O(n/^2).
















§ EXPECTATION BOUND IMPLIES HIGH PROBABILITY BOUND (I AM STILL WRITING)


In this section, we assume that

Given ∑_i=1^mv_iv_i^⊤= and _1,_2,_3,⋯,_k.
Let S=∑_imax_jv_i^⊤_j^-1v_i and β=min_iλ_min_i.
The algorithm 𝚂𝚙𝚊𝚛𝚜𝚒𝚏𝚢(v_i) outputs α_i such
that (α)≤ c_1β S. Also, the matrix Δ=∑_i=1^mα_iv_iv_i^⊤
satisfies


  * c_2^-1β≼Δ≼ c_2β

  * Δ≼1/2_i for all i=1,2,⋯,k.


Note that if we sample v_iv_i^⊤ with probability max_jv_i^⊤_j^-1v_i/∑_imax_jv_i^⊤_j^-1v_i
and reweigh them by 1/max_jv_i^⊤_j^-1v_i.
Then, we have that 𝔼sampled=/S. Hence,
with β S sample, we have ∑sampled∼β.
Next, we note that

    λ_max(^-1/2sampled^-1/2)≤1.

The chernoff bound gives

    ∑sampled≼βλ_max(^-1)

if βλ_max(^-1)≥1. In particular,
we can pick

    β=1/λ_max(^-1)=λ_min().

This explains the statement up to Õ(1).

Recall the lemma

Suppose that Δ≼δ(u-) and Δ≼δ(-ℓ),
then we have that

    Φ_u,ℓ(+Δ)≤Φ_u,ℓ()+1/1-δ⟨ (u-)^-2,Δ⟩ -1/1+δ⟨ (-ℓ)^-2,Δ⟩



Although it seems that Δ≼δ(-ℓ) is not
necessary, it is. In particular, if Δ is very large with small
probability and Δ is 0 otherwise. Then, Δ is
not useful and we are not making too much progress. Hence, both condition
is important.

Now, we use the algorithm 𝚂𝚙𝚊𝚛𝚜𝚒𝚏𝚢(v_i). Hence, we have
that

    c_2^-1β≼Δ≼ c_2β

and Δ≼1/2(u-) and Δ≼1/2(-ℓ).
Hence, we have that

    Φ_u,ℓ(+Δ)    ≤   Φ_u,ℓ()+2⟨ (u-)^-2,Δ⟩ -2/3⟨ (-ℓ)^-2,Δ⟩
       ≤   Φ_u,ℓ()+2c_2β(u-)^-2-2/3c_2^-1β(-ℓ)^-2.

Assuming Φ_u,ℓ()≤1/3, we have that

    Φ_u+u',ℓ+ℓ'(+Δ)    ≤   Φ_u,ℓ(+Δ)-2u'/3(u--Δ)^-2+2ℓ'(+Δ-ℓ)^-2
       ≤   Φ_u,ℓ(+Δ)-2u'/3(u-)^-2+4ℓ'(-ℓ)^-2.

Hence, by picking 2u'/3=2c_2β and 4ℓ'=2/3c_2^-1β,
we have that the potential function is decreasing.

Hence, we have that

    u'=3c_2β and ℓ'=2/12c_2^-1β.

Hence, if δ is small and c_2 close to 1, the increase
of u and ℓ is very close.

Next, we note that whenever we add β S many vectors, 
is increased by Δ∼β. In general, we need S many
vectors to increase  by identity. Now, we need to increase 
by n. Hence, we need n many vectors.

Given the algorithm 𝚂𝚙𝚊𝚛𝚜𝚒𝚏𝚢(v_i), we can give you a
two sided algorithm. Also, if c_2∼1, then we can get a better
sparsifier.

Note that there is no way we can make it better if c_2 is not
close to 1.




§ HOW TO GET EXPECTATION BOUND? (I AM STILL WRITING)


Now, we look at this again. We want to achieve

Given ∑_i=1^mv_iv_i^⊤= and _1,_2,_3,⋯,_k.
Let S=∑_imax_jv_i^⊤_j^-1v_i and β=min_iλ_min_i.
The algorithm 𝚂𝚙𝚊𝚛𝚜𝚒𝚏𝚢(v_i) outputs α_i such
that (α)≤ c_1β S. Also, the matrix Δ=∑_i=1^mα_iv_iv_i^⊤
satisfies

c_2^-1β≼Δ≼ c_2β

Δ≼_i for all i=1,2,⋯,k.

————————

Okay, every term here is important...

———————–

I believe if I can handle one term, I can handle all terms

In that case we have

——————-

Given ∑_i=1^mv_iv_i^⊤= and . The algorithm
outputs the matrix Δ=∑_i=1^mα_iv_iv_i^⊤
such that


  * (α)≤ c_1(λ_min)^-1

  * c_2^-1(λ_min)≼Δ≼ c_2(λ_min)

  * Δ≼

——————————

“



§ NEW ANALYSIS


We look at the following problem without loss of generality: given ∑_i=1^m v_iv_i^⊤=𝐈 for m=Θ(nlog n),  find scalers {α_i}_i=1^m with Θ( n/log^cn) non-zeros, for some constant c, such that the matrix Δ=∑_i=1^m α_i v_iv_i^⊤ satisfies the following conditions:


  * (1-ε)ε·/Θ(polylog n)≼𝔼Δ≼ (1+ε)ε·/Θ(polylog n);

  * Δ≼ε.


We sample every vector v_i with probability 1/log^cn, and define

    Δ=∑_v_iε· v_iv_i^⊤.


It is straightforward to see that

    Δ≼ε∑_i=1^m v_iv_i^⊤=ε·.

Moreover,

    𝔼Δ=∑_i=1^m 1/log^c n·ε v_iv_i^⊤ = ε/log^cn·.






§ NEW SECTION 1


In this section we present the One-sided BSS problem, and show how such problem can be solved by SDP.

 Let  {M_i}_i=1^m be m matrices,  such that M_i≽0 and ∑_i=1^m M_i=I. The goal is to find {α_i}_i=1^m such that the matrix A=∑_i=1^m α_i M_i satisfies the following conditions:


  * A≼ (1+O(1))I;

  * (A)=nnz(α).



Our algorithm proceeds for T iterations, and there is a parameter u_j associated with every iteration j. The formal description of our algorithm is as follows:




 Let A=∑_i=1^mα_i  M_i be the output of 𝖲𝗂𝗆𝗉𝗅𝖾𝖡𝖲𝖲({M_i}, T) with ∑_i=1^m M_i=I. Then, the following two statements hold:


  * nnz(α)≤(A)=T;

  * A≼(u_0+ 4/1-4/u_0T/n)I. In particular, A≼(1+√(T/n))^2 I, if we set u_0=4(1+√(T/n));

  * Every sampled matrix Δ in round j satisfies Δ∙ (u_jI-A_j-1)^-2≤ 4/n·(u_jI-A_j-1)^-2 and Δ≼ 4/u_0· (u_jI-A_j-1) with probability at least 1/2.






By the algorithm description, the resulting matrix A consists of T sampled  matrices, and the trace of every added matrix Δ equals to 1,  hence (A)=T. Since there are at most T different sampled matrices, it holds that nnz(α)≤(A)=T.

For proving the second statement, we define a function Ψ_j=(u_jI-A_j)^-1. By the description of the algorithm, we know that every picked matrix Δ in iteration j satisfies Δ≼ 4/u_0(u_jI-A_j-1). By the Woodbury matrix identity and the assumption that Δ≼ 4/u_0· (u_j-A_j-1), it holds that

    Ψ_j    = (u_jI-A_j)^-1
        =  (u_jI- A_j-1 -Δ)^-1
       ≤(u_jI- A_j-1)^-1 +1/1-4/u_0Δ∙(u_jI-A_j-1)^-2
       ≤(u_j-1I- A_j-1)^-1 + 1/1-4/u_0Δ∙(u_jI-A_j-1)^-2  - 4/1-4/u_01/n(u_j I - A_j-1)^-2,

where the last inequality holds from the fact that the following function

    f(t)≜((u_j-1 + t·4/1-4/u_01/n)I - A_j-1)^-1

is convex and

    df(t)/dt|_t=1≥ f(1)-f(0) = (u_jI-A_j-1)^-1 - (u_j-1I-A_j-1)^-1.

Now we apply the assumption that Δ∙ (u_jI-A_j-1)^-2≤ 4/n·(u_jI-A_j-1)^-2, and obtain

    Ψ_j ≤(u_j-1I- A_j-1)^-1 + 1/1-4/u_0·4/n(u_jI-A_j-1)^-2  - 4/1-4/u_01/n(u_j I - A_j-1)^-2 = Ψ_j-1.

Hence, by induction we have that Ψ_j ≤Ψ_0=n/u_0. This implies that, for any iteration j, it holds that
λ_min(u_jI-A_j) ≥ u_0/n, and it holds for the output matrix A after T iterations that

    λ_max(A)≤ u_0+ 4/1-4/u_0T/n-u_0/n≤ u_0+ 4/1-4/u_0T/n,

which implies that

    A≼(u_0 + 4/1-4/u_0T/n) I.

By setting u_0=4(1+√(T/n)), we have A≼(1+√(T/n))^2 I by direct calculation.

Now for the third statement. Let j≥ 0 be any fixed iteration. Notice that

    Δ= ∑_i=1^m p_i Δ_i = ∑_i=1^m (M_i)/nM_i/(M_i) = 1/n· I,

and Δ∙ (u_jI-A_j-1)^-2 = 1/n· (u_jI-A_j-1)^-2.
By Marokv's Inequality, it holds with probability at least 3/4 that Δ∙ (u_jI-A_j-1)^-2≤4/n· (u_jI-A_j-1)^-2. Also, notice that

    (Δ(u_jI-A_j-1)^-1)=1/n·(u_jI-A_j-1)^-1≤1/nΨ_j≤1/u_0,

and the condition Δ≼ 4/u_0·(u_jI-A_j-1) holds if and only if (Δ· (u_jI-A_j-1)^-1)≤ 4/u_0. By Markov's inequality a random sampled matrix Δ satisfies  Δ≼ 4/u_0·(u_jI-A_j-1)  with probability at least 3/4. Hence, every sampled matrix satisfies the required two conditions with probability at least 1/2.



The lemma above shows that, by sampling O(T) matrices with probability p_i, with high probability we can construct a matrix A from sampled matrices with desired properties. The following theorem shows that,  our final matrix A can be constructed by a semidefinite programming.

Let S be a random O(T) random coordinates sampled from [m]. Let A=∑_i∈ Sα_i M_i be the solution of the following semidefinite programming

    min_α≥ 0( ∑_i∈ Sα_i M_i )      ∑_i∈ Sα_iM_i≼ 4(
    1+√(T/n))^2I.

Then, with high probability in T we have that A≤ 4T.
?Sample each coordinator wst p_i


to be added





