PaperlayerView exact revision on arXiv

Source-linked AI summary

Randomized low-rank approximation of monotone matrix functions

David Persson, Daniel Kressner

arXiv:2209.11023v3math.NAcs.NA

TL;DR

Randomized low-rank approximation of f(A) can be costly because it requires products with f(A), rather than only with A. funNyström derives an approximation from a Nyström approximation of A, and the paper reports fewer products with A without sacrificing accuracy or reliability.

  • Problem

    Randomized low-rank approximation of f(A) requires products with f(A), which can be significantly more expensive than products with A.

    Source: §1 · p. 2

  • Method

    funNyström approximates f(A) by applying f to a Nyström approximation of A, using products with A rather than f(A).

    Source: §2 · pp. 5–6

  • Results

    Numerical experiments report that funNyström needs fewer products with A than existing methods for low-rank approximation or trace estimation, without sacrificing accuracy or reliability.

    Source: §6 · p. 26

  • Takeaways & Limitations

    Low-rank approximation can also benefit estimation of the effective dimension and diagonal entries of f_μ(A), called Ridge leverage scores.

    Source: Effective dimension · p. 4

    The stated probabilistic error bounds assume f is operator monotone and f(0) = 0.

    Source: §1 · p. 3, §3 · p. 6

Abstract

from arXiv · show

This work is concerned with computing low-rank approximations of a matrix function $f(A)$ for a large symmetric positive semi-definite matrix $A$, a task that arises in, e.g., statistical learning and inverse problems. The application of popular randomized methods, such as the randomized singular value decomposition or the Nystr\"om approximation, to $f(A)$ requires multiplying $f(A)$ with a few random vectors. A significant disadvantage of such an approach, matrix-vector products with $f(A)$ are considerably more expensive than matrix-vector products with $A$, even when carried out only approximately via, e.g., the Lanczos method. In this work, we present and analyze funNystr\"om, a simple and inexpensive method that constructs a low-rank approximation of $f(A)$ directly from a Nystr\"om approximation of $A$, completely bypassing the need for matrix-vector products with $f(A)$. It is sensible to use funNystr\"om whenever $f$ is monotone and satisfies $f(0) = 0$. Under the stronger assumption that $f$ is operator monotone, which includes the matrix square root $A^{1/2}$ and the matrix logarithm $\log(I+A)$, we derive probabilistic bounds for the error in the Frobenius, nuclear, and operator norms. These bounds confirm the numerical observation that funNystr\"om tends to return an approximation that compares well with the best low-rank approximation of $f(A)$. Furthermore, compared to existing methods, funNystr\"om requires significantly fewer matrix-vector products with $A$ to obtain a low-rank approximation of $f(A)$, without sacrificing accuracy or reliability. Our method is also of interest when estimating quantities associated with $f(A)$, such as the trace or the diagonal entries of $f(A)$. In particular, we propose and analyze funNystr\"om++, a combination of funNystr\"om with the recently developed Hutch++ method for trace estimation.

1 Introduction

Randomized low-rank methods applied to f(A) require matrix-vector products with f(A), which can be more expensive than products with A. funNyström instead builds an approximation from a Nyström approximation of A, with theoretical bounds and experiments supporting its accuracy and lower product cost.

1 Introduction

  • Randomized low-rank approximation of f(A) requires products with f(A), which can be significantly more expensive than products with A.Lanczos may converge slowly for difficult functions, while rational Krylov methods require shifted linear solves with A.Source: §1 · p. 2

Matrix square root

  • Matrix square roots support sampling from Gaussian distributions, and a low-rank approximation of A^(1/2) enables cheaper approximate sampling.For ψ = μ + Gω with A = GG^T, one possible choice is G = A^(1/2).Source: Matrix square root · p. 3

Effective dimension

  • The effective dimension and diagonal entries of f_μ(A), called Ridge leverage scores, are quantities for which low-rank approximation can benefit estimation.The effective dimension appears in kernel learning and inverse problems and can be estimated using trace estimation.Source: Effective dimension · p. 4

1.1 Contributions

  • funNyström constructs f(A-hat) from a Nyström approximation of A, avoiding matrix-vector products with f(A).It applies when f is increasing and f(0) = 0; it can also be single-pass, accessing A's entries once.Source: §1.1 · p. 4
  • For operator-monotone f with f(0) = 0, the paper derives probabilistic error bounds in the Frobenius, nuclear, and operator norms.The bounds predict accuracy comparable to randomized low-rank approximation applied directly to f(A).Source: §1 · p. 3, §1.1 · p. 4
  • Numerical experiments report fewer products with A for low-rank approximation or trace estimation, without sacrificing accuracy or reliability.The paper also describes using funNyström to improve Hutch++ trace estimation for matrix functions.Source: §1.1 · pp. 4–5

2 funNystr¨om: Randomized low-rank approximation of matrix functions

funNyström forms a low-rank approximation of f(A) by applying f to a Nyström approximation of A. The construction uses matrix-vector products with A, and at q = 1 requires a single pass over A.

2 funNyström: Randomized low-rank approximation of matrix functions

  • At q = 1, Algorithm 1 requires a single pass over A, and its matrix-vector products with A can be carried out entirely in parallel.QR decompositions at each subspace iteration improve numerical stability, and small eigenvalues may be truncated to mitigate ill-conditioning.Source: §2 · p. 5
  • A variant uses randomized SVD of A to construct a truncated spectral decomposition of Q(Q^T A Q)Q^T.This variant first computes a basis for range(A^qΩ); the passage states that it requires (q + 1)k matrix-vector products with A.Source: §2 · p. 5
  • funNyström returns the spectral decomposition of f(Â_q,k) from a Nyström approximation of A, for increasing f with f(0) = 0.The Nyström approximation is formed using a Gaussian sketch and subspace iterations; the method computes a spectral decomposition of the approximation before applying f.Source: §2 · pp. 5–6
  • For the same number of matrix-vector products, the Nyström approximation can use an additional subspace iteration and is often significantly better than randomized SVD for SPSD A.The passage therefore notes that Nyström is often preferred in this setting.Source: §2 · p. 6

3 Error bounds for operator monotone functions

The section establishes general relations between a Nyström approximation of A and the corresponding approximation of f(A). Its bounds assume an SPSD matrix and an operator-monotone function that vanishes at zero.

3 Error bounds for operator monotone functions

  • The main results assume an SPSD A, a rank-(k+p) Nyström approximation with q ≥ 1, and operator-monotone f satisfying f(0) = 0.The setting orders A’s eigenvalues and assumes λ_k > 0 for convenience; if λ_k = 0, the approximation is exact almost surely.Source: §3 · pp. 6–7

3.1 General results

  • Operator monotonicity and the Nyström construction imply that f(A) is approximated from below.The matrix approximation itself is also bounded in norm by A’s matrix function, while a further lower bound relates it to the projected matrix function.Source: §3.1 · pp. 7–8
  • The nuclear-norm error simplifies to the trace error because f(A) dominates the matrix function of the Nyström approximation.This relation connects approximation of the full matrix function with approximation of tr(f(A)).Source: §3.1 · p. 8
  • Under the setting’s assumptions, f(x)/x is decreasing, providing a bound relating errors for the matrix and its matrix function.The section also identifies Frobenius, nuclear, and operator norms as relevant error measures for applications.Source: §3 · p. 6, §3.1 · p. 8

3.2 Frobenius norm error bounds

For operator-monotone f with f(0) = 0, the section develops Frobenius-norm error bounds for funNyström. It first relates the error to projection errors, then derives probabilistic bounds and sharper bounds for the matrix square root.

3.2.2 Probabilistic bounds

  • Theorem 3.5 is stated for k,p > 4, q ≥ 2, and any u,t ≥ 1; the cited proof passage gives a probability lower bound of 1 − 2t^(−p) − e^(−u²/2).Source: §3.2 · p. 9, §3.2.2 · p. 12

3.2 Frobenius norm error bounds

  • For q ≥ 2, funNyström’s error bounds are comparable to bounds for applying randomized SVD directly to f(A), despite that approach being much more expensive.The bounds also allow Algorithm 1 to replace randomized SVD in Hutch++ for trace estimation; the general bounds here require q ≥ 2.Source: §3.2 · p. 9

3.2.1 Structural bound

  • The structural bound sandwiches funNyström’s approximation error between projection errors of f(A) onto ranges generated by A^qΩ and A^(q−1/2)Ω.Consequently, its error never exceeds that of randomized SVD using the same number of matrix-vector products with A.Source: §3.2.1 · p. 10

3.2.3 Improved bounds for the matrix square root

  • For the matrix square root, the improved bound applies already for q ≥ 1 and improves the exponent of γ and the constants relative to the structural bound.Expectation and deviation bounds can then be derived using Lemma 3.8.Source: §3.2.3 · p. 13

3.3 Nuclear norm error bounds

The section derives nuclear-norm error bounds for funNyström, using a structural comparison and probabilistic results. It also notes improved bounds for the matrix square root.

3.3 Nuclear norm error bounds

  • Theorems 3.9 and 3.10 provide nuclear-norm error bounds, which also bound the error in approximating tr(f(A)) by tr(f(Â_q,k)).This follows because f(A) dominates f(Â_q,k), so the nuclear norm of their difference simplifies to the trace difference.Source: §3.3 · p. 13

3.3.1 Structural bound

  • The structural nuclear-norm error is bounded by the trace of f(A) minus the trace of its projection onto the sampled subspace.The argument applies the structural comparison to the square root of f, which is also operator monotone.Source: §3.3.1 · p. 14

3.3.2 Probabilistic bounds

  • Lemma 3.8 is used to turn the structural bound into the probabilistic bounds of Theorems 3.9 and 3.10.The deviation bound is stated to hold with probability at least 1 − 2t^(−p) − e^(−u²/2).Source: §3.3.2 · pp. 14–15

3.3.3 Improved bounds for the matrix square-root

  • For f(x) = √x, the section states that the exponent on γ can be improved.Expectation and deviation bounds are then obtained using the cited bounds in the preceding analysis.Source: §3.3.3 · p. 15

3.4 Operator norm error bounds

This section develops operator-norm error bounds for the approximation, beginning with a structural bound and then deriving probabilistic bounds. It also gives improved bounds for the matrix square root.

3.4 Operator norm error bounds

  • The section establishes operator-norm error bounds for the approximation under Setting 3.1.Theorem 3.12 assumes k, p ≥ 2 and q ≥ 1; the supplied passage states an expectation bound.Source: §3.4 · p. 15

3.4.1 Structural bound

  • A structural bound holds in any unitarily invariant norm when rank(Ω1) = k.The operator-norm corollary simplifies this structural result.Source: §3.4.1 · p. 15

3.4.2 Probabilistic bounds

  • The expectation bound is obtained from the structural result using a Gaussian-matrix bound and Jensen’s inequality.The passage notes that deviation bounds can also be obtained from a deviation bound on ∥Λ1/2…∥2.Source: §3.4 · p. 15, §3.4.2 · p. 16

3.4.3 Improved bounds for the matrix square-root

  • For f(x) = √x, the authors derive improved operator-norm bounds using the structural corollary.The passage states that an expectation bound follows using a result valid for k, p ≥ 2.Source: §3.4.3 · p. 16

4 Numerical Experiments

The experiments test funNyström on synthetic and application-derived matrices, comparing its accuracy and matrix-vector-product cost with alternatives. They report strong accuracy at substantially lower cost, including an O(N) speed-up in a matrix-function multiplication task.

4 Numerical Experiments

  • The experiments compare funNyström with direct Nyström approximation of f(A), whose products with f(A) may be approximated using Lanczos iterations.Approximation quality is measured by relative error in a selected norm.Source: §4 · p. 17

4.1.1 Synthetic matrices

  • The synthetic test matrices use prescribed algebraic and exponential eigenvalue decays, with n = 5000 unless otherwise specified.Their diagonal spectra are parameterized by s, c, and γ, then combined with an orthogonal matrix.Source: §4.1.1 · p. 17

4.1.2 Gaussian process covariance kernels

  • The Gaussian-process tests use squared-exponential and Matérn covariance kernels; estimating tr(log(I + A)) is relevant to Bayesian optimization and maximum-likelihood estimation.The experiments generate 5000 independent standard-normal data points.Source: §4.1.2 · pp. 17–18

4.1.3 Bayesian inverse problem

  • The Bayesian inverse-problem test links tr(log(I + Apde)) to expected information gain, while fine grids require implicit products with Apde.The matrix arises from a discretized parabolic PDE and measurements of its solution.Source: §4.1.3 · p. 18

4.2 Comparing number of mvps

  • Algorithm 1 needs fewer matrix-vector products with A, by up to three orders of magnitude in the q = 1 experiments.In practice, the comparison method’s Lanczos iteration count is selected heuristically because its exact reference approximation is unavailable.Source: §4.2 · p. 19, Figure 1 · p. 20

4.3 Comparing accuracy

  • Algorithm 1 is never less accurate than the more expensive Nyström approximation of f(A), and is sometimes more accurate.For q = 1, the error bounds recover empirical error within a factor of 10; for q = 2, observed performance can substantially exceed the bounds’ predictions.Source: §4.3 · p. 19, Figure 2 · p. 21

4.4 Fast computation of mvps

  • At similar relative error of approximately 0.01, Figure 3 confirms an O(N) speed-up for approximating N products with A1/2.The experiment uses k = 14, q = 1, and p = 21, with N ranging from 10 to 100.Source: §4.4 · pp. 19, 22, Figure 3 · p. 22

5 Application to trace estimation

The paper uses funNyström low-rank approximations to estimate traces, deriving probabilistic error bounds and comparing them with prior bounds. It also combines funNyström with a stochastic correction in funNyström++ for trace estimation.

5.1 Trace estimation via low-rank approximation

  • Approximating tr(f(A)) by the trace of the low-rank transformed Nyström approximation has probabilistic error bounds from Theorems 3.9 and 3.10.The trace error equals the nuclear-norm error because f(A) ⪰ f(Âq,k+p).Source: §5.1 · pp. 22–23
  • For log(1 + x), the authors report tighter trace-estimation bounds than, and their numerical comparison favors funNyström under the stated matrix-vector-product budget.With budget m, Algorithm 1 can choose k + p = m, whereas the method from can choose k + p = m/2.Source: §5.1 · p. 23

5.2 funHutch++

  • funNyström++ combines a low-rank approximation with a stochastic trace correction, avoiding products with f(A) in its first stage.The second stage still requires or approximates such products; Algorithm 2 extends Hutch++ theory for operator-monotone f with f(0) = 0.Source: §5.2 · pp. 24–25
  • In the experiments, Algorithm 2 with q = 1 or q = 2 can significantly outperform Nyström++ at matched matrix-vector-product counts.The reported q = 1 performance is not covered by Theorem 5.1.Source: §5.2 · pp. 25–26, Figure 5 · p. 26

6 Conclusion

The conclusion presents funNyström as a simple randomized low-rank method for f(A), returning f of a Nyström approximation of A and requiring products only with A. The authors report that it can also cheaply estimate associated quantities such as tr(f(A)).

6 Conclusion

  • For operator-monotone f, the paper derives probabilistic bounds across Frobenius, nuclear, and operator norms, and reports fewer products with A without sacrificing accuracy or reliability.The conclusion describes the method as a potential method of choice when f is operator monotone and f(0) = 0.Source: §6 · p. 26

A Proof of Lemma 3.13

The proof expresses the Nyström residual as a projected form of A, then uses matrix inequalities and auxiliary results to establish the structural bound.

A Proof of Lemma 3.13

  • An auxiliary concavity result for partitioned positive-semidefinite matrices is invoked in proving the structural bound.The cited lemma assumes f is concave.Source: A Proof of Lemma 3.13 · pp. 30–31
  • A − Â_q,k+p is written as A^1/2(I − P_Y)A^1/2, with Y defined from A and Ω.The proof also expresses the residual in the eigenbasis of A using the transformed matrix Ỹ.Source: A Proof of Lemma 3.13 · p. 31
  • The proof compares the residual through f using Lemma 3.1 and a bound involving f(A − Â_q,k+p).This reduces the target norm to a function of Λ^1/2(I − P_Ỹ)Λ^1/2.Source: A Proof of Lemma 3.13 · p. 31
  • A range inclusion yields positive-semidefinite order between the projected residuals, supporting the subsequent norm bound.The proof then exploits the 2×2 block structure of I − P_Z and an auxiliary result.Source: A Proof of Lemma 3.13 · p. 31

B Structural bounds for general Schatten norms

For general Schatten norms, the structural bound uses a finite right derivative of f at zero. The proof treats q = 1 directly and handles q ≥ 2 by following the earlier argument.

B Structural bounds for general Schatten norms

Related papers

Loading 2209.11023v3…