A generalization of the randomized singular value decomposition
Authors: Nicolas Boulle
Alex Townsend
Venue: ICLR 2022 Poster (2022)
Keywords: Low rank approximation
Randomized SVD
Hilbert--Schmidt operators
Gaussian processes

Abstract: The randomized singular value decomposition (SVD) is a popular and effective algorithm for computing a near-best rank $k$ approximation of a matrix $A$ using matrix-vector products with standard Gaussian vectors. Here, we generalize the theory of randomized SVD to multivariate Gaussian vectors, allowing one to incorporate prior knowledge of $A$ into the algorithm. This enables us to explore the continuous analogue of the randomized SVD for Hilbert--Schmidt (HS) operators using operator-function products with functions drawn from a Gaussian process (GP). We then construct a new covariance kernel for GPs, based on weighted Jacobi polynomials, which allows us to rapidly sample the GP and control the smoothness of the randomly generated functions. Numerical examples on matrices and HS operators demonstrate the applicability of the algorithm.

Published as a conference paper at ICLR 2022
A GENERALIZATION OF THE RANDOMIZED SINGULAR VALUE DECOMPOSITION
Nicolas Boullé
Alex Townsend Department of Mathematics Cornell University Ithaca, NY 14853, USA
Mathematical Institute University of Oxford
Oxford, OX2 6GG, UK boulle@maths.. ox ac. uk
townsend@correell.co edu
ABSTRACT
The randomized singular value decomposition (SVD) popular and effective algorithm for computing near-best rank & approximation of a matrix A using matrix-vector products with standard Gaussian vectors. Here, we generalize the randomized SVD to multivariate Gaussian vectors, allowing one to incorporate prior knowledge of A into the algorithm. This enables us to explore the contin- uous analogue of the randomized SVD for Srrbert-Schhiitt (HS) operators using operator-function products with functions drawn from Gaussian process (GP). We then construct new covariance kernel for GPs, based on weighted Jacobi polynomials, which allows us to rapidly sample the GP and control the smooth- ness of the randomly generated functions. Numerical examples on matrices and HS operators demonstrate the applicability of the algorithm.
INTRODUCTION
Computing the singular value decomposition (SVD) is a fundamental linear algebra task in machine learning [Paterek] 2007), statistics (Wold et al. 1987). and signal processing (Alter al. 2000 Van
Der Veen et al.) 1993). The SVD of an m x n real matrix A with m 2 is factorization of the form A- UEV where U is an m X m orthogonal matrix of left singular vectors, Σ is an m. diagonal matrix with entries 0 (A) > ... 2 on (A) 2 O, and V is an
orthogonal matrix of right singular vectors (Golub & Van Loan 2013). The SVD plays a central role because truncating it after l terms provides a best rank l approximation to A in the spectral and Frobenius norm (Eckart& & Young 1936; Mirsky 1960). Since computing the SVD of large matrix can be computationally infeasible, there are various alternative algorithms that compute near-best rank k matrix approximations from matrix-vector vector products (Halko et al.) 2011 Martinsson & Tropp 2020) Nakatsukasa 2020 Nyström 1930 Williams & Seeger 2001). The randomized SVD uses matrix- vector products with standard Gaussian random vectors and is one of the most popular algorithms
for constructing low-rank approximation to (alkoo et al. 2011 Martinsson & Tropp 2020).
Currently, the randomized SVD is used and theoretically justified when it uses matrix-vector prod- uets with standard Gaussian random vectors. In this paper, we consider the following generaliza- tions.
Generalization 1. We generalize the randomized SVD when the matrix-vector products are with multivariate Gaussian random vectors. Our theory allows for multivariate Gaussian random input vectors that have general symmetric positive semi-definite covariance matrix. A key novelty of our work is that prior knowledge of A can be exploited to design covariance matrices that achieve lower approximation errors than the randomized SVD with standard Gaussian vectors.
Generalization 2. We generalize the randomized SVD to Hilbert-Schmidt (HS) operators (Hsing & Eubank, (0115). We design a practical algorithm for learning HS operators using random input functions, sampled from Gaussian process (GP). Examples of applications include learning inte- gral kernels such as Green's 's functions associated with linear partial differential equations (Boullé & Townsend 2022 Boullé et al. (021).
The choice of the covariance kernel in the GP is crucial and impacts both the theoretical bounds and numerical results of the randomized SVD. This leads us introduce new covariance kernel

---

Published as a conference paper at ICLR 2022
based weighted Jacobi polynomials for learning HS operators. One of the main advantages of this kernel is that it is directly expressed as Karhunen- Loève expansion (Karhunen 1946 Loeve 1946) that it is faster to sample functions from the associated GP than using a standard squared- exponential kernel. In addition, we show that the smoothness of the functions sampled from a GP with the Jacobi kernel can be controlled as it is related to the decay rate of the kemel''s eigenvalues.
Contributions. We summarize our novel contributions as follows:
1. We provide new theoretical bounds for the randomized SVD for matrices or HS operators when using random input vectors generated from any multivariate Gaussian distribution. This shows when beneficial to use nonstandard Gaussian random vectors in the ran- domized SVD for constructing low-rank approximations.
2. We generalize the randomized SVD to HS operators and provide numerical examples to learn integral kemels.
3. We propose covariance kernel based on weighted Jacobi polynomials and show that one can select the smoothness of the sampled random functions by choosing the decay rate of the kernel eigenvalues.
BACKGROUND: THE RANDOMIZED SVD FOR MATRICES
The randomized SVD computes nnaa-best rank & approximation matrix A. First, one performs the matrix-vector products 1/1 A. Mk+pp = A. where Z1..... Tk+pp are standard Gaus- sian random vectors with identically and independently distributed entries and 2 1 is oversam- pling parameter. Then, one computes the economized QR factorization [У/ ... 3k+ppl = QR, can form QQ' A by computing Q(AQ) via matrix-vector products involving A. The quality of the rank S k approximant QQ" characterized by the following theorem.
before forming the rank < +++ approximant QQ' A. Note that if A is symmetric,
Theorem (Halko et al. (2011) Let A be an X matrix, k Z integer, and choose oversampling parameter Z 4. S E pgx((kpp) is a standard Gaussian random matrix and
QR AΩ is the economized QR decomposition of AΩ, then for all u,t 1.
A - QQ'AIFF (1+4
3k
Σ j=k+1
of(A) +
Vk + p+1 (4+1(A),
(I)
p+1
with failure probability most P eu?
The squared tail of the singular values of A, i.e., E (A), gives the best rank k approxi- mation eиror to A the Frobenius norm. This result shows that the randomized SVD can compute near-best low-rank approximation to A with high probability.
3 GENERALIZED RANDOMIZED SVD FOR MATRICES AND OPERATORS
The theory behind the randomized SVD has been recently extended to nonstandard covariance ma- trices and HS operators (Boullé & Townsend 2022). However, the probability bounds, generalizing Theorem D are not sharp enough to emphasize the improved performance of covariance matrices with prior information over the standard randomized SVD. We provide new bounds for GPs with nonstandard covariance matrices in Theorem 2
An upper bound the expectation is also available in the Appendix (see Proposition 6). While Theorem 2 formulated with matrices, the same result holds for HS operators in infinite dimensions.
For fixed target rank 1 < k < containing the first and last n right singular vectors of A, respectively, and denote by E, R(n--)x(x-B) the diagonal matrix with entries a++1 (A), (A)....,
define V, Bank and V2 E (==((n=)) to be the matrices
(A). ..).)....AA))) We consider a symmetric positive semi-definite covariance matrix K E BRnx, with keth largest eigenvalue X > 0.
Theorem Let be an m x n matrix, 2 1 an integer, and choose an oversampling parameter p 2 4. 1f Z TRa is Gaussian random matrix, where each column is sampled from
2

---

Published as a conference paper at ICLR 2022
multivariate Gaussian distribution with covariance matrix K R=x", and QR AΩ is the economized QR decomposition of AQ, then for all u, 4,t 2
|A QQ'AII F (1+ uty
3k
Bk
(k p)
Σ o3(A):
(2)
++ Ykk
j=k++1
with failure probability most t- (w2-11//// Here, 7 (A denotes the covariance quality factor, and By. = "K(33V=KV2)/(x1 E2ll)) where A is the largest eigenvalue of K.
(((()KV))-)))
The factors "/k and 3k. measuring the quality of the covariance matrix to learn A in Theorem can be respectively bounded (Boullé & Townsend 2022 Lem. 2: Lemma 9 using the eigenvalues A 2 ... 2 An of the covariance matrix K and the singular values of as:
A
A j-k 0)(A) A1
BR
Σ j=++1
E))(A) (A)) j=k+1
7k
j=n-k+1
This shows that the performance of the generalized randomized SVD depends on the decay rate of the sequence { The quantities 7/k and 3g depend on how much prior information of the +, 1...., right singular vectors of A encoded in K. In the ideal situation where these singular vectors are known, then one can define K such that 3g 0 for A6+1 =... An 0. In particular, this highlights that suitably chosen covariance matrix can outperform the randomized SVD with
standard Gaussian vectors (see Section 5.1 for numerical example).
3.1 RANDOMIZED SVD FOR HILBERT-SCHHIIN OPERATORS
We now describe the randomized SVD for learning HS operators (see Algorithm T}. The algorithm implemented in the Chebfun software system (Driscoll et al.) (0044), which is MATLAB package for computing with functions. The Chebfun implementation of the randomized SVD for HS oper- uses Chebfun's is capabilities, which offer continuous analogues of several matrix operations like the QR decomposition and numerical integration. Indeed, the continuous analogue of matrix- vector multiplication AΩ for an HS integral operator F (see Hsing Eubank 2013 for definitions and properties of HS operators), with kernel G DXD-1 R.
(z)))) = - G(I, ())( EED, f L(())
where C Rdd with 2 1, and (D) is the space of square-integrable functions on D.
Algorithm Randomized SVD for HS operators
Input: HS integral operator & with kernel G(z, v), number of samples k > 0
Output: Approximation G2 of G
1: Define a GP covariance kernel K
2: Sample the GP & times to generate a quasimatrix of random functions Ω = [f1... fxl
3: Evaluate the integral operator at Ω Y
((x))
4: Orthonormalize the columns of Y, Q = orth(Y) [q (n-..q
5: Compute an approximation
by evaluating the adjoint of &
6: Initialize GR(z, y) to 0
7: for = do
8: Gz v) + Gx(x, v) + q:(x) So ((2)))
9: end for
The algorithm takes as input an integral operator that we aim to approximate. Note that we focus here learning integral operator, but other HS operators would work similarly. The first step of the randomized SVD for HS operators consists of generating an oo x quasimatrix by sampling GP & times, where lc is the target rank (see Section 4). An set & functions (Townsend & Trefethen 2015), and generalizes the notion of matrix to infinite dimensions. Therefore, each column of Q is object, consisting of polynomial approximation of
quasimatrix is an ordered
3

---

Published as a conference paper at ICLR 2022
smooth random function sampled from the GP in the Chebyshev basis. After evaluating the HS operator at to obtain a quasimatrix Y, we use the QR algorithm (Townsend Trefethen, 2015) to obtain orthonormal basis Q for the range of the columns of Y. Then, the randomized SVD for HS operators requires the left-evaluation of the operator & or, equivalently, the evaluation of its adjoint Ft satisfying:
(f)(()) - G(y, (f)) dy,
x x D.
We evaluate the adjoint of & at each column vector of Q to construct an approximation G3 of G. Finally, the approximation enror between the operator kernel G and the learned kernel Gg can be computed in the -,2-norm, corresponding to the HS norm of the integral operator.
4 COVARIANCE KERNELS
To generate the random input functions fi.....fa f& for the randomized SVD for HS operators, we draw them from GP, denoted by GO K), for a certain covariance kernel K. widely employed covariance kernel is the squared-exponent
function KSE (Rasmussen & Williams/ 2006) given by
(sF y) = exp (-|x - v2?)( y E D,
(3)
where l > O a parameter controlling the length-scale of the GP. This kernel isotropic as it only depends |2: y/. is infinitely differentiable, and its eigenvalues decay supergeometrically to 0. Since the bound in Theorem degrades as the ratio Å A increases for j 2 k + 1, the randomized SVD for learning HS operators prefers covariance kernels with slowly decaying eigenvalues. Our randomized SVD cannot hope to learn HS operators where the range of the operator has rank greater than k, where is such that the kthh eigenvalue of KSE reaches machine precision.
Other popular kernels for GPs include the Matérn kernel (Rasmussen Williams/ 2006 Chapt. 4) and Brownian bridge (Nelsen & Stuart 2021). Prior information on the HS operator can also be enforced through the choice of the covariance kernel. For instance, one can impose the periodicity of the samples by using the following squared-ext -etaured-xppnnntilll periodic kernel:
(-
2 82 sin² ( (²=)):
Kre(F, y)
exp
E. D,
where & > 0 is the length-scale parameter.
4.1 SAMPLE RANDOM FUNCTIONS FROM A GAUSSIAN PROCESS
In finite dimensions, a random vector u ~ (O, K), where K E R a covariance matrix with Cholesky factorization K LL' can be generated from the matrix-vector product Lc. Here, E R" is vector whose entries follow the standard Gaussian distribution. We now detail how this process extends to infinite dimensions with a continuous covariance kernel. Let K be continuous symmetric positive-definite covariance function defined on the domain [a, b] [a, b] C IR with Do. We consider the continuous analogue of the Cholesky factorization to write K as (Townsend & Trefethen 2015)
-00
K(a,y)
Σ (L))) (v) x,y E (a 1=1
where ry is the jth row of Le, which - in Chebfun's terminology- is a lower-triangula triangular quasi- matrix. In practice, we truncate the series after terms, either arbitrarily or when the nihh largest
kernel eigenvalue, An, falls below machine precision. Then, if E RD follows the standard Gaussian distribution, function can be sampled from GP K) as u T.c. That is,
u u<x) Σ gry((), x b)
The continuous Cholesky factorization is implemented Chebfun2 (Townsend & Trefethen 2013). which is the extension of Chebfun for computing with two-dimensional functions. As example,
4

---

Published as a conference paper at ICLR 2022
the polynomial approximation, which is accurate up essentially machine precision, of the squared- exponential covariance kernel KSE with parameter 0.01 on [-1, yields a numerical rank of = 503. The functions sampled from GP(O, KSE) become more oscillatory the length-scale parameter é decreases and hence the numerical rank the kernel increases or, equivalently, the associated eigenvalues sequence (Az) decays slower zero.
4.2 INFLUENCE OF THE KERNEL'S 'S EIGENVALUES AND MERCER'S ' REPRESENTATION
The covariance kernel can also be defined from its Mercer's representation
(xx,))= Σ (1 T, D,
(4)
where { is orthonormal basis of (D) and À 2 2 Thm. 4.6.3). We prefer to construct K directly from Mercer's representations for several reasons:
(Hsing & Eubank 2015
1. One can impose prior knowledge of the kernel of the HS operator on the eigenfunctions of K (such as periodicity or smoothness), One can often generate samples from GP K) efficiently using Equation
4 and 3. One can control the decay rate of the cigenvalues of K,
Hence, the quantity "Ak in the probability bound of Theorem 2 measures the quality of the covariance kernel K in GP(O, K) generate random functions that learn the HS operator F. To minimize 1 1/7h we would like to select the eigenvalues A 2 X 2 ... > 0 of K so that they have the slow- possible decay rate while maintaining 2,1=1
>j oc. One needs { ld to guarantee that w GP(O, K') has finite expected squared 7.--norm, i.e., B|||w|||
25=1 X < cp. The best sequence of eigenvalues we know that satisfies this property is called the Rissanen sequence (Rissa- nen. 1983) and is given by Xj- = Rg == 2 L(j), where
L(j) logo (co) + logg(i), logg (i)
Σ
max
(1) (i),0),
Co
2 log
and log2 (1) log2 o... lo) is the composition of log2(-) i times. Other options for the choice of eigenvalues include any sequence of the form X3 = = for > 1.
4.3 JACOBI COVARIANCE KERNEL
If D 1-1, then a natural choice of orthonormal basis of (D) define the Mercer's repre- sentation of the kernel are weighted Jacobi polynomials (Deheuvels & Martynov 2008 Olver et al.)
2010 That is, for weight function Wa (x) = (1 - x1 2)) a, 11, and any positive
eigenvalue sequence {}i}. consider the Jacobi kernel
k<(a..) (z. v) Jnc
Σ Niti (xx) pla (l)) (v), 2, E,ye[-1,11, i=0
(5)
where p(o..) is the scaled Jacobi polynomial of degree j. The polynomials are normalized such that [| 1/2 p(a.3)
p((03)) [--(11111)) and k/aa..) [[[(__11)]3) In this case, random function can be
sampled as
u(x) Σ VAttic 1/2 p(a
1/2 1/2p(a
(x), 1-1-11
1=0
where oy ~ N(0, for Sjjs..
A desirable property of covariance kernel is to be unbiased towards one spatial direction, i.c., K (z, y) = ((-y.- -z) ) for x, y E [-,, which motivates us to always select a = 3. Moreover, desirable to have the eigenfunctions of K|o to be polynomial so that one generate samples from (P K) efficiently. This leads us to choose a and 3 to be even integers. The choice of
8 0 gives the Legendre kernel (Foster et al. 2020 Habermann 2021). In Section 5 use Equation (5) with a = 3 = 2 (see Figure[] to ensure that functions sampled from the associated GP satisfy homogeneous Dirichlet boundary conditions. In this case, one must have ΣØ= jλح < oo
so that the series of functions in Equation (5) converges uniformly and K (2,2) continuous kernel
5

---

Published as a conference paper at ICLR 2022
(see Appendix B). Under this additional constraint, the best choice of eigenvalues is given by scaled Rissamen sequence: Xy = =Rg/j. for j 2 1 (cf. Section 4.2). Covariance kernels on higher dimensional domains of the form = [-1, 1] for d 2 2, can defined using tensor products of weighted Jacobi polynomials.
4.4 SMOOTHNESS OF FUNCTIONS SAMPLED FROM A GP WITH JACOBI KERNEL.
We now connect the decay rate of the eigenvalues of the Jacobi covariance kernel K (2.2) to the smoothness of the samples from GP(O, K k22-2)). Hence, the Jacobi covariance function allows the control of the decay rate of the eigenvalues { as well as the smoothness of the resulting randomly generated functions. First, Lemma 3 asserts that if the coefficients of an infinite polynomial series have sufficient decay. then the resulting series is smooth with regularity depending on the decay rate. This result can be seen as a reciprocal to (Trefethen 2019) Thm. 7.1).
Lemma Let {py}
P be a family of polynomials such that deg(p;) ( S and max (P)(x)) = 1. If (x) = Σド= (x) with layl S j" for v > 1, then fn converges uniformly to f(x) = Ej=o aypp and is y times continuously difserentiable for uny integer p such that < (1-1)//2
Note that the main application of this lemma occurs when deg(p;) j for all 2 0. The proof of Lemma 3 is deferred the supplementary material. We now state the following theorem about the regularity of functions sampled from GP(0, K (2.2)) which guarantees that if the eigenvalues Jac
are chosen such that A = (11j"*) with v > 3, then f f~ GP(O, (2.22) is almost surely continu- ous. Moreover, faster decay of the eigenvalues of K (2.2) implies higher regularity of the sampled functions, in an almost sure sense.
Inc
Theorem Let { E (IR be positive sequence such that X3 = =O(=") ") for v > 3. If f is sampled from GP(O, (2:22), then f E C=([-1, almost surely for any integer μ < (v (33))/2 32/
This theorem can be seen as a particular case of Driscoll's zero-one law (Driscoll, 1973 which
characterizes the regularity functions samples from GPs (see also Kanagawa et 2018
X, 1/j4
A 1/j3
A R,/j
1.00
1.0
1.00
1.0
1.00
0.75
0.5
0.5
0.75
0.0
0.00
0.0
0.00
0.0
0.50
0.25
0.5
-1.0
0.00
1.0
0.00
1.0
0.00
-1.0
2-
-1
-1-
-2
-1.0
-0.5
0.0
0.5
1.0
-1.0
-0.5
0.0
0.5
1.0
-1.
-0.5
0.0
0.5
1.0
Figure 1: Covariance kernel (222)
constructed using Jacobi polynomials of type (2, 2) with A = 1/ja, 1/j3. and R,/j. where R, the Rissanen sequence (top). The bottom panels illus- trate functions sampled from (P( K-(221 with the different eigenvalue sequences. The series for generating the random functions are truncated to n = 500.
Jac
In Figure I we display the Jacobi kernel type (2 with functions sampled from the correspond-
ing GP. We selected eigenvalue sequences of different decay rates: from the faster 1/j* to the slower Rissanen sequence Rí/j (Section 4.2p. For X = 1/j and A = R,/j. we observe large variation +1, indicating potential discontinuity of the sam-
of the randomly generated functions near ples at these two points as n + oc. This in agreement with Theorem 4 which only guarantees continuity (with probability one) of the randomly generated functions if X 1/j" with 3.

---

Published as a conference paper at ICLR 2022
5 NUMERICAL EXPERIMENTS
5.1 APPROXIMATION OF MATRICES USING NON-STANDARD COVARIANCE FUNCTIONS
The approximation error bound in Theorem 2 depends the eigenvalues of the covariance matrix, which dictates the distribution of the column vectors of the input matrix Ω. Roughly speaking, the more prior knowledge of the matrix A that can be incorporated into the covariance matrix, the better. In this numerical example, we investigate whether the standard randomized SVD, which uses the identity its covariance matrix, can improved by using a different covariance matrix. We then attempt to learn the discretized 2000 x 2000 matrix, i.e., the discrete Green's function, of the inverse of the following differential operator:
Lu d 100 sin(5
[0, 1].
We vary the number of columns (i.e. samples from the GP) in the input matrix from 2000.
0.5
Error Beed approcimation error
Standard covarionce
Standard
Prrior covariance
Best approcimation
Computational time (s)
2.5
2
0.2
0.1
500
1,000
1,500
2,000
500
1,000
1,500
2,000
Number of samples
Number of samples
Figure 2: Left: Ratio between the average randomized SVD approximation error (over 10 runs) of the 2000 2000 matrix of the inverse of the differential operator Lu d/d/dz" 100 (5mrr on [o, 1]. and the best approximation error. The error bars in light colour (blue and red) illustrate one standard deviation. Right: Average computational time of the algorithm (over 10 runs). The cigenvalue decomposition of the covariance matrix has been precomputed before.
In Figure left), we compare the ratios between the relative error the Frobenius norm given by the randomized SVD and the best approximation error, obtained by truncating the SVD of A. The prior covariance matrix K consists of the discretized 2000 x 2000 matrix of the Green's function of the negative Laplace operator Lu -dau u/dx" [0, 1] to incorporate knowledge of the diffusion term in the matrix A. We see that nonstandard covariance matrix leads to a higher approximation accuracy, with a reduction of the eлror by factor of 1.3-1.6 compared the standard randomized SVD. At the same time, the procedure is only 20% slower") on average (Figure 2 (right)) as one can precompute the eigenvalue decomposition of the covariance matrix (see Appendix C of interest to maximize the accuracy of the approximation matrix from a limited number of samples in applications where the sampling time much higher than the numerical linear algebra costs.
5.2 APPLICATIONS OF THE RANDOMIZED SVD FOR HILBERT- HILBERT-SCHMIDT SCHMIDT OPERATORS
We now apply the randomized SVD for HS operators to learn kernels of integral operators. In this first example, the kernel defined as [Twwsendd] 2013)
G(æ, y) = cos(10(2² + v)) (10 + +))) x, y [-1,1]
and is displayed in Figure 3 a). We employ the equared-epponnntia covariance kernel KSE with parameter ( = 0.01 and A 100 samples (see Equation (3) to sample random functions from the associated GP. The learned kernel G₁ is represented on the bottom panel of Figure Bla) a) and has an approximation error around machine precision.
'Timings were performed an Intel Xeon CPU RS-2667 v2 @ 3.30GHz using MATLAB R2020b without explicit parallelization.
7

---

Published as a conference paper at ICLR 2022
(a)
10
(b)
0.50
(c)
1.0
Kernel
y
0.5
0.25
0.5
-1.0
-1
-0.50
1.0
x
57
10
0.50
Learned Kernel
9
0.5
0.25
0.5
-1.0 -1
-0.50
-1
1.0
x
27
Figure 3: Kernels of three HS operators (top) together with the kernels learned by the randomized SVD for HS operators (bottom), using the squared-exponentin covariance kernel KSE with param- eter 8 = 0.01 and one hundred functions sampled from GP(O, KsF).
As a second application of the randomized SVD for HS operators, we learn the kernel G(r,, y) = Ai(-13 for x, y [-1, 1], where Ai is the Airy function (NIS) Chapt. 9) defined by
00
t
Ai(x)
cos
+ rt dt,
I TER.
π
We plot the kernel and its low-rank approximant given by the randomized SVD for HS operators in Figure and obtain an approximation error (measured the L²-norm) of 5.04 10-14 The two kernels have numerical rank equal to 42.
The last example consists learning the HS operator associated with the kernel G(z, y) Jo (1100(zz+t²)) for z, y E [-1, where Jo is the Bessel function of the first kind NNS Chapt. 10) defined as
Jo(x) = _ / cos(a 1) dt, x E R, π
and plotted in Figure c). The rank of this kernel is equal to 91 while its approximation is of rank 89 and the approximation error equal to 4.88 x 10-13 We observe that in the three numerical examples displayed in Figure 3, the difference between the learned and the original kernel is visually not perceptible.
Finally, we evaluate the influence of the choice of covariance kernel and number of samples in Figure Here, we vary the number of samples from = to = 100 and use the randomized SVD for HS operators with four different covariance kernels: the squared-exponential KSE with parameters e = 0.01, 0 0. 1, and the Jacobi kernel ((2.2) with eigenvalues A; = 1/j3, for j 2 In the left panel of Figure 4 represent the eigenvalue ratio / of the four kernels and observe that this quantity falls below machine precision for the squared-exponentn 13 and
kernel with = 1 and 59, respectively. In Figure [ (right), observe that these two
l = 0.1 kernels fail to approximate kernels of high numerical rank. The other two kernels have a much complicated kernels, We then
slower decay of eigenvalues and can capture (or learn) the right panel of Figure 4 that the relative approximation errors obtained using K and KSE are close to the best approximation error given by the squared tail of the singular values of the integral kernel G(x, v). i.e., (E)2441 0/) The overshoot in the error k 100 compared to the machine precision is due to the decay of the cigenvalues of the covariance kernels. Hence, spatial directions associated with small eigenvalues are harder to learn accurately. This issue does not arise in finite dimensions with the standard randomized SVD because the covariance kernel used there
in
8

---

Published as a conference paper at ICLR 2022
10
10
Kia Ay 1/5²
Kon. == 0.01
10
Rse: 6
10-2
Kyu:
Relutive crror
10-3
10-5
Xy
10
10
10
10-11
10 14
10
100
200 300
400
500
20
40
60
80
100
Number of samples
Figure 4: Left: Scaled eigenvalues of the Jacobi covariance kernel K[2.2) with sequence X; 1/ja and seurred-xpponentill kernels KSE with parameters & = 0.01, 0.1, 1, respectively. Right: Average (over 10 runs) relative approximation error in the L²-norm between the Bessel kernel G(xr, 3) Jo(100 + y})) and its low-rank approximation Ga (2,,)) y). obtained from the randomized SVD by sampling the GPs times. The error bars in light colour (blue and red) illustrate one standard deviation and the black line indicates the best approximation values of G.
given by the tail of the singular
is isotropic, i.e., all its eigenvalues are equal to one. However, this choice is no longer possible for learning HS integral operators as the covariance kernel K must be squared-integrable. The relative approximation errors at k = 100 (averaged over 10 runs) using KL (2.2) and KSE with é = 0.01 are Jac Error( K (2.2)) R 2.6 10-11, and Error (sE) 5.7 x 10-13, which gives a ratio of
Error(K]e (2.2) (Estor(KE)) = 45.6.
(6)
However, the square-rodd root of the ratio of the quality of the two kernels for k = 91 is equal to
791 (KSE)/79) (177999(866
(K((22)) (2.2) 117.8,
which is of the same order of magnitude of Equation as predicted by Theorem 2 Equation 7
791 (KsE) x 5.88 10-2 and "99 (k( 4.24
10
both computed using Chebfun.
In conclusion, this section provides numerical insights to motivate the choice of the covariance kernel to learn HS operators. Following Figure 4 kernel with slowly decaying eigenvalues is preferable and yields better approximation errors or higher learning rate with respect to the number of samples, especially when learning kernel with large numerical rank. The optimal choice from a theoretical viewpoint is to select a covariance kernel whose cigenvalues have a decay rate similar the Rissanen sequence (Rissanen) 1983), but other choices may be preferable in practice to ensure smoothness of the sample functions (cf. Section 4.4).
6 CONCLUSIONS
We have explored practical extensions of the randomized SVD to nonstandard covariance matri- and Hilbert-Schmidt operators. This paper motivates computational and algorithmic ap- proaches for preconditioning the covariance kernel based on prior information to compute a low- rank approximation of matrices and impose properties on the learned matrix and random functions from the GP. Numerical experiments demonstrate that covariance matrices with prior knowledge outperform the standard identity matrix used in the literature and lead to near-optimal approxima- tion enrors. In addition, we proposed a covariance kernel based on weighted Jacobi polynomials, which allows the control of the smoothness the random functions generated and may find prac- tical applications in PDE learning (Boullé et al.] 2020 2021) it imposes prior knowledge of Dirichlet boundary conditions. The algorithm presented in this paper is limited to matrices and HS operators and does not extend unbounded operators such as differential operators. Additionally, the theoretical bounds only offer probabilistic guarantees for Gaussian inputs, while sub-Gaussian distributions (Kahane 1960) of the inputs would be closer realistic application settings.

---

Published as a conference paper at ICLR 2022
ACKNOWLEDGMENTS
We thank Joel Tropp for his suggestions leading Lemma and Daniel Kressner for his com- ments. This work was supported by the EPSRC Centre For Doctoral Training in Industrially Fo- cused Mathematical Modelling (EP/LO15803//) in collaboration with Simula Research Laboratory, and the National Science Foundation grants DMS-1818755 DMS-1952752 and DMS- RMS-2045646.
REFERENCES
NIST Digital Library of Mathematical Functions. http://dlmf..ist.gov... Release 1.1.0 of 2020-12-15. F. W. Olver, A. B. Olde Daalhuis, W. Lozier, B. Schneider, F. Boisvert, C. W. Clark,
B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
O. Alter, P. O. Brown, and D. Botstein. Singular value decomposition for gedo expression data processing and modeling. Proc. Narl. Acad. Sci., 97(18):00011-10000 2000.
N. Boullé and A. Townsend. Learning elliptic partial differential equations with randomized linear algebra. Found. Comput. Math., 2022.
N. Boullé, Y. Nakatsukasa, and A. Townsend. Rational neural networks. In Proc. Advances in Neural Information Processing Systems, volume 33, pp. 14243-14253, 2020.
N. Boullé, C. J. Earls, and A. Townsend. Data-driven discovery of physical laws with human- understandable deep learning. arXiv preprint mxxi:2105.002666 2021.
H. Chernoff. A measure of asymptotic efficiency for tests of hypothesis based on the sum of observations. Ann. Math. Stat. pp. 493-507, 1952.
Deheuvels and G. V. Martynov. Karhunen-Loes decomposition of Gaussian process gener- ated by independent pairs of exponential random variables. J. Funct. Anal., 255(9):2333-2394 2008.
M. F. Driscoll. The reproducing kernel Hilbert space structure of the sample paths of a Gaussian process. Z. Walhsscheinlichkeittthherrr 26(4):309-31) 3366 1973.
T. A Driscoll, N. Hale, and L. N. Trefethen. Chebfin Guide. Pafnuty Publications, 2014. URL http:/ http::/wwwt wwww///wwwwccco//////////ww chebfun..cuu httg:/doo//qquii
R. Durrett. Probability: Theory and Examples. Cambridge University Press, Sth edition, 2019.
C. Eckart and G. Young. The approximation of matrix by another of lower rank. Psychometrika, 1(3):211-21) 1-218, 1936.
J. Foster, T. Lyons, and H. Oberhauser. An optimal polynomial approximation of Brownian motion. SIAM J. Numer. Anal., 5933)119, 2020.
G. H. Golub and C. F. Van Loan. Matrix Computations. JHU Press, 4th edition, 2013.
K. Habermann. A semicircle law and decorrelation phenomena for iterated Kolmogorov loops. J. London Math. Soc., 103(2):558-586, 2021.
N. Halko, P.-G. Martinsson, and J. A. Tropp. Finding strueture with randomness: Probabilistic al-
gorithms for constructing approximate matrix decompositions. SIAM Rev. 53(2):1172288, 2011.
T. Hsing and R. Eubank. Theoretical foundations of functional data analysis, with an introduction to linear operators. John Wiley & Sons, 2015.
J.-P. Kahane. Propriétés locales des fonctions à séries de Fourier aléatoires. Stud. Math., 19(1): 1-25, 1960.
M. Kanagawa, P. Hennig, D. Sejdinovic, and B. K. Sriperumbudur. Gaussian processes and kemnel methods: review on connections and equivalences. arXiv preprint 30XXvv/807.02533 2018.
K. Karhunen. Über lineare Methoden in der Wahnschennlccheissrcchnnng Ann. Acad. Science Fenn., Ser. A. I., 37:3-79, 1946.
10

---

Published as a conference paper at ICLR 2022
T. Kato. Perturbation theory for linear operators, volume 132. Springer Science & Business Media, 2013.
F. C. Leone, L. S. Nelson, and R. B. Nottingham. The folded normal distribution. Technometrics, (4):543-550, 1961.
M. Loève. Fonetions aleatoire de second ordre. Rev. Sci., 84:195-206, 1946.
A. A. Markov, On question by D. I. Mendeleev. Zapiski Imp. Akad. Nauk, 62(12):1-24, 1889.
P.-G. Martinsson and A. Tropp. Randomized numerical linear algebra: Foundations and algo- rithms. Acta Numer. 29:403572, 2020.
L. Mirsky. Symmetric gauge functions and unitarily invariant norms. Q. J. Math., 11(1):00-59, 1960.
Mirsky, A trace inequality of John von Neumann. Monatsh. Math., 79(4):003-366, 1975.
Y. Nakatsukasa. Fast and stable randomized low-rank matrix approximation. arXiv preprint arXiv:2009.11390 2020.
N. H. Nelsen and A. M. Stuart. The random feature model for input-output maps between Banach spaces. SIAM J. Sci. Comput., 43(5):43212-3324, 2021.
E. Nyström. Über die praktische Auflösung von integralgleichungen mit Anwendungen auf randwertaufgaben. Acta Math., 144:185-244 1930.
W. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark. NIST Handbook of Mathematical Functions. Cambridge University Press, 2010.
A. Paterek. Improving regularized singular value decomposition for collaborative filtering. Pruc. KDD cup and workshop, pp. 5-8, 2007.
C. Rasmussen and C. Williams. Gaussian processes for machine learning. MIT Press, 2006.
J. Rissanen. A universal prior for integers and estimation by minimum description length. Ann. Stat., 11:416-433 4311 1983.
W. Rudin. Principles of mathematical analysis. McGraw-Hill, 3rd edition, 1976.
A. Townsend. Pretty functions approximated by Chebfun2. https: htttp::/www.co. www.w...ww...occc/m//o/ org/ wxamp//ee/ppproo///Pettttttteeee htmll 2013.
A. Townsend and L. N. Trefethen. An extension of Chebfun to two dimensions. SIAM J. Sci. Comput., 35(6):C495 C518, 2013.
A. Townsend and L. N. Trefethen. Continuous analogues of matrix factorizations. P. Roy. Soc. A. 471(2173):00145555 2015.
L. Trefethen. Approximation Theory and Approximation Practice. SIAM, extended edition, 2019.
A.-J. Van Der Veen, E. Deprettere, and A. L. Swindlehurst. Subspace- based signal analysis using singular value decomposition. Proc. IEEE, 8179):2777-1388 1993.
J. Von Neumann. Some matrix-inequalittis and metrization of matric space. Tomsk Univ. Rev., 1: 286-300, 1937.
C. Williams and M. Seeger. Using the Nyström method to speed up kernel machines. In Proc. Advances in Neural Information Processing Systems, volume 14, pp. 682-688, 2001.
S. Wold, K. Esbensen, and P. Geladi. Principal component analysis. Chemometr. Intell. Lab., 2(1-3): 37-52, 1987.
11

---

Published as a conference paper at ICLR 2022
A RANDOMIZED SVD WITH MULTIVARIATE GAUSSIAN INPUTS
Let m 2 n 2 and A be an m x n real matrix with singular value decomposition = UVV* where U and V are orthonormal matrices, and Σ be an n diagonal matrix with entries o, (A) 2 (A) > 0. For a fixed target rank 2 1, we define Σn and Σn to be the diagonal matrices, which respectively contain the first k singular values of A: 01 (A) 2 > (A),,, and the remaining singular values. Let 11 be the n x k matrix obtained by truncating 11 after columns and V2 the remainder. In this section, K denotes symmetric positive semi-definite n matrix and Ω E Rnx( a Gaussian random matrix with & 2 k independent columns sampled from a multivariate normal distribution with covariance matrix K. Finally, we define S2 = Vi and Ω2 = ViΩ. We first refine (Boullé & Townsend 2022 Lem. 5).
Lemma Let 8 2 1, S Ruxt be a Gaussian random matrix, where each column is sampled from a multivariate Gaussian distribution with covariance matrix K, and T be an & x k matrix. Then,
E|| Σ2 VASTT = Tr I
Proof. Let K - A be the eigenvalue decomposition of K, where QK is orthonormal and A is a diagonal matrix containing the eigenvalues of K in decreasing order: A1 2 ... 2 An 2 O. We note that 52 can be expressed as Ω QK Ak/2//C where is a standard Gaussian matrix. Prop. A.1), which shows that
Let Σg V AKA2/2 1/2, the proof follows from (Halko al. 2011
B|GT = | ITI
Note that one can bound the term Tr (E by A1 A 1... where A is the largest eigenvalue of K (Boullé & Townsend, 2022). While this provides a simple upper bound, it does not demonstrate that the use of a covariance matrix containing prior information on the singular vectors of can outperform the randomized SVD) with standard Gaussian inputs. Then, combining the proof of (Boullé & Townsend 2022, Thm. I) and Lemma 5 we prove the following proposition, which bounds the expected approximation error of the randomized SVD with multivariate Gaussian inputs.
Proposition 6 Let be an m x matrix, & 2 1 an integer, and choose an oversampling parameter P 2 2. If Q (Rxx(++p)) is a Gaussian random matrix, where each column is sampled from multivariate Gaussian distribution with covariance matrix K E gaxx, and QR AQ is the economized QR decomposition of ASt, then,
BR k(k + p)
E I QQ'QQQQAAAT
A S 1+
Σ (A), j=k+1
Tkk
P
where Yh- k/(A1 (v((KVV,, and 3x = ThE [
We remark that for standard Gaussian inputs, we have 7 = 3x = in Proposition and we recover
the average Frobenius eπrоr of the randomized SVD (Hakko et al. 2011 Thm. 10.5) up to factor of (A: + p) due to the non-independence of 52, and 52, in general.
Lemma 7 With the notations introduced at the beginning of the section, for all 2 0, we have st/2
P PI > E(1 5) Tr (EE s (1 + =2/22
Proof. Let Wj be the jth column of Ω for S j S l and uI,..., en be the n columns of the orthonormal matrix V. We first remark that
n-k
1 f = Σin Zj
2j = Σ o (A) Wz)²
2=1
AL=1
where the Z, are i.i because Wj (0, K) are i.i.d. Let d 2 An 2 be the eigenvalues of K with eigenvectors vi,..., Un
E R". For 1 S S l, we have,
W/ Σ V
where 100 ~N(00, 1) are Lid 515 and S < l. Then,
n Z₁
n-k
n
Σ
(0
Σ
oirm (A) 2= Σ X,
1=1
n==1
12

---

Published as a conference paper at ICLR 2022
where the X; independent. Let N 1 I isis n.
> Em oe (A) 1.), then X, XXX for
Let 0 0 < 1//2:==11 7a), bound the moment generating function of Σi= X4 as
-1/2
E
c -
Xi
II B [
II 20%() 11/2 - 1 - 20 1=1
because the Xi/h are independent and follow a chi-squared distribution. The right inequality is obtained by showing by recurrence that, if ff,....... On 2 0 such that Σ 0g S then I a,) 2 1 En ag For convenience, we define C₁ = Σx=1 7n. we have shown that
E [ S (1 --000) (1-20C,)
Moreover, we find that
7 C
n1 Uk+ni Σ
n-k
Σ df
khty (A) Kik+ni
ktny Vktng
Tr
Let 2 0 and 0 < < 12(2i)), by Chernoff bound (Chernofl 1952 Thm. 1), we obtain
PP ( E | > {(1 + s) e (1+)) 00)
= (1+s))00000 20C;)
We minimize the bound over 0 < < 1)((TT((()) by choosing 0 8/( + s)C,), which gives
IP 12 > {(1 + +)C - (1 8///2
The proof of Theorem 2 will require to bound the term || which we achieve using the following
result (Boullé & Townsend 2022 Lem. 3).
Lemma 8 (Bouléé & Townsend 2022) Let lt f. 2 1 such that k 2 with the notations intro-
duced at the beginning of the section, for all 2 1, we have
P 1 [V 51-(011)
312
k-k
We now prove Theorem 2 which provides a refined probability bound for the performance of the generalized randomized SVD on matrices.
Proof of Theorem 2 Using (Boullé & Townsend, 2022 Thm. 2) and the submultipliaatiity of the
Frobenius norm, we have
AA QQ" E2
Eg
(10)
Let 8 = k + p with p 2 4, combining Lemmas and to bound the terms || EES22113 and ||S2| in
Equation (10 yields the following probability estimate:
E211 + +3/ + 8) k "")
||A - QQ'A <
1
(k+ 3k
< 1 312 + +)
(A))
p+1
7ik
j=k+1
with failure probability at most t-P + (1 (k+ (bbtt)))) Note that we introduced "Kk =
k/(,, T((V,KV, -))) and 3x = Tr Σ V ()/ 12) We conclude the proof by defining /1 + 1.
The following Lemma provides an estimate of the quantity 3g introduced in the statement Theo- rem 2
13

---

Published as a conference paper at ICLR 2022
Lemma Let BR
then the following inequality holds
Bk <
Σ j=
d ))
Σ (A) j=++1
Proof. Let [l ... Ma-k be the eigenvalues of the matrix V2KV2.
Using Von Neumann's
trace inequality (Mirsky 1975
Von Neumann 1937), we have
TV( & E ()) 5=
Then, the matrix V VV is principal submatrix of V*K,, which has the same eigenvalues of K. Therefore, by Kato 2013 Thm. 6.46), the eigenvalues of ViKV2 are individually bounded by the eigenvalues of K, Lee, My I by for S S n k, which concludes the proof.
Finally, we highlight that the statement of Theorem 2 can be simplified by choosing p 5, t 4, and u = 3.
Corollary (Generalized randomized SVD) Let be an x n matrix and 2 an integer. If Ω E Bkx(k+5) is Gaussian random matrix, where each column is sampled from a multivariate Gaussian distribution with symmetric positive semi-definite covariance matrix K E RRaxn, and QR = AS is the economized QR decomposition of AS,, then
P A - QQ"AlE < (1+9 k(k + +5)
BK Th
Σ (A) j=++1
> 0.999.
In contrast, simplification of Theorem by choosing t = and u = gives the following result.
Corollary (Randomized SVD) Let A be an n X matrix and k 2 an integer. If E (k+))) a standard Gaussian random matrix and QR A.t the economized QR decomposition of AΩ, then
P A - QQ' S (01A:=((++66V/))) (1+16Vk+3)
Σ (2 J=++1
0.999.
B CONTINUITY OF THE JACOBI KERNEL
Here, we show that if the kernel's eigenvalue sequence, { is such that ES j < co. then the Jacobi kernel (=(2.2 defined by Equation continuous. First, note that p(2.2) is scaled ultraspherical polynomial 6(5/2) with parameter 5/2 and degree j 2 so can be bounded by the following proposition.
Proposition 10 Let (1(/2) be the ultraspherical polynomial of degree j with parameter 5/2, nor- =2)20 malized such that [_ (1-1122)0(00 (x:² dr = 1. Then,
[11 - 23)0 (z) < max xE[-1,1]
j 0.
(11)
Proof. Let 2 0 and
E [-1, 1). according to NIS Table 18.3.1),
j + 5/2 0/ (x) =3 (j+ 1)(j + 2j + 3) + 4)
C (1)2)
(12)
where c(0/2)(z)) (a:) is the standard ultraspherical polynomial. Using NNS
(18.9.8)), we have
(1 - x (6122) (x)= (33 +)C) (3/2) (x) (j 1)(j 2) (x)
6(++5/2)
14

---

Published as a conference paper at ICLR 2022
By using NNS (18.9.7)), have (0)22((0) (x) 20 (x))/2 (j 5/2) (1/2) (x) and hence,
(+
(1
2 (
(2/2)) (x).
3
j+2
We bound the two terms with
NIS (18.14.4) to obtain the following inequalities:
[( (x) 1i++++++)
1.
2
3+2
Hence, (1- (11- (x)1 s 2( + 1)(j + 2)/3 and following Equation 12)
we obtain
(++ (+10 + +
1(1 (0-=2) (x) /
S 22/3 5/12,
which concludes the proof.
The bound given in Proposition differs by a factor of 4/3 from the numerically observed upper bound (1.5 /j 5 5/12). This result shows that if -ー=1 jAj < co, then the series of functions in
Equation 5 converges uniformly and k/2.2) is continuous.
Proof Lemma By Markov brothers' inequality (Markov 1889 for all j 2 0 and S μ S j. we have max lp) (x)] $ Therefore, 1f% ((()) (x)
Σ
1= a slllp
(4) 1l I Eド=o jak-N
so 1 (x) < do if (v 1)/2. The result follows from (Rudin 1976
Thm. 7.17).
Proof Theorem 4 Since f ~ gP(O, K
x²) p(2.2) (æ), where ej ~ (0, 1) for j 2 0. Let fn denote the truncation after n terms. By letting M be the constant such that A++1 I M( 1)-", we find that
Σ5=u Gj Gi
A5t1(1 CiA2s1111
- Inllo S Sn,
2v M
Σ j=++2
16
where we used max,e(-1.1] (1 - x²) x²) p(2.2) (x)] S 22 +1 (ef. Proposition 10. Thus, we have
P lim |lf fn = =0) 2 IP lim Sn 0):
Here, Sn ~ X, Ej=n+2 Y₁j (2-022 where Y, follows a half-normal distribution (Leone et al.
1961
For with parameter o and the (Yj), are independent. We want to show that X, 4.8. >>> U, using Chebyshev's inequality, have:
P(XX 2 c) <
Σ J=n+2 1
1-121
1
2
R=0
==0
n=1
which is finite if 3. Therefore, using Manee-Cantelli Lemma DDurrett 2019 converges to 0 almost surely and Pliimmn-oe. X, = 0) = Finally,
Chapt. 2.3), X,.
lim 1f - falloc = =0) > lim X2 = =0) =1, n-
P
which proves that } converges uniformly and hence f is continuous with probability one. The statement for higher order derivatives follows the proof of Lemma 3
C EXPERIMENTS WITH NON-STANDARD COVARIANCE MATRICES
In this section, we provide more detailed explanations regarding the numerical experiments con- ducted in Section 5.
Sampling a random vector from a multivariate normal distribution with an arbitrary covariance ma- trix K can be computationally expensive when the dimension, of the matrix is large it requires
15

---

Published as a conference paper at ICLR 2022
the computation of Cholesky factorization, which be done in O operations. We high- light that this step can be precomputed once, such that the overhead of the generalized SVD can be essentially expressed the cost of an extra matrix-vector multiplication. Then, the difference in timings between standard and prior covariance matrices is marginal as shown by the right panel of would like highlight that prior covariance matrices and derived using physical knowledge of the problem, such as diffusion behaviour of the sys- tem, which can also significantly decrease the precomputation
Figure 2 Additionally,
be designed
As an example, in Section 5.1 we employ the discretized Green's s function of the negative Laplacian operator with homogeneous Dirichlet boundary conditions, given by Lu -d dr on [0, 1], for which we know the eigen- value decomposition. Hence, the eigenvalues and normalized eigenfunctions are respectively given by
in (x) = sin(nrz)))
x E
[0, 1], n 1
Therefore, one can employ Mercer's representation (see Equation 4 to sample the random vectors and precompute the covariance matrix in (n²²) operations. For problem of size 2000, takes 0.16s to precompute the matrix.
16


Reviews:
[{"review_id": "o10u7_sv5fU", "reviewer": "Reviewer_uuEP", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "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": "\n0) In abstract: \"Here, we generalize the theory of randomized SVD to to multivariate Gaussian vectors, [...]\"\n\n- A theory of randomized SVD already exists for multivariate Gaussian vectors though? Furthermore, read 1). \n\n1) \"Generalization 1 [...]\"\n\n- This bit is unclear. It seems like the authors are claiming to generalise the randomised SVD for Gaussian vectors with non-identity covariance. However, in the beginning of Sec 3 this was already developed by Boulle & Townsend 2021. I guess what is meant is that new theoretical results pertaining to the approximation error are obtained? I would expect this to be more clear in a  camera-ready version. Also I guess this is the same issue as in abstract (point 0) above). \n\n2) \"In particular, this highlights that a suitably chosen covariance matrix can outperform the randomized SVD with standard Gaussian vectors\"\n\n- Does the converse hold as well? i.e. can an ill-suited covariance matrix make the algorithm worse?\n\n3) \"Therefore, each column of \\Omega is an object, consisting of a polynomial representation [...]\"\n\n- Why polynomials instead of something else? \n\n4) Bad figure ordering / Layout.\n\n-The first Figure to be references is Fig. 4 (rather than 1) on page 4. However, the actual figure does not appear until page 9. \n\n5) Bad notation.\n\n-L is used both to denote the Cholesky factor of a kernel and the space of square integrable functions. E.g. L^2(D) could be interpreted as the square of the Cholesky factor evaluated at D. \n\n6) Bad enumeration notation. In the first paragraph of page 5 items are enuemrated as (1), (2), and (3). However this is the same notation used to reference numbered equations making the paragraph quite hard to read. \n\n7) \"A desireable property of a covariance kernel is to be unbiased towards one spatial direction\"\n\n-This would require further motivation. Would it not depend on application / prior knowledge of the problem? \n\n8) \"[...], we observe a large variation of the randomly generated functions near x = \\pm 1, [...]\"\n\n-Upon examining Figure 1 it seems to me the opposite is true, i.e. the largest variation is found in mid-interval while it is the smallest at the end-points. Please clarify. \n\n9) Figure 2. \n\n- I believe the authors might get their point across better if (additionally) the error is plotted vs time. \n\n10) Proofs. \n\n- The proofs would be easier to follow if the authors actually state what results they are using in their arguments instead of e.g. just referencing (Halko et al, 2011 Prop A.1). \n\n11) Proof of Lemma 7. \n\n-It is not stated what is used how to retrieve the upper bound of the moment generating function. \n\n12) Proof of Theorem 2. \"[...] combining (Boulle & Townsend, 2021, Thm 3.2) and Lemma 7 in Equation (10) [...]\"\n\n- Lemma 7 is not given in Equation (10). Also the proof would be more readable if the result by Boulle and Townsend was stated in the text. \n\n13)  \n\n", "summary_of_the_paper": "This paper studies the randomised SVD algorithm with non-standard Gaussian vectors which is then used to approximate Hilbert-Schmidt operators using a new kernel based on out products of weighted Jacobi polynomials.", "summary_of_the_review": "While the paper appears to adduce some novel results. However, the presentation is unclear which makes it hard to verify the results of the authors and properly situate it in the context of previous literature. \n\nI do however think that the authors have a reasonable chance to address my concerns for a camera-ready version, why I recommend a marginal accept. \n\n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "jOAhojbZcom", "reviewer": "Reviewer_ksBY", "summary": "", "questions": "", "limitations": "", "rating": 8, "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": "Randomized Singular Value Decomposition (SVD) is an important algorithm in the field of numerical linear algebra. Suppose $A \\in \\mathbb{R}^{m \\times n}$ is a large matrix. It is often the case that although $m$ and $n$ are large, the rank of the matrix $A$ is small, and therefore a low rank approximation $A \\approx B C$ for $B \\in \\mathbb{R}^{m \\times k}$ and $A \\in \\mathbb{R}^{k \\times n}$ and $k \\ll \\min(m,n)$ is desired. Randomized SVD computes a rank-$k$ approximation for $A$. The recipe for doing that is to first compute an approximate basis for the range space of $A$. Specifically, we want an orthonormal matrix $Q \\in \\mathbb{R}^{m \\times k}$ such that $A \\approx Q Q^* A$. Randomized SVD does this by letting $\\Omega \\in \\mathbb{R}^{n \\times k}$ be a matrix whose entries are i.i.d. standard Gaussian, and then letting $Q$ be the orthonormal matrix from the QR-decomposition $Q R = A \\Omega$. Once we have $Q$, we have a low rank approximation $A \\approx Q B$ for $B = Q^* A$. The next step is to compute the SVD of the small matrix $B = \\tilde{U} \\Sigma V*$, which finally gives the approximate SVD for $A$ as $U \\Sigma V^*$ for $U = Q \\tilde{U}$.\n\nThis paper starts by first generalizing the above recipe by letting the columns of $\\Omega$ to be sampled from a multivariate Gaussian distribution with a general covariance matrix $K \\in \\mathbb{R}^{n \\times n}$. The error $\\lVert A - Q Q^* A \\rVert_F$ is upper bounded, which is similar to the upper bound in vanilla Randomized SVD. Then this method is also generalized to approximation of Hilbert-Schmidt operators $\\mathscr{F}$, instead of matrices $A$. The samples are now coming from a Gaussian process with a covariance kernel $K$. As is expected, the approximation quality of the randomized algorithms is closely tied to the spectrum of the covariance matrix/kernel $K$. To this end, some intuition is provided for the relation between the choice of the covariance kernel and the approximation ability of the algorithm. Finally, numerical experiments provide justification for these generalizations by showcasing their power.\n\nThe only minor complaint about the paper is that the presentation could be made less terse, especially in Section 3.1 where Randomized SVD for Hilbert-Schmidt operators is discussed. For example, it would be nice to provide some background on how QR algorithm is used on $Y$ to get $Q$, or why each column of $\\Omega$ consists of \"polynomial representation of a smooth random function sampled from the GP in the Chebyshev basis\". It's also unfortunate that no results are provided for the approximation errors in the case of Hilbert-Schmidt operators.\n\nI also didn't understand why the quantity $\\gamma_k$ measures \"the quality of the covariance kernel $K$ in $\\mathcal{GP}(0, K)$ to generate random functions that can learn the H-S operator $\\mathscr{F}$\" claimed in Section 4.2.", "summary_of_the_paper": "This paper proposes and analyzes a generalization of the randomized SVD to incorporate any covariance matrix and to Hilbert-Schmidt operators. Extensive numerical experiments further strengthen the case for this generalization.", "summary_of_the_review": "The paper should be accepted because it is well-written for the most part, tackles an important problem, and provides a novel generalization of a well-known method.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "WPoXthTdjZA", "reviewer": "Reviewer_FGbe", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "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": "Basically, I like the paper's motivation; it is quite a fundamental research problem to reduce computational time by using randomized SVD and improve its approximation quality. Besides, this paper is well structured, and the theoretical background of the proposed approach is well described in the paper. However, I have the following several concerns to the paper:\n\nIn the experiments, it uses only synthetic datasets. Therefore, it is unclear whether the proposed approach is useful for real applications. Typically, if the given matrix is dense, singular values decay rapidly. If the given matrix is sparse, singular values decay slowly. I want to know whether the proposed approach can achieve high approximation accuracy for various types of matrices. \n\nIn the experiment of Figure 2, it seems that the computational time of the Cholesky factorization is not included. I think this is not a fair comparison to the previous approach. Since the computational cost of the Cholesky factorization is cubic according to the size of the covariance matrix, the proposed approach would need a large computational time for the Cholesky factorization. In contrast, the previous approach does not need to compute Cholesky factorization. I want to know the end-to-end computational times of the proposed and previous approach for the specific target rank of SVD. \nIn terms of computational time, the paper should show the theoretical computational cost of the proposed approach. \n\nAs for the approximation quality, the previous approach, as well as the proposed approach, can more accurately perform SVD as the oversampling parameter increases. So, I think the paper should show the experimental results of the approximation quality of each approach against the oversampling parameter. The setting of the oversampling parameter in the experiments is unclear from the descriptions of the paper. ", "summary_of_the_paper": "This paper proposes a generalization approach for the randomized SVD. In the existing approach of the randomized SVD, a standard Gaussian random matrix is used to reduce the size of a matrix to perform SVD. The proposed approach uses a multivariate Gaussian distribution with a covariance matrix instead of the standard Gaussian random matrix in performing SVD. Since the covariance matrix can reflect prior knowledge of the given matrix, the proposed approach can improve the approximation quality of the randomized SVD. By using a synthetic dataset, this paper conducts experiments to show the effectiveness of the proposed approach.", "summary_of_the_review": "This paper is well-written. \nReal datasets should be used in the experiments. \nThe computational time of the Cholesky factorization should be included in Figure 2. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper studies a generalization of the randomized SVD algorithm with non-standard Gaussian vectors, which is then used to incorporate any covariance matrix and to Hilbert-Schmidt operators.  It uses a new kernel related to products of weighted Jacobi polynomials; and extensive numerical experiments further strengthen the case for this generalization.  Reviewers had initial concerns that were addressed, and the method should be of broad interest.
Decision: Accept (Poster)