A GENERALIZATION OF THE RANDOMIZED SINGULAR VALUE DECOMPOSITION
Year: 2022

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 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.

heading: 1 INTRODUCTION
text: 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 et al., 2000; Van Der Veen et al., 1993). The SVD of an m × n real matrix A with m ≥ n is a factorization of the form A = UΣV∗, where U is an m × m orthogonal matrix of left singular vectors, Σ is an m × n diagonal matrix with entries σ1(A) ≥ · · · ≥ σn(A) ≥ 0, and V is an n × n orthogonal matrix of right singular vectors (Golub & Van Loan, 2013). The SVD plays a central role because truncating it after k terms provides a best rank k approximation to A in the spectral and Frobenius norm (Eckart & Young, 1936; Mirsky, 1960). Since computing the SVD of a large matrix can be computationally infeasible, there are various alternative algorithms that compute near-best rank k matrix approximations from matrix-vector products (Halko et al., 2011; Martinsson & Tropp, 2020; Nakatsukasa, 2020; Nyström, 1930; Williams & Seeger, 2001). The randomized SVD uses matrixvector products with standard Gaussian random vectors and is one of the most popular algorithms for constructing a low-rank approximation to A (Halko et al., 2011; Martinsson & Tropp, 2020).
Currently, the randomized SVD is used and theoretically justified when it uses matrix-vector products with standard Gaussian random vectors. In this paper, we consider the following generalizations. 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 a 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, 2015). We design a practical algorithm for learning HS operators using random input functions, sampled from a Gaussian process (GP). Examples of applications include learning integral kernels such as Green’s functions associated with linear partial differential equations (Boullé & Townsend, 2022; Boullé et al., 2021).
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 to introduce a new covariance kernel
based on weighted Jacobi polynomials for learning HS operators. One of the main advantages of this kernel is that it is directly expressed as a Karhunen–Loève expansion (Karhunen, 1946; Loève, 1946) so that it is faster to sample functions from the associated GP than using a standard squaredexponential 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 kernel’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 it is beneficial to use nonstandard Gaussian random vectors in the randomized SVD for constructing low-rank approximations.
2. We generalize the randomized SVD to HS operators and provide numerical examples to learn integral kernels.
3. We propose a 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.

heading: 2 BACKGROUND: THE RANDOMIZED SVD FOR MATRICES
text: The randomized SVD computes a near-best rank k approximation to a matrixA. First, one performs the matrix-vector products y1 = Ax1, . . . , yk+p = Axk+p, where x1, . . . , xk+p are standard Gaussian random vectors with identically and independently distributed entries and p ≥ 1 is an oversampling parameter. Then, one computes the economized QR factorization [y1 · · · yk+p] = QR, before forming the rank ≤ k + p approximant QQ∗A. Note that if A is symmetric, one can form QQ∗A by computing Q(AQ)∗ via matrix-vector products involving A. The quality of the rank ≤ k + p approximant QQ∗A is characterized by the following theorem. Theorem 1 (Halko et al. (2011)) Let A be an m × n matrix, k ≥ 1 an integer, and choose an oversampling parameter p ≥ 4. If Ω ∈ Rn×(k+p) is a standard Gaussian random matrix and QR = AΩ is the economized QR decomposition of AΩ, then for all u, t ≥ 1,
‖A−QQ∗A‖F ≤ ( 1 + t √ 3k
p+ 1 )√√√√ n∑ j=k+1 σ2j (A) + ut √ k + p p+ 1 σk+1(A), (1)
with failure probability at most 2t−p + e−u 2 . The squared tail of the singular values of A, i.e., √∑n
j=k+1 σ 2 j (A), gives the best rank k approxi-
mation error to A in the Frobenius norm. This result shows that the randomized SVD can compute a near-best low-rank approximation to A with high probability.

heading: 3 GENERALIZED RANDOMIZED SVD FOR MATRICES AND OPERATORS
text: The theory behind the randomized SVD has been recently extended to nonstandard covariance matrices and HS operators (Boullé & Townsend, 2022). However, the probability bounds, generalizing Theorem 1, 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 on the expectation is also available in the Appendix (see Lemma 7). While Theorem 2 is formulated with matrices, the same result holds for HS operators in infinite dimensions.
For a fixed target rank 1 ≤ k ≤ n, we define V1 ∈ Rn×k and V2 ∈ Rn×(n−k) to be the matrices containing the first k and last n − k right singular vectors of A, respectively, and denote by Σ2 ∈ R(n−k)×(n−k), the diagonal matrix with entries σk+1(A), . . . , σn(A). We consider a symmetric positive semi-definite covariance matrix K ∈ Rn×n, with kth largest eigenvalue λk > 0. Theorem 2 Let A be an m × n matrix, k ≥ 1 an integer, and choose an oversampling parameter p ≥ 4. If Ω ∈ Rn×(k+p) is a Gaussian random matrix, where each column is sampled from a
multivariate Gaussian distribution with covariance matrix K ∈ Rn×n, and QR = AΩ is the economized QR decomposition of AΩ, then for all u, t ≥ 1,
‖A−QQ∗A‖F ≤ ( 1 + ut √ (k + p) 3k
p+ 1 βk γk )√√√√ n∑ j=k+1 σ2j (A), (2)
with failure probability at most t−p + [ue−(u 2−1)/2]k+p. Here, γk = k/(λ1 Tr((V∗1KV1) −1))) denotes the covariance quality factor, and βk = Tr(Σ22V ∗ 2KV2)/(λ1‖Σ2‖2F), where λ1 is the largest eigenvalue of K.
The factors γk and βk, measuring the quality of the covariance matrix to learn A in Theorem 2, can be respectively bounded (Boullé & Townsend 2022, Lem. 2; Lemma 9) using the eigenvalues λ1 ≥ · · · ≥ λn of the covariance matrix K and the singular values of A as:
1 γk ≤ 1 k n∑ j=n−k+1 λ1 λj , βk ≤ n∑ j=k+1 λj−k λ1 σ2j (A) / n∑ j=k+1 σ2j (A).
This shows that the performance of the generalized randomized SVD depends on the decay rate of the sequence {λj}. The quantities αk and βk depend on how much prior information of the k+ 1, . . . , n right singular vectors of A is encoded in K. In the ideal situation where these singular vectors are known, then one can define K such that βk = 0 for λk+1 = · · · = λn = 0. In particular, this highlights that a suitably chosen covariance matrix can outperform the randomized SVD with standard Gaussian vectors (see Section 5.1 for a numerical example).

heading: 3.1 RANDOMIZED SVD FOR HILBERT–SCHMIDT OPERATORS
text: We now describe the randomized SVD for learning HS operators (see Algorithm 1). The algorithm is implemented in the Chebfun software system (Driscoll et al., 2014), which is a MATLAB package for computing with functions. The Chebfun implementation of the randomized SVD for HS operators uses Chebfun’s capabilities, which offer continuous analogues of several matrix operations like the QR decomposition and numerical integration. Indeed, the continuous analogue of a matrixvector multiplication AΩ for an HS integral operator F (see Hsing & Eubank, 2015 for definitions and properties of HS operators), with kernel G : D ×D → R, is
(Ff)(x) = ∫ D G(x, y)f(y) dy, x ∈ D, f ∈ L2(D),
where D ⊂ Rd with d ≥ 1, and L2(D) is the space of square-integrable functions on D.
Algorithm 1 Randomized SVD for HS operators Input: HS integral operator F with kernel G(x, y), number of samples k > 0 Output: Approximation Gk of G
1: Define a GP covariance kernel K 2: Sample the GP k times to generate a quasimatrix of random functions Ω = [f1 . . . fk] 3: Evaluate the integral operator at Ω, Y = [F (f1) . . .F (fk)] 4: Orthonormalize the columns of Y , Q = orth(Y ) = [q1 . . . qk] 5: Compute an approximation to G by evaluating the adjoint of F 6: Initialize Gk(x, y) to 0 7: for i = 1 : k do 8: Gk(x, y)← Gk(x, y) + qi(x) ∫ D G(z, y)qi(z) dz 9: end for
The algorithm takes as input an integral operator that we aim to approximate. Note that we focus here on learning an integral operator, but other HS operators would work similarly. The first step of the randomized SVD for HS operators consists of generating an∞× k quasimatrix Ω by sampling a GP k times, where k is the target rank (see Section 4). An ∞ × k quasimatrix is an ordered set of k functions (Townsend & Trefethen, 2015), and generalizes the notion of matrix to infinite dimensions. Therefore, each column of Ω is an object, consisting of a polynomial approximation of
a 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 an 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 F or, equivalently, the evaluation of its adjoint Ft satisfying:
(Ftf)(x) = ∫ D G(y, x)f(y) dy, x ∈ D.
We evaluate the adjoint of F at each column vector of Q to construct an approximation Gk of G. Finally, the approximation error between the operator kernel G and the learned kernel Gk can be computed in the L2-norm, corresponding to the HS norm of the integral operator.

heading: 4 COVARIANCE KERNELS
text: To generate the random input functions f1, . . . , fk for the randomized SVD for HS operators, we draw them from a GP, denoted by GP(0,K), for a certain covariance kernel K. A widely employed covariance kernel is the squared-exponential function KSE (Rasmussen & Williams, 2006) given by
KSE(x, y) = exp ( −|x− y|2/(2`2) ) , x, y ∈ D, (3)
where ` > 0 is a parameter controlling the length-scale of the GP. This kernel is isotropic as it only depends on |x − y|, is infinitely differentiable, and its eigenvalues decay supergeometrically to 0. Since the bound in Theorem 2 degrades as the ratio λ1/λj increases for j ≥ 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 a rank greater than k̃, where k̃ is such that the k̃th eigenvalue of KSE reaches machine precision.
Other popular kernels for GPs include the Maté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-exponential periodic kernel:
KPer(x, y) = exp ( − 2 `2 sin2 ( x− y 2 )) , x, y ∈ D,
where ` > 0 is the length-scale parameter.

heading: 4.1 SAMPLE RANDOM FUNCTIONS FROM A GAUSSIAN PROCESS
text: In finite dimensions, a random vector u ∼ N (0,K), where K ∈ Rn×n is a covariance matrix with Cholesky factorization K = LL∗, can be generated from the matrix-vector product u = Lc. Here, c ∈ Rn is a 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 a continuous symmetric positive-definite covariance function defined on the domain [a, b] × [a, b] ⊂ R2 with −∞ < a < b <∞. We consider the continuous analogue of the Cholesky factorization to write K as (Townsend & Trefethen, 2015)
K(x, y) = ∞∑ j=1 rj(x)rj(y) = Lc(x)L ∗ c(y), x, y ∈ [a, b],
where rj is the jth row of Lc, which —in Chebfun’s terminology— is a lower-triangular quasimatrix. In practice, we truncate the series after n terms, either arbitrarily or when the nth largest kernel eigenvalue, λn, falls below machine precision. Then, if c ∈ Rn follows the standard Gaussian distribution, a function u can be sampled from GP(0,K) as u = Lcc. That is,
u(x) = n∑ j=1 cjrj(x), x ∈ [a, b].
The continuous Cholesky factorization is implemented in Chebfun2 (Townsend & Trefethen, 2013), which is the extension of Chebfun for computing with two-dimensional functions. As an example,
the polynomial approximation, which is accurate up to essentially machine precision, of the squaredexponential covariance kernel KSE with parameter ` = 0.01 on [−1, 1]2 yields a numerical rank of n = 503. The functions sampled from GP(0,KSE) become more oscillatory as the length-scale parameter ` decreases and hence the numerical rank of the kernel increases or, equivalently, the associated eigenvalues sequence {λj} decays slower to zero.

heading: 4.2 INFLUENCE OF THE KERNEL’S EIGENVALUES AND MERCER’S REPRESENTATION
text: The covariance kernel can also be defined from its Mercer’s representation as
K(x, y) = ∞∑ j=1 λjψj(x)ψj(y), x, y ∈ D, (4)
where {ψj} is an orthonormal basis of L2(D) and λ1 ≥ λ2 ≥ · · · > 0 (Hsing & Eubank, 2015, Thm. 4.6.5). We prefer to construct K directly from Mercer’s representations for several reasons: 1. One can impose prior knowledge of the kernel of the HS operator on the eigenfunctions of K (such as periodicity or smoothness), 2. One can often generate samples from GP(0,K) efficiently using Equation (4), and 3. One can control the decay rate of the eigenvalues of K,
Hence, the quantity γk in the probability bound of Theorem 2 measures the quality of the covariance kernel K in GP(0,K) to generate random functions that can learn the HS operator F . To minimize 1/γk we would like to select the eigenvalues λ1 ≥ λ2 ≥ · · · > 0 of K so that they have the slowest possible decay rate while maintaining ∑∞ j=1 λj <∞. One needs {λj} ∈ `1 to guarantee that
ω ∼ GP(0,K) has finite expected squared L2-norm, i.e., E[‖ω‖2L2(D)] = ∑∞ j=1 λj <∞. The best sequence of eigenvalues we know that satisfies this property is called the Rissanen sequence (Rissanen, 1983) and is given by λj = Rj := 2−L(j), where
L(j) = log2(c0) + log ∗ 2(j), log ∗ 2(j) = ∞∑ i=2 max(log (i) 2 (j), 0), c0 = ∞∑ i=2 2− log ∗ 2(i),
and log(i)2 (j) = log2 ◦ · · · ◦ log2(j) is the composition of log2(·) i times. Other options for the choice of eigenvalues include any sequence of the form λj = j−ν for ν > 1.

heading: 4.3 JACOBI COVARIANCE KERNEL
text: If D = [−1, 1], then a natural choice of orthonormal basis of L2(D) to define the Mercer’s representation of the kernel are weighted Jacobi polynomials (Deheuvels & Martynov, 2008; Olver et al., 2010). That is, for a weight function wα,β(x) = (1−x)α(1+x)β with α, β > −1, and any positive eigenvalue sequence {λj}, we consider the Jacobi kernel
K (α,β) Jac (x, y) = ∞∑ j=0 λj+1w 1/2 α,β(x)P̃ (α,β) j (x)w 1/2 α,β(y)P̃ (α,β) j (y), x, y ∈ [−1, 1], (5)
where P̃ (α,β)j is the scaled Jacobi polynomial of degree j. The polynomials are normalized such that ‖w1/2α,βP̃ (α,β) j ‖L2([−1,1]) = 1 and K (α,β) Jac ∈ L2([−1, 1]2). In this case, a random function can be sampled as
u(x) = ∞∑ j=0 √ λj+1cjw 1/2 α,βP̃ (α,β) j (x), x ∈ [−1, 1],
where cj ∼ N (0, 1) for 0 ≤ j ≤ ∞. A desirable property of a covariance kernel is to be unbiased towards one spatial direction, i.e., K(x, y) = K(−y,−x) for x, y ∈ [−1, 1], which motivates us to always select α = β. Moreover, it is desirable to have the eigenfunctions of K(α,β)Jac to be polynomial so that one can generate samples from GP(0,K) efficiently. This leads us to choose α and β to be even integers. The choice of α = β = 0 gives the Legendre kernel (Foster et al., 2020; Habermann, 2021). In Section 5, we use Equation (5) with α = β = 2 (see Figure 1) to ensure that functions sampled from the associated GP satisfy homogeneous Dirichlet boundary conditions. In this case, one must have ∑∞ j=1 jλj <∞ so that the series of functions in Equation (5) converges uniformly and K(2,2)Jac is a continuous kernel
(see Appendix B). Under this additional constraint, the best choice of eigenvalues is given by a scaled Rissanen sequence: λj = Rj/j, for j ≥ 1 (cf. Section 4.2). Covariance kernels on higher dimensional domains of the form D = [−1, 1]d, for d ≥ 2, can be defined using tensor products of weighted Jacobi polynomials.

heading: 4.4 SMOOTHNESS OF FUNCTIONS SAMPLED FROM A GP WITH JACOBI KERNEL
text: We now connect the decay rate of the eigenvalues of the Jacobi covariance kernel K(2,2)Jac to the smoothness of the samples from GP(0,K(2,2)Jac ). Hence, the Jacobi covariance function allows the control of the decay rate of the eigenvalues {λj} 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 3 Let {pj} be a family of polynomials such that deg(pj) ≤ j and maxx∈[−1,1] |pj(x)| = 1. If fn(x) = ∑n j=0 ajpj(x) with |aj | ≤ j−ν for ν > 1, then fn converges uniformly to f(x) =∑∞
j=0 ajpj(x) and f is µ times continuously differentiable for any integer µ such that µ < (ν−1)/2.
Note that the main application of this lemma occurs when deg(pj) = j for all j ≥ 0. The proof of Lemma 3 is deferred to the supplementary material. We now state the following theorem about the regularity of functions sampled from GP(0,K(2,2)Jac ), which guarantees that if the eigenvalues are chosen such that λj = O(1/jν) with ν > 3, then f ∼ GP(0,K(2,2)Jac ) is almost surely continuous. Moreover, a faster decay of the eigenvalues of K(2,2)Jac implies higher regularity of the sampled functions, in an almost sure sense. Theorem 4 Let {λj} ∈ `1(R+) be a positive sequence such that λj = O(j−ν) for ν > 3. If f is sampled from GP(0,K(2,2)Jac ), then f ∈ Cµ([−1, 1]) almost surely for any integer µ < (ν − 3)/2.
This theorem can be seen as a particular case of Driscoll’s zero-one law (Driscoll, 1973), which characterizes the regularity of functions samples from GPs (see also Kanagawa et al., 2018).
In Figure 1, we display the Jacobi kernel of type (2, 2) with functions sampled from the corresponding GP. We selected eigenvalue sequences of different decay rates: from the faster 1/j4 to the slower Rissanen sequence Rj/j (Section 4.2). For λj = 1/j3 and λj = Rj/j, we observe a large variation of the randomly generated functions near x = ±1, indicating a potential discontinuity of the samples at these two points as n → ∞. This is in agreement with Theorem 4, which only guarantees continuity (with probability one) of the randomly generated functions if λj ∼ 1/jν with ν > 3.

heading: 5 NUMERICAL EXPERIMENTS

heading: 5.1 APPROXIMATION OF MATRICES USING NON-STANDARD COVARIANCE FUNCTIONS
text: The approximation error bound in Theorem 2 depends on 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 as its covariance matrix, can be improved by using a different covariance matrix. We then attempt to learn the discretized 2000×2000 matrix, i.e., the discrete Green’s function, of the inverse of the following differential operator:
Lu = d2u/dx2 − 100 sin(5πx)u, x ∈ [0, 1]. We vary the number of columns (i.e. samples from the GP) in the input matrix Ω from 1 to 2000.
In Figure 2(left), we compare the ratios between the relative error in 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× 2000 matrix of the Green’s function of the negative Laplace operator Lu = −d2u/dx2 on [0, 1] to incorporate knowledge of the diffusion term in the matrix A. We see that a nonstandard covariance matrix leads to a higher approximation accuracy, with a reduction of the error by a factor of 1.3-1.6 compared to the standard randomized SVD. At the same time, the procedure is only 20% slower1 on average (Figure 2 (right)) as one can precompute the eigenvalue decomposition of the covariance matrix (see Appendix C). It is of interest to maximize the accuracy of the approximation matrix from a limited number of samples in applications where the sampling time is much higher than the numerical linear algebra costs.

heading: 5.2 APPLICATIONS OF THE RANDOMIZED SVD FOR HILBERT–SCHMIDT OPERATORS
text: We now apply the randomized SVD for HS operators to learn kernels of integral operators. In this first example, the kernel is defined as (Townsend, 2013)
G(x, y) = cos(10(x2 + y)) sin(10(x+ y2)), x, y ∈ [−1, 1], and is displayed in Figure 3(a). We employ the squared-exponential covariance kernel KSE with parameter ` = 0.01 and k = 100 samples (see Equation (3)) to sample random functions from the associated GP. The learned kernel Gk is represented on the bottom panel of Figure 3(a) and has an approximation error around machine precision.
1Timings were performed on an Intel Xeon CPU E5-2667 v2 @ 3.30GHz using MATLAB R2020b without explicit parallelization.
As a second application of the randomized SVD for HS operators, we learn the kernel G(x, y) = Ai(−13(x2y + y2)) for x, y ∈ [−1, 1], where Ai is the Airy function (NIS, Chapt. 9) defined by
Ai(x) = 1
π ∫ ∞ 0 cos ( t3 3 + xt ) dt, x ∈ R.
We plot the kernel and its low-rank approximant given by the randomized SVD for HS operators in Figure 3(b) and obtain an approximation error (measured in the L2-norm) of 5.04×10−14. The two kernels have a numerical rank equal to 42.
The last example consists of learning the HS operator associated with the kernel G(x, y) = J0(100(xy+y
2)) for x, y ∈ [−1, 1], where J0 is the Bessel function of the first kind (NIS, Chapt. 10) defined as
J0(x) = 1
π ∫ π 0 cos(x sin t) dt, x ∈ R,
and plotted in Figure 3(c). The rank of this kernel is equal to 91 while its approximation is of rank 89 and the approximation error is equal to 4.88 × 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 4. Here, we vary the number of samples from k = 1 to k = 100 and use the randomized SVD for HS operators with four different covariance kernels: the squared-exponential KSE with parameters ` = 0.01, 0.1, 1, and the Jacobi kernel K(2,2)Jac with eigenvalues λj = 1/j
3, for j ≥ 1. In the left panel of Figure 4, we represent the eigenvalue ratio λj/λ1 of the four kernels and observe that this quantity falls below machine precision for the squared-exponential kernel with ` = 1 and ` = 0.1 at j = 13 and j = 59, respectively. In Figure 4 (right), we observe that these two kernels fail to approximate kernels of high numerical rank. The other two kernels have a much slower decay of eigenvalues and can capture (or learn) more complicated kernels. We then see in the right panel of Figure 4 that the relative approximation errors obtained using K(2,2)Jac and KSE are close to the best approximation error given by the squared tail of the singular values of the integral kernel G(x, y), i.e., ( ∑ j≥k+1 σ 2 j )
1/2. The overshoot in the error at k = 100 compared to the machine precision is due to the decay of the eigenvalues 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
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 K(2,2)Jac and KSE with ` = 0.01 are Error(K(2,2)Jac ) ≈ 2.6× 10−11, and Error(KSE) ≈ 5.7× 10−13, which gives a ratio of
Error(K(2,2)Jac )/Error(KSE) ≈ 45.6. (6) However, the square-root of the ratio of the quality of the two kernels for k = 91 is equal to√
γ91(KSE)/γ91(K (2,2) Jac ) ≈ 117.8, (7)
which is of the same order of magnitude of Equation (6) as predicted by Theorem 2. In Equation (7), γ91(KSE) ≈ 5.88× 10−2 and γ91(K(2,2)Jac ) ≈ 4.24× 10−6 are 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, a 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 a kernel with a large numerical rank. The optimal choice from a theoretical viewpoint is to select a covariance kernel whose eigenvalues have a decay rate similar to the Rissanen sequence (Rissanen, 1983), but other choices may be preferable in practice to ensure smoothness of the sample functions (cf. Section 4.4).

heading: 6 CONCLUSIONS
text: We have explored practical extensions of the randomized SVD to nonstandard covariance matrices and Hilbert–Schmidt operators. This paper motivates new computational and algorithmic approaches for preconditioning the covariance kernel based on prior information to compute a lowrank 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 approximation errors. In addition, we proposed a covariance kernel based on weighted Jacobi polynomials, which allows the control of the smoothness of the random functions generated and may find practical applications in PDE learning (Boullé et al., 2020; 2021) as 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 to 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 to realistic application settings.

heading: ACKNOWLEDGMENTS
text: We thank Joel Tropp for his suggestions leading to Lemma 7 and Daniel Kressner for his comments. This work was supported by the EPSRC Centre For Doctoral Training in Industrially Focused Mathematical Modelling (EP/L015803/1) in collaboration with Simula Research Laboratory, and the National Science Foundation grants DMS-1818757, DMS-1952757, and DMS-2045646.

heading: A RANDOMIZED SVD WITH MULTIVARIATE GAUSSIAN INPUTS
text: Let m ≥ n ≥ 1 and A be an m × n real matrix with singular value decomposition A = UΣV∗, where U and V are orthonormal matrices, and Σ be anm×n diagonal matrix with entries σ1(A) ≥ · · · ≥ σn(A) > 0. For a fixed target rank k ≥ 1, we define Σ1 and Σ2 to be the diagonal matrices, which respectively contain the first k singular values of A: σ1(A) ≥ · · · ≥ σk(A), and the remaining singular values. Let V1 be the n × k matrix obtained by truncating V1 after k columns and V2 the remainder. In this section, K denotes a symmetric positive semi-definite n× n matrix and Ω ∈ Rn×` a Gaussian random matrix with ` ≥ k independent columns sampled from a multivariate normal distribution with covariance matrix K. Finally, we define Ω1 := V∗1Ω and Ω2 := V ∗ 2Ω. We first refine (Boullé & Townsend, 2022, Lem. 5). Lemma 5 Let ` ≥ 1, Ω ∈ Rn×` be a Gaussian random matrix, where each column is sampled from a multivariate Gaussian distribution with covariance matrix K, and T be an `× k matrix. Then,
E[‖Σ2V∗2ΩT‖2F] = Tr(Σ22V∗2KV2)‖T‖2F. (8)
Proof. Let K = QKΛQ∗K be the eigenvalue decomposition of K, where QK is orthonormal and Λ is a diagonal matrix containing the eigenvalues of K in decreasing order: λ1 ≥ · · · ≥ λn ≥ 0. We note that Ω can be expressed as Ω = QKΛ1/2G, where G is a standard Gaussian matrix. Let S = Σ2V∗2QKΛ
1/2, the proof follows from (Halko et al., 2011, Prop. A.1), which shows that E‖SGT‖2F = ‖S‖2F‖T‖2F. Note that one can bound the term Tr(Σ22V ∗ 2KV2) by λ1‖Σ2‖2F, where λ1 is the largest eigenvalue of K (Boullé & 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 A can outperform the randomized SVD with standard Gaussian inputs. Then, combining the proof of (Boullé & Townsend, 2022, Thm. 1) 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 A be an m×n matrix, k ≥ 1 an integer, and choose an oversampling parameter p ≥ 2. If Ω ∈ Rn×(k+p) is a Gaussian random matrix, where each column is sampled from a multivariate Gaussian distribution with covariance matrix K ∈ Rn×n, and QR = AΩ is the economized QR decomposition of AΩ, then,
E [‖A−QQ∗A‖F] ≤ ( 1 + √ βk γk k(k + p) p− 1 )√√√√ n∑ j=k+1 σ2j (A),
where γk = k/(λ1 Tr((V∗1KV1) −1))) and βk = Tr(Σ22V ∗ 2KV2)/(λ1‖Σ2‖2F).
We remark that for standard Gaussian inputs, we have γk = βk = 1 in Proposition 6, and we recover the average Frobenius error of the randomized SVD (Halko et al., 2011, Thm. 10.5) up to a factor of (k + p) due to the non-independence of Ω1 and Ω2 in general. Lemma 7 With the notations introduced at the beginning of the section, for all s ≥ 0, we have
P { ‖Σ2Ω2‖2F > `(1 + s) Tr(Σ22V∗2KV2) } ≤ (1 + s)`/2e−s`/2.
Proof. Let ωj be the jth column of Ω for 1 ≤ j ≤ ` and v1, . . . , vn be the n columns of the orthonormal matrix V. We first remark that
‖Ω2‖2F = ∑̀ j=1 Zj , Zj := n−k∑ n1=1 σ2k+n1(A)(v ∗ k+n1ωj) 2, (9)
where the Zj are i.i.d because ωj ∼ N (0,K) are i.i.d. Let λ1 ≥ · · · ≥ λn ≥ 0 be the eigenvalues of K with eigenvectors ψ1, . . . , ψn ∈ Rn. For 1 ≤ j ≤ `, we have,
ωj = n∑ i=1 (c (j) i ) 2 √ λiψi,
where c(j)i ∼ N (0, 1) are i.i.d for 1 ≤ i ≤ n and 1 ≤ j ≤ `. Then,
Zj = n∑ i=1 (c (j) i ) 2λi n−k∑ n1=1 σ2k+n1(A)(v ∗ k+n1ψi) 2 = n∑ i=1 Xi
where the Xi are independent. Let γi = λi ∑n−k n1=1 σ2k+n1(A)(v ∗ k+n1 ψi) 2, then Xi ∼ γiχ2 for 1 ≤ i ≤ n. Let 0 < θ < 1/(2 ∑n i=1 γi), we can bound the moment generating function of ∑n i=1Xi as
E [ eθ ∑n i=1Xi ] = n∏ i=1 E [ eθXi ] = n∏ i=1 (1− 2θγi)−1/2 ≤ ( 1− 2θ n∑ i=1 γi )−1/2 because the Xi/γi are independent and follow a chi-squared distribution. The right inequality is obtained by showing by recurrence that, if a1, . . . , an ≥ 0 are such that ∑n i=1 ai ≤ 1, then ∏n i=1(1−
ai) ≥ 1− ∑n i=1 ai. For convenience, we define C1 := ∑n i=1 γn, we have shown that
E [ eθZj ] ≤ (1− 2θC1)−1/2.
Moreover, we find that
C1 = n−k∑ n1=1 σ2k+n1v ∗ k+n1 ( n∑ i=1 ψ∗i λiψi ) vk+n1 = n−k∑ n1=1 σ2k+n1(A)v ∗ k+n1Kvk+n1
= Tr(Σ22V ∗ 2KV2).
Let s ≥ 0 and 0 < θ < 1/(2C1), by Chernoff bound (Chernoff, 1952, Thm. 1), we obtain P { ‖Σ2Ω2‖2F > `(1 + s)C1 } ≤ e−(1+s)C1`θE [ eθZj ]` = e−(1+s)C1`θ(1− 2θC1)−`/2.
We minimize the bound over 0 < θ < 1/(2 Tr(K)) by choosing θ = s/(2(1 + s)C1), which gives P { ‖Σ2Ω2‖2F > `(1 + s)C1 } ≤ (1 + s)`/2e−`s/2.
The proof of Theorem 2 will require to bound the term ‖Ω†1‖2F, which we achieve using the following result (Boullé & Townsend, 2022, Lem. 3). Lemma 8 (Boullé & Townsend 2022) Let k, ` ≥ 1 such that ` − k ≥ 4, with the notations introduced at the beginning of the section, for all t ≥ 1, we have
P { ‖Ω†1‖2F > 3t2 Tr((V∗1KV1) −1)
`− k + 1
} ≤ t−(`−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 (Boullé & Townsend, 2022, Thm. 2) and the submultiplicativity of the Frobenius norm, we have
‖A−QQ∗A‖2F ≤ ‖Σ2‖2F + ‖Σ2Ω2‖2F‖Ω † 1‖2F. (10)
Let ` = k + p with p ≥ 4, combining Lemmas 7 and 8 to bound the terms ‖Σ2Ω2‖2F and ‖Ω † 1‖2F in Equation (10) yields the following probability estimate:
‖A−QQ∗A‖2F ≤ ‖Σ2‖2F + 3t2(1 + s) k + p p+ 1 Tr((V∗1KV1) −1) Tr(Σ22V ∗ 2KV2)
≤ ( 1 + 3t2(1 + s) (k + p)k
p+ 1 βk γk ) n∑ j=k+1 σ2j (A),
with failure probability at most t−p + (1 + s)(k+p)/2e−s(k+p)/2. Note that we introduced γk := k/(λ1 Tr((V ∗ 1KV1) −1))) and βk := Tr(Σ22V ∗ 2KV2)/(λ1‖Σ2‖2F). We conclude the proof by
defining u = √
1 + s ≥ 1. The following Lemma provides an estimate of the quantity βk introduced in the statement of Theorem 2.
Lemma 9 Let βk = Tr(Σ22V∗2KV2)/(λ1‖Σ2‖2F), then the following inequality holds
βk ≤ n∑
j=k+1
λj−k λ1 σ2j (A) / n∑ j=k+1 σ2j (A).
Proof. Let µ1 ≥ · · · ≥ µn−k be the eigenvalues of the matrix V∗2KV2. Using Von Neumann’s trace inequality (Mirsky, 1975; Von Neumann, 1937), we have
Tr(Σ22V ∗ 2KV2) ≤ n∑ j=k+1 µj−kσ 2 j (A).
Then, the matrix V∗2KV2 is a principal submatrix of V ∗KV, which has the same eigenvalues of K. Therefore, by (Kato, 2013, Thm. 6.46), the eigenvalues of V∗2KV2 are individually bounded by the eigenvalues of K, i.e., µj ≤ λj for 1 ≤ j ≤ 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 1 (Generalized randomized SVD) Let A be an n × n matrix and k ≥ 1 an integer. If Ω ∈ Rn×(k+5) is a Gaussian random matrix, where each column is sampled from a multivariate Gaussian distribution with symmetric positive semi-definite covariance matrix K ∈ Rn×n, and QR = AΩ is the economized QR decomposition of AΩ, then
P ‖A−QQ∗A‖F ≤ (1 + 9 √ k(k + 5)
βk γk )√√√√ n∑ j=k+1 σ2j (A)  ≥ 0.999. In contrast, a simplification of Theorem 1 by choosing t = 6 and u = 4 gives the following result.
Corollary 2 (Randomized SVD) Let A be an n×n matrix and k ≥ 1 an integer. If Ω ∈ Rn×(k+5) is a standard Gaussian random matrix and QR = AΩ is the economized QR decomposition of AΩ, then
P ‖A−QQ∗A‖F ≤ (1 + 16√k + 5) √√√√ n∑ j=k+1 σ2j (A)  ≥ 0.999.

heading: B CONTINUITY OF THE JACOBI KERNEL
text: Here, we show that if the kernel’s eigenvalue sequence, {λj}, is such that ∑∞ j=1 jλj < ∞, then the Jacobi kernel K(2,2)Jac defined by Equation (5) is continuous. First, note that P̃ (2,2) j is a scaled ultraspherical polynomial C̃(5/2)j with parameter 5/2 and degree j ≥ 0 so it can be bounded by the following proposition.
Proposition 10 Let C̃(5/2)j be the ultraspherical polynomial of degree j with parameter 5/2, normalized such that ∫ 1 −1(1− x 2)2C̃ (5/2) j (x) 2 dx = 1. Then,
max x∈[−1,1]
|(1− x2)C̃(5/2)j (x)| ≤ 2 √ j + 5/12, j ≥ 0. (11)
Proof. Let j ≥ 0 and x ∈ [−1, 1], according to (NIS, Table 18.3.1),
C̃ (5/2) j (x) = 3
√ j + 5/2
(j + 1)(j + 2)(j + 3)(j + 4) C
(5/2) j (x), (12)
where C(5/2)j (x) is the standard ultraspherical polynomial. Using (NIS, (18.9.8)), we have
(1− x2)C(5/2)j (x) = (j + 3)(j + 4)C
(3/2) j (x)− (j + 1)(j + 2)C (3/2) j+2 (x)
6(j + 5/2) .
By using (NIS, (18.9.7)), we have (C(3/2)j+2 (x)− C (3/2) j (x))/2 = (j + 5/2)C (1/2) j+2 (x) and hence,
(1− x2)C(5/2)j (x) = 2
3 C
(3/2) j (x)−
(j + 1)(j + 2)
3 C
(1/2) j+2 (x).
We bound the two terms with (NIS, (18.14.4)) to obtain the following inequalities:
|C(3/2)j (x)| ≤ (j + 1)(j + 2)
2 , |C(1/2)j+2 (x)| ≤ 1.
Hence, |(1− x2)C(5/2)j (x)| ≤ 2(j + 1)(j + 2)/3 and following Equation (12) we obtain
|(1− x2)C̃(5/2)j (x)| ≤ 2
√ (j + 1)(j + 2)(j + 5/2)
(j + 3)(j + 4) ≤ 2 √ j + 5/12,
which concludes the proof.
The bound given in Proposition 10 differs by a factor of 4/3 from the numerically observed upper bound (1.5 √ j + 5/12 ). This result shows that if ∑∞ j=1 jλj < ∞, then the series of functions in Equation (5) converges uniformly and K(2,2)Jac is continuous.
Proof of Lemma 3. By Markov brothers’ inequality (Markov, 1889), for all j ≥ 0 and 0 ≤ µ ≤ j, we have maxx∈[−1,1] |p (µ) j (x)| ≤ j2µ. Therefore, |f (µ) n (x)| ≤ ∑n j=0 |aj |‖p (µ) j ‖∞ ≤ ∑n j=0 j 2µ−ν so |f (µ)n (x)| <∞ if µ < (ν − 1)/2. The result follows from (Rudin, 1976, Thm. 7.17).
Proof of Theorem 4. Since f ∼ GP(0,K(2,2)Jac ), f ∼ ∑∞ j=0 cj √ λj+1(1 − x2)P̃ (2,2)j (x), where cj ∼ N (0, 1) for j ≥ 0. Let fn denote the truncation of f after n terms. By letting M > 0 be the constant such that λj+1 ≤M(j + 1)−ν , we find that
‖f − fn‖∞ ≤ Sn, Sn := 2 √ M ∞∑ j=n+2 |cj−1|j(1−ν)/2,
where we used maxx∈[−1,1] |(1− x2)P̃ (2,2) j (x)| ≤ 2
√ j + 1 (cf. Proposition 10). Thus, we have
P (
lim n→∞
‖f − fn‖∞ = 0 ) ≥ P ( lim n→∞ Sn = 0 ) .
Here, Sn ∼ Xn = ∑∞ j=n+2 Yjj
(1−ν)/2, where Yj follows a half-normal distribution (Leone et al., 1961) with parameter σ = 1 and the (Yj)j are independent. We want to show that Xn
a.s.−−→ 0. For > 0, using Chebyshev’s inequality, we have:
∞∑ n=0 P(|Xn| ≥ ) ≤ 1 2 ∞∑ n=0 ( 1− 2 π ) ∞∑ j=n+2 1 jν−1 ≤ 1 2 ( 1− 2 π ) 1 ν − 2 ∞∑ n=1 1 nν−2 ,
which is finite if ν > 3. Therefore, using Borel–Cantelli Lemma (Durrett, 2019, Chapt. 2.3), Xn converges to 0 almost surely and P(limn→∞Xn = 0) = 1. Finally,
P (
lim n→∞
‖f − fn‖∞ = 0 ) ≥ P ( lim n→∞ Xn = 0 ) = 1,
which proves that {fn} converges uniformly and hence f is continuous with probability one. The statement for higher order derivatives follows the proof of Lemma 3.

heading: C EXPERIMENTS WITH NON-STANDARD COVARIANCE MATRICES
text: In this section, we provide more detailed explanations regarding the numerical experiments conducted in Section 5.1.
Sampling a random vector from a multivariate normal distribution with an arbitrary covariance matrix K can be computationally expensive when the dimension, n, of the matrix is large as it requires
the computation of a Cholesky factorization, which can be done in O(n3) operations. We highlight that this step can be precomputed once, such that the overhead of the generalized SVD can be essentially expressed as 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 Figure 2. Additionally, we would like to highlight that prior covariance matrices can be designed and derived using physical knowledge of the problem, such as a diffusion behaviour of the system, which can also significantly decrease the precomputation cost. As an example, in Section 5.1, we employ the discretized Green’s function of the negative Laplacian operator with homogeneous Dirichlet boundary conditions, given by Lu = −d2u/dx2 on [0, 1], for which we know the eigenvalue decomposition. Hence, the eigenvalues and normalized eigenfunctions are respectively given by
λn = 1
π2n2 , ψn(x) =
√ 2 sin(nπx), x ∈ [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 O(n2) operations. For a problem of size n = 2000, it takes 0.16s to precompute the matrix.

Review summaries:
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.