On Solving Linear Systems in Sublinear Time Alexandr Andoni, - - PowerPoint PPT Presentation

On Solving Linear Systems in Sublinear Time

Alexandr Andoni, Columbia University Robert Krauthgamer, Weizmann Institute Yosef Pogrow, Weizmann Institute  Google WOLA 2019

Solving Linear Systems



Input: 𝐵 ∈ ℝ𝑜×𝑜 and 𝑐 ∈ ℝ𝑜



Output: vector 𝑦 that solves 𝐵𝑦 = 𝑐



Many algorithms, different variants:



Matrix 𝐵 is sparse, Laplacian, PSD etc.



Bounded precision (solution 𝑦 is approximate) vs. exact arithmetic



Significant progress: Linear system in Laplacian matrix 𝑀𝐻 can be solved approximately in near-linear time ෨ 𝑃(nnz 𝑀𝐻 ⋅ log

1 𝜗) [Spielman-

Teng’04, …, Cohen-Kyng-Miller-Pachocky-Peng-Rao-Xu’14]

On Solving Linear Systems in Sublinear Time

Our focus: Sublinear running time

2

Sublinear-Time Solver



Input: 𝐵 ∈ ℝ𝑜×𝑜, 𝑐 ∈ ℝ𝑜 (also 𝜗 > 0) and 𝑗 ∈ [𝑜]



Output: approximate coordinate ො 𝑦𝑗 from (any) solution 𝑦∗ to 𝐵𝑦 = 𝑐



Accuracy bound ො 𝑦 − 𝑦∗ ∞ ≤ 𝜗 𝑦∗ ∞



Formal requirement: There is a solution 𝑦∗ to the system, such that ∀𝑗 ∈ 𝑜 , Pr ො 𝑦𝑗 − 𝑦𝑗

∗ ≤ 𝜗 𝑦∗ ∞ ≥ 3 4



Follows framework of Local Computation Algorithms (LCA), previously used for graph problems [Rubinfeld-Tamir-Vardi-Xie’10]

On Solving Linear Systems in Sublinear Time 3

Motivation



Fast quantum algorithms for solving linear systems and for machine learning problems [Harrow-Hassidim-Lloyd’09, …]



Can we match their performance classically?



Recent success story: quantum  classical algorithm [Tang’18]



New direction in sublinear-time algorithms



“Local” computation in numerical problems



Compare computational models (representation, preprocessing), accuracy guarantees, input families (e.g., Laplacian vs. PSD)



Known quantum algorithms have modeling requirements (e.g., quantum encoding of 𝑐)

On Solving Linear Systems in Sublinear Time 4

Algorithm for Laplacians



Informally: Can solve Laplacian systems of bounded-degree expander in polylog(n) time



Key limitations: sparsity and condition number



Notation:



𝑀𝐻 = 𝐸 − 𝐵 is the Laplacian matrix of graph 𝐻



𝑀𝐻

+ is its Moore-Penrose pseudo-inverse 

Theorem 1: Suppose the input is a 𝑒-regular 𝑜-vertex graph 𝐻, together with its condition number 𝜆 > 0, 𝑐 ∈ ℝ𝑜, 𝑣 ∈ 𝑜 and 𝜗 > 0. Our algorithm computes ො 𝑦𝑣 ∈ ℝ such that for 𝑦∗ = 𝑀𝐻

+𝑐,

∀𝑣 ∈ 𝑜 , Pr ො 𝑦𝑣 − 𝑦𝑣

∗ ≤ 𝜗 𝑦∗ ∞ ≥ 3 4,

and runs in time ෨ 𝑃(𝑒𝜗−2𝑡3) for 𝑡 = ෨ 𝑃(𝜆 log 𝑜).

On Solving Linear Systems in Sublinear Time

More inputs? Faster?

5

Some Extensions



Can replace 𝑜 with 𝑐 0



Example: Effective resistance can be approximate (in expanders) in constant running time! 𝑆eff(𝑣, 𝑤) = 𝑓𝑣 − 𝑓𝑤 𝑈𝑀𝐻

+(𝑓𝑣 − 𝑓𝑤) 

Improved running time if



Graph 𝐻 is preprocessed



One can sample a neighbor in 𝐻, or



Extends to Symmetric Diagonally Dominant (SDD) matrix 𝑇



𝜆 is condition number of 𝐸−1/2𝑇𝐸−1/2

On Solving Linear Systems in Sublinear Time 6

Lower Bound for PSD Systems



Informally: Solving “similar” PSD systems requires polynomial time



Similar = bounded condition number and sparsity



Even if the matrix can be preprocessed



Theorem 2: For certain invertible PSD matrices 𝑇, with bounded sparsity 𝑒 and condition number 𝜆, every randomized algorithm must query 𝑜Ω(1/𝑒2) coordinates of the input 𝑐.



Here, the output is ො 𝑦𝑣 ∈ ℝ for a fixed 𝑣 ∈ 𝑜 , required to satisfy ∀𝑣 ∈ 𝑜 , Pr ො 𝑦𝑣 − 𝑦𝑣

∗ ≤ 1 5 𝑦∗ ∞ ≥ 3 4,

for 𝑦∗ = 𝑇−1𝑐.



In particular, 𝑇 may be preprocessed

On Solving Linear Systems in Sublinear Time 7

Dependence on Condition Number



Informally: Quadratic dependence on 𝜆 is necessary



Our algorithmic bound ෩ O(𝜆3) is near-optimal, esp. when matrix 𝑇 can be preprocessed



Theorem 3: For certain graphs 𝐻 of maximum degree 4 and any condition number 𝜆 > 0, every randomized algorithm (for 𝑀𝐻) with accuracy 𝜗 =

1 log 𝑜 must probe ෩

Ω(𝜆2) coordinates of the input 𝑐.



Again, the output is ො 𝑦𝑣 ∈ ℝ for a fixed 𝑣 ∈ 𝑜 , required to satisfy ∀𝑣 ∈ 𝑜 , Pr ො 𝑦𝑣 − 𝑦𝑣

∗ ≤ 1 log 𝑜 𝑦∗ ∞ ≥ 3 4,

for 𝑦∗ = 𝑀𝐻

+𝑐.



In particular, 𝐻 may be preprocessed

On Solving Linear Systems in Sublinear Time 8

Algorithmic Techniques



Famous Monte-Carlo method of von Neumann and Ulam: Write matrix inverse by power series ∀ 𝑌 < 1, 𝐽 − 𝑌 −1 = σ𝑢≥0 𝑌𝑢 then estimate it by random walks (in 𝑌) with unbiased expectation



Inverting a Laplacian 𝑀𝐻 = 𝑒𝐽 − 𝐵 corresponds to summing walks in 𝐻



For us: view 𝑓𝑣

𝑈 σ𝑢≥0 𝐵𝑢𝑐 as sum over all walks, estimate it by sampling

(random walks)



Need to control: number of walks and their length



Large powers 𝑢 > 𝑢∗ contribute relatively little (by condition number)



Estimate truncated series (𝑢 ≤ 𝑢∗) by short random walks (by Chebyshev’s inequality)

On Solving Linear Systems in Sublinear Time 9

Related Work – All Algorithmic



Similar techniques were used before in related contexts but under different assumptions, models and analyses:



Probabilistic log-space algorithms for approximating 𝑀𝐻

+ [Doron-Le Gall-

Ta-Shma’17]



Asks for entire matrix, uses many long random walks (independent of 𝜆)



Local solver for Laplacian systems with boundary conditions [Chung- Simpson’15]



Solver relies on a different power series and random walks



Local solver for PSD systems [Shyamkumar-Banerjee-Lofgren’16]



Polynomial time nnz 𝑇 2/3 under assumptions like bounded matrix norm and random 𝑣 ∈ 𝑜



Local solver for Pagerank [Bressan-Peserico-Pretto’18, Borgs-Brautbar- Chayes-Teng’14]



Polynomial time O(𝑜2/3) and O( nd 1/2) for certain matrices (non-symmetric but by definition are diagonally-dominant)

On Solving Linear Systems in Sublinear Time 10

Lower Bound Techniques



PSD lower bound: Take Laplacian of 2𝑒-regular expander but with:



high girth,



edges signed ±1 at random, and



𝑃( 𝑒) on the diagonal (PSD but not Laplacian)



The graph looks like a tree locally



Up to radius Θ log 𝑜 around 𝑣



Set 𝑐𝑥 = ±1 for 𝑥 at distance 𝑠, and 0 otherwise



Signs have small bias 𝜀 ≈ 𝑒−𝑠/2



Recovering it requires reading Ω(𝜀−2) entries



Using inversion formula, 𝑦𝑣 ≈ average of 𝑐𝑥‘s



Condition number lower bound: Take two 3-regular expanders connected by a matching of size 𝑜/𝜆



Let 𝑐𝑥 = ±1 with slight bias inside each expander

On Solving Linear Systems in Sublinear Time

𝑠 𝑐𝑥 = ±1

𝑐𝑥 = 0 𝑐𝑥 = 0

𝑣

11

Further Questions



Accuracy guarantees



Different norms?



Condition number of 𝑇 instead of 𝐸−1/2𝑇𝐸−1/2?



Other representations (input/output models)?



Access the input 𝑐 via random sampling?



Sample from the output 𝑦?



Other numerical problems?

On Solving Linear Systems in Sublinear Time

Thank You!

12