RANGE-NET: A HIGH PRECISION NEURAL SVD
Year: 2021

Abstract:
For Big Data applications, computing a rank-r Singular Value Decomposition (SVD) is restrictive due to the main memory requirements. Recently introduced streaming Randomized SVD schemes work under the restrictive assumption that the singular value spectrum of the data has an exponential decay. This is seldom true for any practical data. Further, the approximation errors in the singular vectors and values are high due to the randomized projection. We present Range-Net as a low memory alternative to rank-r SVD that satisfies the lower bound on tail-energy given by Eckart-Young-Mirsky (EYM) theorem at machine precision. Range-Net is a deterministic two-stage neural optimization approach with random initialization, where the memory requirement depends explicitly on the feature dimension and desired rank, independent of the sample dimension. The data samples are read in a streaming manner with the network minimization problem converging to the desired rank-r approximation. Range-Net is fully interpretable where all the network outputs and weights have a specific meaning. We provide theoretical guarantees that Range-Net extracted SVD factors satisfy EYM tail-energy lower bound with numerical experiments on real datasets at various scales that confirm these bounds. A comparison against the state-of-the-art streaming Randomized SVD shows that Range-Net is six orders of magnitude more accurate in terms of tail energy while correctly extracting the singular values and vectors.

text: For Big Data applications, computing a rank-r Singular Value Decomposition (SVD) is restrictive due to the main memory requirements. Recently introduced streaming Randomized SVD schemes work under the restrictive assumption that the singular value spectrum of the data has an exponential decay. This is seldom true for any practical data. Further, the approximation errors in the singular vectors and values are high due to the randomized projection. We present Range-Net as a low memory alternative to rank-r SVD that satisfies the lower bound on tail-energy given by Eckart-Young-Mirsky (EYM) theorem at machine precision. Range-Net is a deterministic two-stage neural optimization approach with random initialization, where the memory requirement depends explicitly on the feature dimension and desired rank, independent of the sample dimension. The data samples are read in a streaming manner with the network minimization problem converging to the desired rank-r approximation. Range-Net is fully interpretable where all the network outputs and weights have a specific meaning. We provide theoretical guarantees that Range-Net extracted SVD factors satisfy EYM tail-energy lower bound with numerical experiments on real datasets at various scales that confirm these bounds. A comparison against the state-of-the-art streaming Randomized SVD shows that Range-Net is six orders of magnitude more accurate in terms of tail energy while correctly extracting the singular values and vectors.

heading: 1 INTRODUCTION
text: Singular Value Decomposition (SVD) is pivotal to exploratory data analysis in identifying an invariant structure under a minimalistic representation (assumptions on the structure) containing the span of resolvable information in the dataset. Finding a low rank structure is a fundamental task in applications including Image Compression (de Souza et al., 2015), Image Recovery (Brand, 2002), Background Removal (Wang et al., 2018), Recommendation Systems (Zhang et al., 2005) and as a pre-processing step for Clustering (Drineas et al., 2004) and Classification (Jing et al., 2017). With the advent of digital sensors and modern day data acquisition technologies, the sheer amount of data now requires that we revisit the solution scheme with reduced memory consumption as the target. In this work, we reformulate SVD with special emphasis on the main memory requirement, with no loss in accuracy, that precludes it’s use for big data applications.
It is well known that natural data matrices have a decaying spectrum wherein saving the data in memory in its original form is either redundant or not required from an application point of view. However, any assumption on the decay rate can only be validated if the singular value decomposition is known a priori, which is seldom the case in exploratory data analysis (see Fig. 3). Visually assessing a rank-r approximation for image processing applications might seem correct qualitatively (see Fig. 4) but are still prone to large errors due to limited human vision acuity (see Fig. 6). This is further exacerbated when the application at hand is associated with scientific computations wherein the anomalies or unaccounted phenomena are still being explored from large scale datasets. The reader is preemptively referred to Fig. 19 where the high frequency features related to turbulence cannot be disregarded. Furthermore, for classification and clustering problems where feature dimension reduction is desirable it is imperative that a low-rank approximation of a dataset contains most (≥ 90%) of the original information content without altering the subspace information. In this case, an over-sampled rank can exceed the feature dimension of the data itself (see Section 4.2).

heading: 1.1 PROBLEM STATEMENT
text: Let us denote the raw data matrix as X ∈ Rm×n of rank f ≤ min(m,n) and its approximation as Xr ∈ Rm×n. The SVD of X is X = UΣV T , where U ∈ Rm×f = [u1, · · · , uf ] and V ∈ Rn×f = [v1, · · · , vf ] are its left and right singular vectors respectively, and Σ ∈ Rf×f = diag(σ1, · · · , σf ) are its largest non-zero singular values. The rank r (r ≤ f) truncation of X is then Xr = UrΣrV Tr , where Σr = Σ[1:r] are the top r singular values, and Ur = U[1:r] and Vr = V[1:r] are the left and right singular vectors. In other words, X = UΣV T = UrΣrV Tr + Uf\rΣf\rV T f\r = Xr + Xf\r. Here, Uf\r, Vf\r are the trailing f − r left and right singular vectors. Theorem 1. Eckart-Young-Mirsky Theorem (Eckart & Young, 1936; Mirsky, 1960): Let X ∈ Rm×n be a real, rank-f, matrix with m ≥ n with the singular value decomposition as X = UΣV T , where the orthonormal matrices U, V contain the left and right singular vectors of X and Σ is a diagonal matrix of singular values. Then for an arbitrary rank-r, r ≤ f matrix Br ∈ Rm×n,
‖X −Br‖F ≥ ‖X −Xr‖F where Xr = UrΣrVr with Σr is the diagonal matrix of the largest r singular values and Ur, Vr are the corresponding left and right singular vector matrices.
The problem statement is then: Given X ∈ Rm×n find X̂ such that, arg min
X̂∈Rm×n, rank(X̂)≤r ‖X − X̂‖F (1)
In effect, the minimizer X̂∗ of the above problem gives us the rank-r approximation of X such that Xr = X̂∗. In this work we utilize the minimizer of the above problem to extract the top rank-r SVD factors of X without loading the entire data matrix into the main memory. Note that the minimizer naturally gives the lower bound on this tail energy in addition to being a rank-r approximation.

heading: 1.2 MAIN CONTRIBUTIONS
text: Data and Representation Driven Neural SVD: The representation driven network loss terms ensures that the data matrix X is decomposed into the desired SVD factors such that X = UΣV T . In the absence of the representation enforcing loss term, the minimizer of Eq. 1 results in an arbitrary decomposition such that X = ABC different from SVD factors. A Deterministic Approach with GPU Bit-precision Results: The network can be initialized with weights drawn from a random distribution where the iterative minimization is deterministic. The streaming order of the samples is of no consequence and the user is free to choose the order in which the samples are processed in a batch-wise manner (indexed or randomized). First Streaming Architecture with Exact Low Memory Cost: Range-net requires an exact memory specification based upon the desired rank-r and data dimensions X ∈ Rm×n given by r(n+ r), independent of the sample dimension m. This is the first streaming algorithm that does require the user to wait until the streaming step is complete, contrary to randomized streaming algorithms. Layer-wise Fully Interpretable: Range-Net is a low-weight, fully interpretable, dense neural network where all the network weights and outputs have a precise definition. The network weights are placeholders for the right (or left) orthonormal vectors upon convergence of the network minimization problems (see Appendix D). The user can explicitly plug a ground truth solution to verify the network design and directly arrive at the tail energy bound.

heading: 2 RELATED WORKS
text: The core idea behind randomized matrix decomposition is to make one or more passes over the data and compute efficient sketches. They can be broadly categorized into four branches: 1) Sampling based methods (Subset Selection (Boutsidis et al., 2014) and Randomized CUR (Drineas et al., 2006)); 2) Random Projection based QR (Halko et al., 2011); 3) Randomized SVD (Halko et al., 2011); and 4) Power iteration methods (Musco & Musco, 2015). The sketches can represent any combination of row space, column space or the space generated by the intersection of rows and columns (core space). However, all of these methods require loading the entire data in memory. Readers are referred to Kishore Kumar & Schneider (2017); Ye et al. (2019) for an expanded survey.
Conventional SVD although deterministic and accurate, becomes expensive when the data size increases and requires r passes over the data (see Table 1). The two branches of interest to us are the randomized SVD and Power iteration Methods for extracting SVD factors. Randomized SVD algorithms (Halko et al., 2011) are generally a two stage process: 1) Randomized Sketching uses random sampling to obtain a reduced matri(x/ces) which covers any combination of the row, column and core space of the data; and 2) Deterministic Post-processing performs conventional SVD on the reduced system from Randomized Sketching stage. These approaches make only one pass over the data assuming that the singular value spectrum decays rapidly.
Power iteration based approach (Musco & Musco, 2015) requires multiple passes over the data and are used when the singular value spectrum decays slowly. This class of algorithm constructs a Krylov matrix inspired by Block Lanczos (Golub & Underwood, 1977) to obtain a polynomial series expansion of the sketch. Although these algorithms achieve lower tail-energy errors, they cannot be used in big-data applications when X itself is too large to be retained in the main memory. Here, constructing a Krylov matrix with higher order terms such as (AAT or ATA) is not feasible1.
Due to main-memory restrictions on remote compute machines, streaming (Clarkson & Woodruff, 2009; Liberty, 2013) algorithms became popular. For low-rank SVD approximations these involve streaming the data and updating low-memory sketches covering the row, column and core spaces. Existing randomized SVD capable of streaming include Halko et al. (2011); Upadhyay (2016); Tropp et al. (2017a; 2019), each with different sketch sizes and upper bounds on approximation errors (Table 1).
SketchySVD (Tropp et al., 2019) is the state of the art streaming randomized SVD, with sketch sizes comparable to it’s predecessors and tighter upper bounds on the tail energy and lower errors. As a two stage approach, SketchySVD (Alg. 1) constructs an overestimated rank-(k, s) sketch of the data based on row, column and core projections. A QR decomposition on the row and column sketches gives an estimate of the rank-k subspace. This is followed by a conventional SVD on the core matrix to extract it’s singular values and vectors. Finally, the singular vectors are returned after projecting them back to the original row and column space. The time cost of SketchySVD is O(k2(m + n)) with memory cost O(k(m+ n) + s2) with oversamling parameters k = 4r + 1 and s = 2k + 1.

heading: 2.1 LIMITATIONS OF RANDOMIZED SVD APPROACHES
text: We discuss a few limitations that led us to reformulate the problem in the spirit of EYM theorem. The reader is referred to Appendix A for a detailed discussion and supporting numerical examples to further enunciate these limitations.
Tall and Skinny Matrices: For a rank-r approximation of X ∈ Rm×n of rank-f , Randomized SVD methods rely upon rank-k (k > r) sketches of X . However, these methods are useful only when k ≥ f but for practical datasets f ≤ min(m,n) and therefore the memory requirement can still be overbearing. The reader is referred to Section 4.2 for a practical example. Exponential Decay of Singular Values: Assuming exponential decay implies the rank of X itself is such that f min(m,n). For real world applications, the data matrices are almost full rank f ≤ min(m,n), where a rank r truncation is chosen such that the desired dominant features are accounted for. Appendix A shows a synthetic case with non-exponential decay of singular values where sketching accrues substantial errors. Further, it is difficult to assume that the decay rate will follow a strict functional form: mixture of linear, exponential and others (see Fig. 3). Upper Bound on Tail Energy: The problem statement in Eq. 1 suggests finding a minimum with the minimizer providing a lower bound on the tail-energy. Even if the solution scheme is upper bounded (Halko et al., 2011; Upadhyay, 2016; Tropp et al., 2017a; 2019), the minimizer X̂∗ in Eq. 1
1Readers are referred to Fig. 5 (c) and Fig. 8 for a performance comparison between power iteration schemes and Range-Net for a mid-sized real and a synthetic dataset that can be loaded on our compute machine.
or equivalently achieving the lower bound is necessary. Approximation Errors: A low-rank SVD solver that does not iteratively compute the projection (left or right) while solving Eq. 1 cannot extract SVD factors with low errors even with multiple passes over the data matrix or multiple runs. As shown later in Theorem 2, any subspace projection (left or right) of the original data matrix that does not correspond to the minimizer in Eq. 1 increases the tail energy and therefore results in incorrect low-rank SVD factors (singular values and vectors). Memory Requirement: Randomized SVD requires an optimal choice of hyper-parameters (sketch sizes etc.) that are subjective to the dataset being processed. In a practical, limited memory scenario, this entails tuning the hyper-parameters for optimal trade-off between memory requirement, compute time and an approximation error that does not violate the upper bound on the tail energy.
Remark. A low relative error in tail energy does not imply the extracted singular values and vectors will have similar relative errors at scale. The issue has been raised by Musco & Musco (2015), that merely upper bounding the tail energy equipped with Frobenius or spectral norm does not bound the approximation errors in the extracted singular values or vectors. Therefore, for a fair comparison we show error metrics on the extracted singular factors for all our numerical experiments.

heading: 3 RANGE-NET: A 2-STAGE SVD SOLVER
text: We present Range-Net that explicitly relies upon solving the minimization problem in Eq. 1 to achieve the lower bound on the tail-energy for a desired rank-r approximation of a data matrix X under a streaming setting. The readers are referred to Appendix B for the preliminaries followed by proofs of theorems and lemmas associated with each of the two stages.

heading: 3.1 NETWORK ARCHITECTURE
text: The proposed network architecture is divided into two stages: (1) Projection, and (2) Rotation, each containing only one dense layer of neurons and linear activation functions with no biases. Fig. 2 shows an outline of the this two-stage network architecture where all the weights and outputs have a specific meaning enforced using representation and data driven loss terms. Contrary to randomized SVD algorithms the subspace projection (Stage 1) is not specified preemptively (consequently no assumptions) but is computed by solving an iterative minimization problem following EYM theorem corresponding to Eq. 1. The rotation stage (Stage 2) then reuses the EYM theorem again in a modified form to extract the singular vectors and values.
Stage 1: Rank-r Sub-space Identification: The projection stage constructs an orthonormal basis that spans the r-dimensional sub-space of a data matrix X ∈ Rm×n of an unknown rank f ≤ min(m,n). This orthonormal basis (Ṽ ) is extracted as the stage-1 network weights once the network minimization problem converges to a fixed-point. The representation loss ‖Ṽ T Ṽ − Ir‖F in stage-1 enforces the orthonormal requirement on the projection space (even when r > f ) while the datadriven loss ‖X − XṼ Ṽ T ‖F minimizes the tail energy. Although the minimization problem is non-convex, the tail-energy is guaranteed to converge to the minimum at machine precision. The reader is referred to Appendix C for a discussion on the minimization problem (bi-quadratic loss function with 2r global minima) for details regarding the convergence behavior.
Theorem 2. For any r, f ∈ Z+, 0 < r ≤ f , if the tail energy of a rank-f matrix X ∈ Rm×n, f ≤ min(m,n), with respect to an arbitrary rank-r matrix Br = XṼrṼ Tr is bounded below by the tail
energy of X with respect to it’s rank-r approximation Xr = XVrV Tr as, ‖X−Br‖F ≥ ‖X−Xr‖F where, Vr = span{v1, v2, · · · , vr} and vis are the right singular vectors corresponding to the largest r singular values then the minimizer of arg min
Ṽ ∈R(n×r) ‖X −XṼ Ṽ T ‖F is V∗ such that V∗V T∗ = VrV Tr .
As per Theorem 2, the equality holds true only when span{Br} = span{Xr}. Further, we define Br as, Br = XṼ Ṽ T then, V∗V T∗ = VrV T r and V T ∗ V∗ = Ir where Vr is a rank-r matrix with column vectors as top-r right singular vectors of X . The minimization problem then reads,
min Ṽ
‖X −XṼ Ṽ T ‖F s.t. Ṽ T Ṽ = Ir (2)
with a minimum at the fixed point V∗ = span{v1, . . . , vr} where vi=1,...,r are the right singular vectors of Xr. This minimization problem describes the Stage 1 loss function of our network architecture. Upon convergence, the minimizer Ṽ∗ is such that V∗V T∗ = VrV T r following Theorem 2 where Vr is the matrix with columns as right singular vectors of X corresponding to the largest r singular values of X . Lemma 2.1. If V Tr Vr = Ir and VrV Tr = V∗V T∗ then V T∗ V∗ = Ir. Lemma 2.2. If X ∈ Rm×n is a rank f matrix, then for any rank r > f , where {r, f} ≤ min(m,n), if V T∗ V∗ = Ir and V∗V T ∗ = VrV T r then V T r Vr = Ir. Remark. Note that for r ≤ f , the orthonormality constraint is trivially satisfied as shown in Lemma 2.1. However for r > f , the orthonormality constraint ensures that the column vectors in V∗ are orthonormal (see Lemma 2.2) allowing us to extract orthonormal right singular column vectors of Vr from the Stage 2 minimization problem.
Stage 2: Singular Value and Vector Extraction: The rotation stage then extracts the singular values by rotating the orthonormal vectors (V∗) to align with the right singular vectors (Vr = V∗Θr). From the fixed point of the Stage-1 minimization problem Eq. 2 we have V∗V T∗ = VrV T r . According to the EYM theorem the tail energy of a rank-r matrix XV∗Cr, where Cr is an arbitrary rank-r, real valued, square matrix, with respect to XV∗ is now bounded below by 0 or ‖XV∗ −XV∗Cr‖F ≥ 0. Theorem 3. Given a rank-r matrixXV∗ ∈ Rm×r and an arbitrary, rank-r matrixC ∈ Rm×r, following Theorem 1, the tail energy of XV∗ with respect to XV∗C is bounded as, ‖XV∗ −XV∗C‖F ≥ 0, where the equality holds true if and only if C = Ir. Lemma 3.1. If C = ΘrΘTr , where Θr ∈ Rr×r is a rank-r matrix such that C = Ir, then Θr is a real-valued unitary matrix in an r-dimensional Euclidean space. Theorem 4. Given a rank-r matrix XV∗ ∈ Rm×r, such that V∗V T∗ = VrV Tr where Vr is a matrix with column vectors as the top-r right singular vectors of X , and a real-valued unitary matrix Θr ∈ Rr×r then (XV∗Θr)T (XV∗Θr) is a diagonal matrix Σ2r where Σ2r = diag(σ21 , σ22 , · · · , σ2r) and σis are the top-r singular values of X if and only if V∗Θr = Vr.
From Theorem 3 and Lemma 3.1 we know that, Cr = ΘrΘTr , where Θr is a rank-r, unitary matrix in an r-dimensional Euclidean space. Further, from Theorem 4 we have that (XV∗Θr)T (XV∗Θr) is a diagonal matrix Σ2r = diag(σ 2 1 , · · · , σ2r), where σis are the top-r singular values of X if and only if V∗Θr = Vr. Assuming Y = XV∗Θr, the minimization problem now reads:
min Θr
‖YΘTr −XV∗‖F s.t. Y TY − diag(Y TY ) = 0 (3)
Remark. Note that stage 1 can be verified numerically independently of stage 2 by checking whether the orthonormality condition is met in addition to minimization problem converging to the tail-energy bound. Similarly, stage 2 minimization problem will return a rotation matrix Θr (ΘTr Θr = ΘrΘ T r = Ir, det(Θr) = ±1) upon convergence that can again be verified numerically.
As discussed previously, this choice of loss terms equipped with a Frobenius norm ensures a rank-r approximation in accord with the Eckart-Young-Mirsky (EYM) theorem. We therefore state that the expected values of the stage-1 loss term at the minimum must correspond to the rank (n − r) tail energy. Further, the stage-2 loss is expected to reach a machine precision zero at the minimum. Once, the network minimization problems converge, the singular values are extracted from Stage 2 network weights Θr as Σ2r = (XV∗Θr)
T (XV∗Θr). The right singular vectors can now be extracted using Stage 2 layer weights given by Vr = V∗Θr. Once Vr and Σr are known, left singular vectors Ur = XV∗ΘrΣ−1r . Please note that for r > f , f − r singular values are zero and therefore Σ−1r implies inverting the non-singular values using a threshold of 10−8.
Implementation Details: A detailed discussion and justification of our implementation choices are given in Appendix E. These include a) Activation: all the activation functions are linear with no biases, since SVD requires linearly separable orthogonal features; b) Data Split: we do not perform any split of the training data, since SVD factors are unique to the full data only; c) Data Streaming: we stream data from HDD into main memory to avoid data sample load; d) Training and Setup; e) Error Metrics; and f) Loss Profiles. Note that Range-Net has no hyper-parameters and therefore does not require any post-hoc tuning or adjustments.

heading: 4 RESULTS
text: We present numerical experiments for three datasets: (a) Parrot image, (b) MNIST, and (c) hurricane Sandy. The reader is referred to Appendix E.5 for the definitions of the error metrics used. Note that, these error metrics rely upon conventional SVD as the baseline for a fair comparison. For additional numerical experiments on sparse graph datasets and other low rank approximations see Appendix F. SketchySVD’s algorithmic implementation (Alg. 1) can be found in Appendix G.

heading: 4.1 IMAGE COMPRESSION: PARROTS (SVD)
text: As an example for SVD of natural images, we use the well known Parrots image from the image processing domain. The original image is in an RGB format, converted to a gray scale for demonstration purposes followed by normalization between [0, 1]. This 1024× 1536 data matrix is then used to compute a rank r = 20 approximation for comparison and numerical analysis. Fig. 4 shows the result of the low rank reconstruction for SVD, SketchySVD and Range-Net. Visually one can verify that Fig. 4 (b), (d) are similar while Fig. 4 (c) is different. To make the error in approximation more clear, we plot the absolute difference of SketchySVD and our net from the truncated rank image. Fig. 4 (e), (f) shows the the corresponding plots with heatmaps imposed for clarity. Notice that while the reconstruction error for our network (≈ 10−7) is close to the GPU precision, SketchySVD has significantly higher error scale (≈ 10−1), validating the artifacts in the approximated image. The singular value spectrum does not decay exponentially (Fig. 3) and the data matrix is near-full rank (f ≈ 1024). Fig. 5 (a, b) shows the scree error as the absolute difference between the predicted and the true singular values. For SketchySVD, the error fluctuates across the top r = 20 values, but also the scale of fluctuations is around 1. Comparably, Range-Net incurs significantly lower errors in singular values (scale of 10−4). Fig. 5 (c) shows the reconstruction errors in Frobenius norm for SketchySVD (Tropp et al., 2019) (red line), Block Lanczos with Power Iteration (Musco & Musco, 2015) (black line), Sklearn’s randomized SVD (skr) implementation Halko et al. (2011) with (solid
cyan line) and without power iteration (dashed blue line), and Range-Net (green line) over 1000 runs on the Parrot data. This shows that in order to gain lower reconstruction errors a power iteration is necessary that quickly becomes intangible in a big-data setting. Further, note that the expected error (upper bound) over multiple runs of Randomized SVD algorithms does not contract (reduce).
Fig. 6 shows the cross-correlation between extracted right singular vectors from SketchySVD (left) and Range-Net (right) against conventional SVD for a rank-20 approximation. SketchySVD oversampled rank is k = 81 and still the extracted right singular vectors deviate substantially. This implies that the extracted vectors from SketchySVD do not span the top rank-20 subspace of X as opposed to RangeNet where stage-1 explicitly ensures this span without any oversampling. The higher the vector index, the higher the spread owing to random projections. Range-Net has a near-perfect cross-correlation with the true vectors, indicated by the solid diagonal and zero off-diagonal (near GPU-precision) entries.
We tabulate error metrics (see Appendix E for definitions) for SketchySVD and Range-Net for various ranks in Table. 2. Note that all errors are reported w.r.t. the true SVD, where errfr and errsp denote the frobenius and spectral errors respectively. One can easily see that as the rank increases, SketchySVD’s performance keeps on deteriorating for the χ2err metric (measure of correlation between true and estimated vectors), also evident from Fig. 6. Range-Net on the other hand has consistent lower errors in all metrics while simultaneously being memory efficient.

heading: 4.2 DIMENSION REDUCTION: MNIST (EIGEN / PCA)
text: Principal Component Analysis (PCA) is a special variant of Eigen decomposition, where the samples are mean corrected before constructing a feature covariance matrix followed by Eigen decomposition. Note that Range-Net does not require construction of the feature covariance matrix and can directly extract the eigenvectors and values without any modification. This is due to the fact that the Eigenvalues are the square of the singular values for any non-square data matrix X ∈ Rm×n where right singular vectors are exactly the same as the eigenvectors.
MNIST has 60k images of size 28× 28 in the training set. We reshape each image into a 784-dim vector to obtain the data matrix X ∈ R60000×784, as a tall skinny matrix. In a streaming setting, the mean feature vector computation requires one pass over the data matrix. This can be subsequently used during the network training (Stage-1) to mean correct streamed input vectors.
For this dataset, it is well known that r = 200 captures ≥ 90% variance in the dataset. For SketchySVD, this results in projection matrices of ranks k = 4r + 1 = 801 and s = 2k + 1 = 1602. Since MNIST only has n = 784 features, SketchySVD (Alg. 1) is almost equal if not more memory intensive than conventional SVD for such tall and skinny matrices. Table. 3 shows the error metrics under different rank setting, where even with oversampling, SketchySVD’s errors are high. Range-Net on the other hand with an exact memory requirement (r(n+ r)) can handle much larger full rank tall and skinny matrices without incurring extraneous memory load. As discussed before, since this data matrix is tall and skinny (60k× 784) we already know that for SketchySVD any rank-r s.t. (r ≥ 196) will result in the oversampling parameters k ≥ 784 and s ≥ 1569. SketchySVD will now extract lower-error SVD factors since the oversampled rank redundantly exceeds the feature dimension.

heading: 4.3 SCIENTIFIC COMPUTING: SANDY BIG DATA (SVD)
text: Satellite data gathered by NASA for Hurricane Sandy over the Atlantic ocean represents the big data counter-part for scientific computations. The data-set is openly available2 and comprises of RGB snapshots captured at approximately one-minute interval. The full data-set consists of 896 × 719 pixel images for 1208 time-instances is of size 24 GB. We chose this particular big data so that a conventional SVD can be performed on our machine (16 GB RAM) for benchmarking. Please note that this restriction is imposed by conventional SVD method due to its high main memory requirements. In contrast, our neural SVD solver can handle data sets that are orders of magnitude larger in size with the same hardware specification.
Range-Net can not only handle larger datasets than SketchySVD, but also ensures lower errors in approximating the SVD factors. Tab. 4 shows the error metrics for Range-Net with a comparison of peak main-memory load between SketchySVD and Range-Net for ranks r = [10, 50, 100]. Appendix F.3 shows a comparison of dynamic mode reconstructions obtained from SketchySVD, conventional full rank SVD, and Range-Net and the associated scree errors in the computed singular values.
Randomized SVD extracted factors deviate quite substantially when the user specified rank r is such that the oversampled rank k is much lower than the unknown rank f of a given data matrix. The reader is also referred to the additional experiment in Appendix F.4 where a low rank (r = 10) approximation is extracted. Here, SketchySVD deviates quite substantially after rank-3 while RangeNet still remains in excellent agreement with the baseline singular vectors and values.
2https://www.nasa.gov/mission_pages/hurricanes/archives/2012/h2012_Sandy.html

heading: 4.4 STORAGE COMPLEXITY ANALYSIS
text: To estimate the memory efficiency of Range-Net, let us consider the peak main memory (RAM) requirement for the compuation of SVD factors. Range-Net has two layers in succession, one corresponding to the low-rank projector Ṽ r×n and the rotation matrix Θr×r. For Sketchy SVD, the peak memory load occurs during the construction of a core matrix Cs×s (see Alg. 1). This requires that the two projection matrices Φm×s,Ψn×s, one projected data matrix Zs×s, and two rank-k decomposition Qm×k, Pn×k and the core matrix Cs×s, be present in the memory simultaneously. The memory efficiency factor (seff ) for a rank-r approximation with k = 4r + 1, s = 2k + 1 is:
seff = SketchySVD Range-Net = ns+ms+ 2s2 +mk + nk rn+ r2 = (m+ n)(k + s) + 2s2 r(n+ r) ≈ 12(m+ n) + 128r (n+ r)
≈ 7.67e2 for MNIST(m = 60k, n = 784, r = 200)
To validate the ratio, we constructed a synthetic dataset of m = 50k rows and the number of columns were varied starting at n = 10k with increments of 10k. The expected rank was held at r = 200. Fig. 7 shows the memory allocation (in MegaBytes (MB)) of SketchySVD vs. Range-Net, while n varies between [10k − 150k]. When m = 50k, n = 150k and r = 200, SketchySVD has a peak memory consumption of 14GB due to oversampling parameters of k = 801, s = 1603, while Range-
Net only requires 916MB. Since Range-Net has an exact memory requirement of r(n + r) for a rank-r SVD, it always occupies main lesser memory than oversampling of SketchySVD (or any other randomized method) by two orders of magnitude, and scalable to larger datasets.

heading: 5 DISCUSSION
text: From the numerical experiments it is evident that Range-Net performs at par with conventional SVD with GPU bit precision results, while being extremely memory efficient. Range-Net constructs the rank-r projector V V T iteratively to reach the tail energy lower bound given by EYM. Alternatively, the rank-r minimizer V of Eq. 1 gives the correct rank-r right projector V V T of X . Any arbitrary (random) projection of X onto an oversampled rank-k subspace (k > r), rather than the tail energy minimizing subspace V V T , can inadvertently annihilate the desired top ranking singular features resulting in irreducible approximation errors.
This issue is specially important in exploratory data analysis for scientific computing, where one is not only interested in the top-most singular values of the dataset but also the dominant phenomena (singular vectors). Furthermore, with increasing digital sensor resolution the focus now is to isolate and study the lower energy (high frequency/ small spatial scale) features as in the case of Hurricane Sandy, where turbulence manifests as low rank features (see Figs. 19, 24 in Appendix).
We again point out that accuracy is gained by achieving (or finding) this minimizer or lower bound of the rank-r tail-energy of X and therefore the upper bound on this tail-energy is of no consequence to Range-Net. As pointed out by Musco & Musco (2015), upper bounding the tail energy does not ensure that the approximation errors in the extracted singular values or vectors will reduce. Therefore, small relative tail energy errors can mislead the reader that the singular factors are equivalently accurate wherein the absolute errors can still be substantially large.

heading: 6 CONCLUSION
text: We present Range-Net as a low-weight, high-precision, fully interpretable neural SVD solver for big data applications that is independently verifiable without performing a full SVD. We show that our solution approach achieves lower errors metrics for the extracted singular vectors and values compared to Randomized SVD methods. A discussion is also provided on the limiting assumptions and practical consequences of using Randomized SVD schemes for big data applications. We also verify that our network minimization problems converges to the EYM tail energy bound in Frobenius norm at machine precision. A number of practical problems are considered, where SVD or Eigen decompositions are required, that demonstrate the applicability of Range-Net to large scale datasets. A fair comparison is provided against the state of the art randomized, streaming SVD algorithm with conventional SVD solution as the baseline for benchmarking and verification.

heading: A NEED FOR RANGE-NET
text: To illustrate the limitations of current streaming randomized SVD approaches, we consider a synthetic data matrix X with slow decay in the singular value. The numerical results section later presents a number of these singular value spectra for different practical datasets (Fig. 3) to demonstrate that the decay rates are subjective to the problem at hand.
X = diag(450, 449, · · · , 2, 1︸ ︷︷ ︸ f=450 , 0, · · · , 0︸ ︷︷ ︸ n−f )
Here X ∈ Rm×n is a strictly diagonal matrix with m = n = 500 with rank f = 450, where the singular value spectrum decays linearly. Fig. 8 shows a comparison of reconstruction errors (see Metrics in Appendix E.5 for the definition) for SketchySVD (Tropp et al., 2019) (red line), Block Lanczos with Power Iteration (Musco & Musco, 2015) (black line), Sklearn’s randomized SVD (skr) implementation (Halko et al., 2011) with (solid cyan line) and without power iteration (dashed blue line), and Range-Net (green line) over 1000 runs for this synthetic dataset.
Please note that although power iteration improves the reconstruction error for both Block Lanczos (Musco & Musco, 2015) and Sklearn’s RandSVD (Halko et al., 2011), power iteration itself requires a persistent presence of the data matrix X in the main memory. For a practical big data scenario, power iteration is therefore not a feasible alternative when the data matrix X or it’s sketch is itself too big to be loaded into the main memory. Note that the error expectation (upper bound) over multiple runs of Randomized SVD algorithms does not reduce. We further identify the following requirements for SketchySVD to return SVD factors with relatively lower approximation errors:
1. Decay rate of singular values of a dataset must be exponential.
2. For a rank-f matrix, the desired rank r must be chosen such that the oversampled rank k is strictly greater than f (k ≥ f) to achieve lower errors at scale compared to other runs.
We suggest that the reader also attempt the case where all the diagonal entries are strictly ones and zeros under a high rank setting.
Fig. 9 shows the singular values extracted by SketchySVD for a linearly decaying spectrum with corresponding errors in absolute and relative tail energies. The reader is referred to Appendix
B.1 for the definitions of tail energy and relative tail energy. Note that the synthetic data is a diagonal matrix chosen specifically so that the exact tail energies can be computed using Frobenius norm as (
k∑ i=r+1 x2ii
) 1 2
. For a rank-r approximation, SketchySVD suggests oversampling by a
factor of k = 4r + 1 to extract the rank-r factors correctly. Hence for an oversampled rank k = [41, 81, 121, 161, 201] the corresponding top rank r = [10, 20, 30, 40, 50] extracted singular values and vectors will have the lowest approximation errors. However as shown in the Fig. 9 (a) the extracted singular values have an order of magnitude difference w.r.t. the ground truth. Consequently, Fig. 9 (b) shows that the absolute tail-energies of the extracted features deviate quite substantially from the true tail-energy. Furthermore, we also notice that the deviations remains large as long as the oversampled rank-k is such that k < f . For a practical dataset f is either unknown or almost full rank f ≤ min(m,n) or both and can only be detected by performing a full SVD on the dataset. This poses a serious restriction on SketchySVD’s reliability for a realistic big data application, due to an exponential decay assumption.
We also notice that for smaller values of r, the accrued error in both the extracted singular values and tail energy error is worse. Fig. 9 (b) shows that for different rank approximations SketchySVD tail energy deviates from the truth quite substantially. This is due to the fact that the oversampled rank k = 4r + 1 f , as pointed out before. Note that for oversampling parameter k < f , although the error decreases as k → f the memory requirement increases correspondingly asO(k) for extracting a low rank-r approximation. This implies that for slow decaying spectrum optimal values of k are such that k ≥ f even when a rank-1 approximation is desired. In Fig. 9 (c), one can observe relative errors between 10−2 and 100. Although this implies that the actual rank-r tail-energy approximation error is off by 1% at the best, the extracted singular values and vectors are off by one order of magnitude. As a consequence, the extracted singular vectors no longer represent the features of the dataset. Remark. For a randomized SVD algorithm to converge (without power iterations) to a rank-r approximation Xr over multiple runs, we posit that a rank-r sketch matrix X̂ for a given rank-f dataset X , for f ≥ r, be such that P (span{X̂} ∩ span{Xr} = span{Xr}) ≥ 0.5. However, ensuring this requires substantial amount of prior knowledge or intelligent sampling (a multipass iterative algorithm).
Range-Net with it’s explicit minimization of tail energy is capable of intelligent sampling on an arbitrary matrix without requiring any prior information. The key point to note here is that Range-Net relies upon an iterative computation of a near optimal projector instead of arbitrary/user-specified projectors used in Randomized SVD schemes. Even if the tail energy is theoretically upper bounded for some of the Randomized SVD schemes, the target is to find the lower-bound (minimizer) on the tail energy as discussed in Section 1.1. Furthermore, since none of the randomized SVD schemes construct the projector in an iterative manner while minimizing Eq. 1, the relative error in the tail energy remains high. Even if multiple runs of SketchySVD or Sklearn’s RandSVD are performed, the reconstruction errors in tail energies remain the same at scales shown in Fig. 8. We would also like to point out that although one must strive for lower errors (relative or otherwise) and tighter theoretical upper and lower bounds, in practice we should also closely monitor if these theoretical bounds deliver us the desired solution.

heading: B THEORETICAL GUARANTEES

heading: B.1 PRELIMINARIES
text: The Frobenius norm of a matrix A is given by,
‖A‖F =
 ∑
i
∑
j
a2ij
  1 2 = ( Tr(ATA) ) 1 2 = ( Tr(AAT ) ) 1 2
Further, the Frobenius norm can be used to bound Trefethen & Bau III (1997) the norm of a matrix product as,
‖AB‖F ≤ ‖A‖F ‖B‖F
For a Frobenius norm we have that,
‖A+B‖F ≤ ‖A‖F + ‖B‖F ‖A−B‖F ≥ ∣∣‖A‖F − ‖B‖F ∣∣
Also the Frobenius norm of a rank-f matrix A,
‖A‖F = ‖Σf‖F where Σf = diag(σ1, σ2, · · · , σf ) and σis are the f non-zero, singular values of A. Let A, B, C be matrices such that the following matrix products are feasible. The cyclic property of the linear trace operator is,
Tr(ABC) = Tr(BCA) = Tr(CAB)
Definition 1. The tail energy of an arbitrary matrix B ∈ Rm×n with respect to a given matrix X ∈ Rm×n equipped with a Frobenius norm is defined as,
τ = ‖X −B‖F Definition 2. Let r, f ∈ Z+ be positive integers such that 0 < r ≤ f , then a rank-r truncation Xr of a rank-f matrix X is defined as,
Xr = UrΣrV T r = XVrV T r
where, Σr = diag(σ1, σ2, · · · , σr) and σis are the top r singular values of X and Vr = [v1, v2, · · · , vr], and Ur = [u1, u2, · · · , ur] are matrices such that vis and uis are the corresponding right and left singular vectors, respectively.
The relative tail energy of a rank-r matrix Br with respect to a rank-f matrix X (r ≤ f ) is then defined as,
τrel,r = ‖X −Br‖F ‖X −Xr‖F − 1

heading: B.2 STAGE 1
text: Proof of Theorem 2. For any r, f ∈ Z+, 0 < r ≤ f , if the tail energy of a rank-f matrix X ∈ Rm×n, f ≤ min(m,n), with respect to an arbitrary rank-r matrix Br = XṼrṼ Tr is bounded below by the tail energy of X with respect to it’s rank-r approximation Xr = XVrV Tr as,
‖X −Br‖F ≥ ‖X −Xr‖F where, Vr = span{v1, v2, · · · , vr} and vis are the right singular vectors corresponding to the largest r singular values then the minimizer of arg min
Ṽ ∈R(n×r) ‖X −XṼ Ṽ T ‖F is V∗ such that V∗V T∗ = VrV Tr .
From Theorem 1 we have,
‖X −Br‖F − ‖X −Xr‖F ≥ 0 (4)
Let Vr = {v1, v2, . . . , vr} be the top-r, right-singular vectors of X corresponding to the largest singular values,
Xr = XVrV T r = UΣV TVrV T r = UΣV T r (5)
Also let Br = XṼrṼ Tr where Ṽr is an arbitrary rank-r matrix. From triangle inequality we have that,
‖X(VrV Tr − ṼrṼ Tr )‖F ≥ ‖X(In − ṼrṼ Tr )‖F − ‖X(In − VrV Tr )‖F (6)
Combining Eq. 4 and Eq. 6 we get,
‖X(VrV Tr − ṼrṼ Tr )‖F ≥ 0 (7)
Additionally, ‖X‖F ‖VrV Tr − ṼrṼ Tr ‖F ≥ ‖X(VrV Tr − ṼrṼ Tr )‖F (8)
Using the above two inequalities we arrive at,
‖X‖F ‖VrV Tr − ṼrṼ Tr ‖F ≥ 0 (9) Since ‖X‖F > 0, equality is achieved when ṼrṼ Tr = V∗V T∗ = VrV Tr . In other words, span{Vr} = span{V∗} since (V∗Θr)(V∗Θr)T = V∗V T∗ for any rank-r, real valued, unitary matrix Θr spanning the top rank-r subspace of X .
Remark. Theorem 2 also implies that any matrix ṼrṼ Tr that does not span the same rank-r subspace of X as VrV Tr will result in a higher tail-energy than given by the EYM tail-energy bound equipped with a Frobenius norm.
Proof of Lemma 2.1. If V Tr Vr = Ir and VrV Tr = V∗V T∗ then V T∗ V∗ = Ir.
Let us assume V Tr Vr = Ir then,
‖V Tr Vr − Ir‖2F = 0 Tr ( V Tr VrV T r Vr ) + Tr (Ir)− 2 Tr ( V Tr Vr ) = 0
Using the cyclic property of the trace operator we have,
Tr ( VrV T r VrV T r ) + Tr (Ir)− 2 Tr ( VrV T r ) = 0
Using Theorem 2, V∗V T∗ = VrV Tr ,
Tr ( V∗V T ∗ V∗V T ∗ ) + Tr (Ir)− 2 Tr ( V∗V T ∗ ) = 0
Again using the cyclic property of the trace operator we now get,
Tr ( V T∗ V∗V T ∗ V∗ ) + Tr (Ir)− 2 Tr ( V T∗ V∗ ) = 0
Hence, V T∗ V∗ = Ir. This shows that following Theorem 2, the matrix V∗ comprises of orthonormal column vectors spanning the same top rank-r subspace of X as the orthonormal column vectors Vr.
Proof of Lemma 2.2. If X ∈ Rm×n is a rank f matrix, then for any rank r > f , where {r, f} ≤ min(m,n), if V T∗ V∗ = Ir and V∗V T∗ = VrV Tr then V Tr Vr = Ir.
‖V T∗ V∗ − Ir‖2F = 0 Tr ( V T∗ V∗V T ∗ V∗ ) + Tr (Ir)− 2 Tr ( V T∗ V∗ ) = 0
Using the cyclic property of the trace operator we have,
Tr ( V∗V T ∗ V∗V T ∗ ) + Tr (Ir)− 2 Tr ( V∗V T ∗ ) = 0
Using Theorem 2, V∗V T∗ = VrV Tr ,
Tr ( VrV T r VrV T r ) + Tr (Ir)− 2 Tr ( VrV T r ) = 0
Again using the cyclic property of the trace operator,
Tr ( V Tr VrV T r Vr ) + Tr (Ir)− 2 Tr ( V Tr Vr ) = 0
Tr ( V Tr VrV T r Vr + Ir − 2V Tr Vr ) = 0
‖V Tr Vr − Ir‖F = 0 Hence V Tr Vr = Ir. Remark. Lemma 2.2 shows that for a rank-r approximation of a rank-f matrix X such that r > f , the extracted right singular vectors are orthonormal when V T∗ V∗ = Ir. This justifies the constraint Ṽ T Ṽ = Ir for the stage-1 minimization problem in Eq. 2 and is numerically verified in Appendix D. Remark. Note that Lemma 2.1 and 2.2 does not imply that V∗ = Vr instead V∗Θr = Vr, where Θr is any real valued unitary matrix for the equality to hold true.

heading: B.3 STAGE 2
text: Proof of Theorem 3. Given a rank-r matrix XV∗ ∈ Rm×r and an arbitrary, rank-r matrix C ∈ Rm×r, following Theorem 1, the tail energy of XV∗ with respect to XV∗C is bounded as, ‖XV∗ −XV∗C‖F ≥ 0 where the equality holds true if and only if C = Ir.
‖XV∗(Ir − C)‖F ≥ 0 ‖XV∗‖F ‖Ir − C‖F ≥ ‖XV∗(Ir − C)‖F
‖XV∗‖F ‖Ir − C‖F ≥ 0 Since XV∗ > 0, this implies equality is achieved if and only if C = Ir.
Proof of Lemma 3.1 If C = ΘrΘTr , where Θr ∈ Rr×r is a rank-r matrix such that C = Ir, then Θr is a real-valued unitary matrix in an r-dimensional Euclidean space.
ΘrΘ T r = Ir
‖ΘrΘTr − Ir‖2F = 0 Tr ( ΘrΘ T r ΘrΘ T r + Ir − 2ΘrΘTr ) = 0
Tr ( ΘrΘ T r ΘrΘ T r ) + Tr (Ir)− 2 Tr ( ΘrΘ T r ) = 0
Using the cyclic property of the trace operator, Tr ( ΘTr ΘrΘ T r Θr ) + Tr (Ir)− 2 Tr ( ΘTr Θr ) = 0
Tr ( (ΘTr Θr − Ir)(ΘTr Θr − Ir)T ) = 0
‖ΘTr Θr − Ir‖2F = 0 This implies ΘTr Θr = Ir. Since Θ T r Θr = ΘrΘ T r = Ir this implies that Θr is a real-valued unitary matrix in the r-dimensional Euclidean space.
Proof of Theorem 4. Given a rank-r matrix XV∗ ∈ Rm×r, such that V∗V T∗ = VrV Tr where Vr is a matrix with column vectors as the top-r right singular vectors of X , and a real-valued unitary matrix Θr ∈ Rr×r then (XV∗Θr)T (XV∗Θr) is a diagonal matrix Σ2r where Σ2r = diag(σ21 , σ22 , · · · , σ2r) and σis are the top-r singular values of X if and only if V∗Θr = Vr.
(XV∗Θr) T (XV∗Θr) = Σ 2 r
ΘTr V T ∗ (X TX)V∗Θr = Σ 2 r
V T∗ X TXV∗ = ΘrΣ 2 rΘ T r
V∗V T ∗ X TXV∗V T ∗ = V∗ΘrΣ 2 rΘ T r V T ∗
Using V∗V T∗ = VrV T r from Theorem 2,
VrV T r X TXVrV T r = V∗ΘrΣ 2 rΘ T r V T ∗
VrΣ 2 rV T r = (V∗Θr)Σ 2 r(V∗Θr) T
(Vr − V∗Θr)Σ2r(Vr − V∗Θr)T = 0 ‖(Vr − V∗Θr)Σ2r(Vr − V∗Θr)T ‖F = 0
Using Frobenius norm to bound the matrix product, ‖Σ2r‖F ‖Vr − V∗Θr‖2F ≥ 0 Since ‖Σ2r‖F > 0, equality is achieved if and only if V∗Θr = Vr. Remark. Note that Θr is a rank-r unitary matrix wherein both rotation (det(Θr) = 1) and reflection (det(Θr) = −1) are valid since the order of the orthonormal vectors in the matrix V∗Θr = Vr do not alter ‖Vr − V∗Θr‖F . In practice, Θr manifests itself predominantly as a rotation matrix during the iterative minimization using gradient descent.

heading: C ENERGY MINIMIZATION, LOSS SURFACE GEOMETRY, AND CONVERGENCE
text: In this section, we consider the energy minimization problem that constructs the projection space spanning the rank-r sub-space of a given data matrix. For ease of visualization, we consider a 2× 2 matrix X = diag(5, 1) with singular values 5 and 1 corresponding to right singular vectors v1 = [1, 0]T and v2 = [0, 1]T , respectively. Our objective here is to extract a rank r = 1 approximation of this rank f = 2 matrix X . Certainly, this corresponds to identifying the right singular vector v1 with singular value 5. The tail-energy surface (log-scale) corresponding to the ‖XṽṽT −X‖F is shown in Fig. 10. Here, ṽ is the test vector for a rank 1 approximation of X. The tail-energy is a bi-quadratic function in ṽ with 1 maximum, 2r minima and 2r saddle points, where r is the desired low-rank approximation of a given data matrix. Furthermore, all minima have the same tail energy: a property of bi-quadratic functions. For the current specific example, the two equal tail-energy minima correspond to v1 and −v1, respectively.
Although the minimization problem is non-convex, convergence is guaranteed since any perturbedgradient descent approach converges to either of the stable fixed points (minima). In other words, any test vector ṽ other than v1 or −v1 will increase the tail-energy and hence will not be the solution. The same argument applies for a high-dimensional dataset X where a low rank-r approximation is desired with the number of equal tail-energy minima corresponding to 2r for all possible negative and positive combinations of the r right singular vectors vi=1,··· ,r. In effect, the stage 1 minimization problem constructs a right projection space Ṽ = span{v1, . . . , vr} that spans the top rank-r subspace of a given dataset X . A similar line of argument then applies to our stage 2 minimization problem as well. A mild limitation, that will be addressed in our future work, occurs when X = diag(5, 5 + ), 0 ≤ << 1, wherein the two right singular values cannot be resolved accurately (still better than Randomized SVD methods) without further considerations. This latter case, with near algebraic multiplicity in singular values is a special case for conventional SVD as well.

heading: D NETWORK INTERPRETABILITY
text: As described before in Fig. 2, our network weights and outputs are strictly defined and incorporated as losses in the network minimization problem. We therefore refer to the problem informed (SVD) restrictions on the network weights as representation driven losses. This is in contrast to kernel regularization loss often considered to impose a weak requirement on the network weights to be small. The representation driven, orthonormality loss term, in Stage 1 enforces that the weights Ṽ must be orthonormal or (V T∗ V∗ = Ir) for a desired rank-r which is greater than the rank-f of the data matrix X . We numerically verify the interpretability of the layer outputs and weights by considering two networks: (1) with, and (2) without the aforementioned orthonomality loss. For each of these two cases, two synthetic datasets are considered corresponding to r ≤ f and r > f . Please note that in a practical scenario f is an unknown and can be determined only by performing a full SVD of X . Therefore, numerically testing this aspect for our solver is necessary.
For the first case, we consider a synthetic data matrix X15×15 where the top 5 singular values are positive (f = 5) while the rest are zero. The objective is to extract the top 10 (r = 10) singular vectors where the desired rank is higher the the rank of the system itself. A total of four training runs
are considered: one run for a network with the orthonormality condition imposed and three runs for another network without this additional constraints. Fig. 11 shows the correlation map between the recovered vectors V∗ for each of the four runs. Notice that only when the orthonormality criteria is not imposed, we get scatter away from the diagonal matrix, although all four runs converged to the same tail energy. Since the true rank of X is 5, the null space of X is of dimension 10. The absence of this orthonormality imposing, representation loss results in non-orthonormal vectors V∗ spanning the low-rank range space.
For the second case, we consider a full-rank, synthetic data matrix X15×15 with f = 15. As before, we extract the top 10 singular vectors (r = 10) using four training runs: one with and three without imposing the orthonormality loss. Since the desired rank 10 system now itself is full rank, this additional loss term does not contribute, as desired. Fig. 12 shows that the extracted vectors V∗ remain orthonormal, for all the four runs, so as to minimize tail energy (‖X(Ṽ Ṽ T − I)‖F ), as described in Section 3.1 above. In fact, for any rank r approximation of a rank f system such that r ≤ f , an arbitrary non-orthonormal matrix V∗ will increase the tail energy and hence will not be a fixed point (solution) of our network minimization problem.
E IMPLEMENTATION DETAILS

heading: E.1 CHOICE OF ACTIVATION FUNCTIONS
text: All the activation functions in both stages are Linear with no biases. One might argue that this choice is not a neural approach, since all the activations are linear. However, please note that singular vectors are linearly separable orthogonal features of a dataset, and therefore any other choice of activation function will result in approximation errors. Since the elements of singular vectors are in [−1, 1], a choice of relu activation is problematic. A simple verification is to approximate a straight line segment in [−1, 1] with tanh activation, only to realize that the approximation error→ 0 as the number of neurons→ ∞. These arguments can also be verified by replacing linear activation in Stage 1 by any non-linear activation only to find that the tail energy bound cannot be satisfied. Note that given a small matrix, one can calculate the right singular vectors and substitute them directly as our network weights to confirm this tail energy bound.

heading: E.2 TRAINING, VALIDATION, AND TESTING SPLIT
text: An issue with training and validation split in matrix decomposition problems is that the error norm cannot be bounded in a deterministic manner or computationally verified. For Singular Value Decomposition of a given data matrix X differs from SVD on a truncated dataset X̂ in it’s singular triplets (singular values and vectors). Ensuring these triplets do not change over an arbitrary split is a non-trivial computational task.
Remark. For datasetX , an arbitrary training/validation split results in a varying dataset X̂ wherein the norm ‖Xpred − X̂‖F changes according to the split. Since the desired features are unknown a priori, a consistent truncated dataset X̂c that spans the same space as the full data X cannot be obtained using an arbitrary split.
This results in a large variance in extracted features over multiple training/testing splits since the span of X̂ itself is changing with each split. Fig. 13 shows a numerical experiment with a synthetic (rank-2) dataset containing 100 samples drawn (red dots) from a 2D ellipse with major and minor axes 2 and 1, respectively. The ground truth right singular vectors over the unsplit dataset are computed using conventional SVD (solid red lines). Next the dataset was split with an 80/20 ratio as the training/testing split for 200 different realizations per run. The mean/expectation of the extracted right singular vectors (dashed green line) over the 200 realizations are then reported. One can easily see that the expected right singular vectors vary quite substantially compared to the ground truth over the three runs indicating that an arbitrary split of the dataset does not ensure accuracy.

heading: E.3 DATA STREAMING
text: Given a data matrix Xm×n, we stream the data along the smaller dimension assuming the user prescribed rank-r is such that r ≤ min(m,n). For the sake of simplicity we assume that the data matrix has m samples and n features, where m > n, and consequently feature vectors of samples are streamed in batches. We rely upon the built-in Keras fit_generator class for data streaming from the secondary memory (HDD). For a big data matrix X that cannot be loaded into the main memory, this allows us to mimic the modality of data residing on an external server. Thus, given a pointer to the data, the function yields a batch of specified size for the network to train on for specific epochs. This ability saves main memory load and allows us to process bigger datasets on smaller main-memory machines than reported in prior works.
Note that for the stage-1 network to converge to a desired tolerance, we require multiple passes (empirically ≤ 5) over the original data streamed batchwise. Therefore for Stage 1, the input data is never persistently present in the main memory of the remote machine. The output data is dumped onto the secondary memory assuming that storing a low rank approximation is still main memory intensive. For Stage 2, this low rank approximation in the secondary memory is streamed as input, and the extracted singular values and vectors are saved in main memory.

heading: E.4 SETUP AND TRAINING
text: All experiments were done on a setup with Nvidia 2060 RTX Super 8GB GPU, Intel Core i7-9700F 3.0GHz 8-core CPU and 16GB DDR4 memory. We use Keras (Chollet, 2015) running on a Tensorflow 2.0 backend with Python 3.7 to train the networks presented in this paper. For optimization, we use AdaMax (Kingma & Ba, 2014) with parameters (lr= 0.001) and 1000 steps per epoch.

heading: E.5 ERROR METRICS
text: As discussed before in Section 2.1, since relative errors in tail energies do not imply similar errors at scales in extracted singular factors, we rely upon additional error metrics on the extractor factors for
performance comparison and benchmarking. In the following X and X̂ are used to denote the true and the reconstructed data matrices.
• Scree Error: Absolute difference between true and approximated singular values.
screeerr(r) = |σi(X)− σ̂i(X̂)| ∀ i ∈ [1, r] • Reconstruction Error: Frobenius norm error of the true data and its rank-r approximation.
froberr(r) = ‖X − X̂‖2F − ‖X −Xr‖2F • Spectral Error: 2-norm of the singular value of the true data and its rank-r approximation.
spectralerr(r) = ‖X − X̂‖2 − ‖X −Xr‖2 • Chi Square Statistic: Deviation between the true and approximated singular vectors.
χ2err(r) = 1− 1
r ‖Corr(v[:r](X), v̂[:r](X̂))‖F
Here, σis are the true singular values and Xr is the desired rank-r approximation of X using conventional SVD as the baseline for benchmarking. Under perfect recovery, all error metrics are expected to approximately achieve zero at machine precision. All of our numerical experiments were performed on a GPU using single (32-bit) precision floating point operations. Therefore, the tail energies are expected to be correct to up to 8 significant digits in all the subsequent calculations. In the following sections, X̂ is replaced by approximations from SketchySVD and Range-Net.

heading: E.6 LOSS PROFILE

heading: F ADDITIONAL EXPERIMENTS

heading: F.1 FEATURE EXTRACTION: GRAPH (EIGEN)
text: Large scale networks occur in many applications where SVD is primarily used to identify the most important nodes or as a pre-processing step for community detection. For these kind of graph based datasets, we either perform SVD or Eigen decomposition on the graph, depending on the format in which the data arrives. We demonstrate results on the following graphs of varying size, tabulated in Table 5. If the data arrives directly in the form of an adjacency matrix, we can perform SVD or Eigen decomposition on it directly. For cases, where an adjacency list is provided, a pre-processing step is required to convert the list representation in a sparse vector.
Since an Eigen decomposition problem is a special case of SVD, where the data matrix is symmetric positive semi-definite, Range-Net is directly applicable. The benchmark was generated for smaller graphs using a conventional SVD solver. For larger graphs, a similar benchmark was constructed using the irlba routine by Baglama & Reichel (2005). Table. 5 shows the error metrics for all the graphs, where consistently low values are observed for Frobenius and Spectral error metrics.

heading: F.2 NAVIER-STOKES SIMULATED DATA (SVD)
text: For our next numerical experiment, we rely upon synthetic data generated using a Navier Stokes flow simulator for an incompressible fluid. The flow data is available on tensor-product grid on two-spatial and one temporal dimensions of size (w, h, t) = (100 × 50 × 200). For each point on the grid, velocity vector values are available in both x and y spatial dimensions for 200 time instances. The fully-developed, flow pattern exhibits a periodicity in the time dimension at approximately every∼ 60 time step that can be identified using SVD as characteristic modes. The data is therefore reshaped into a spatial vector for each time instance resulting in a spatio-temporal matrix X ∈ R5000×200. For comparison purposes, we use only the x-direction stream velocity.
Fig. 15 shows the evolution of the stream velocity over three time instances and the inherent timeperiodic nature of the data. Notice the central-left region of primary flow across all three images. Fig. 16 shows the reshaped U vectors for SketchySVD, conventional SVD and Range-Net. Notice that the images corresponding to the first left singular mode (also called dynamic modes) captures a notion of the primary flow in the left-center part. The second one captures spatial variations of the flow as time progresses. For all the three methods, all the modes have similar solution visually. Note that
Range-Net computes a rank r = 5 approximation without oversampling, wherein SketchySVD relies upon memory intensive sketchy projections of ranks k = 21, s = 43 to arrive at the solution. Overall, for this low-rank dataset both SketchySVD and Range-Net perform reasonably well qualitatively looking at the features in Fig. 16 and the singular value spectrum in Fig. 17, because this synthetic dataset is extemely low-rank (f = 10).
Remark. For a given data matrixX ∈ Rm×n of rank-f (f ≤ min(m,n)) SketchySVD generates low error approximations if the oversampled rank k is such that k ≥ f . For full rank tall skinny matrices, this implies that k ≥ min(m,n). For cases where k < f (full rank or otherwise), SketchySVD accrues large approximation errors resulting in incorrect SVD factors.
From a use-case point of view, randomized SVD generates low-error factors for full-rank, tall skinny matrices (X ∈ Rm×n) only when the oversampled rank k ≥ n. This poses a serious limitation for all applications where this requirement is not met and consequently randomized SVD algorithms accrue large approximation errors in SVD factors as shown in the Sandy Big Data case study below and Appendix F.4. Note that given a big data matrix X determining the rank f of X is unknown and therefore selecting an oversampled rank k such that k ≥ f is impractical in such cases. Remark. Once the SVD factors are extracted, SketchySVD cannot be independently verified without performing a full SVD. In contrast, Range-Net is independently verifiable since Stage-1 of Range-Net cannot return orthonormal vectors if the vectors do not span the rank-r subspace of a given data X .

heading: F.3 HIGH RANK APPROXIMATION: SANDY BIG DATA
text: This section provides an addendum to the numerical results presented in Section 4.3 of the main text. Fig. 18 shows evolution of Hurricane Sandy for two time instances. Similar to the Navier-Stokes simulation data, Fig. 19 shows three dynamic modes corresponding to rank 1, 20, 50, 100 singular values. As shown, our results are in good agreement with conventional SVD whereas, Sketchy SVD shows substantial deviations after the first 50 dynamic modes. Fig. 20 shows the scree-error in the singular values extracted by SketchySVD and Range-Net.
Fig. 21 shows the cross correlation of the right singular vectors and scree-error in the corresponding singular values extracted by SketchySVD and Range-Net. Note that for a rank-100 approximation, SketchySVD extracted right singular vectors start deviating after rank-10 as shown in Figs. 21 while the singular values deviate quite substantially from rank-1. Range-Net on the other hand is in excellent agreement with the right singular vectors and values for all desired 100 indices.
We point out that accuracy is a matter of special concern in scientific computations. Any compression that results in substantial loss of information or obscuring an otherwise identifiable feature in recorded observations directly culls our capacity to make scientific improvements. Consequently, any exploratory data analysis, however big or small, must accurately identify the underlying features. Range-Net achieves the lower bound on tail-energy given by EYM theorem to ensure an accurate resolution in big data setting. Note that increasing the sensor resolution implies that we are interested in exploring and understanding the high-frequency features (lower singular values) of the data.

heading: F.4 LOW RANK APPROXIMATION: SANDY BIG DATA
text: In this experiment, we extract the rank-10 SVD factors using SketchySVD and Range-Net for the Sandy dataset. The oversampled ranks for sketchy SVD are k = 4r + 1 = 41 and s = 2k + 1 = 83 where k, s min(m,n). As before, Fig. 22 show the cross-correlation between the extracted and true (conventional SVD) right singular vectors using SketchySVD and Range-Net. Fig. 23 shows the scree error in the extracted singular values for the two methods with singular values from conventional SVD as the baseline. Finally, Fig. 24 shows a comparison between extracted dynamic modes corresponding to indices i = [1, 4, 7, 10] from SketchySVD, conventional SVD, and Range-Net. One can easily see
that Sketchy SVD extracted dynamic modes/right singular vectors deviate quite substantially for i = [4, 7, 10].

heading: F.5 PARROT: ADDENDUM CORRELATION OF LEFT SINGULAR VECTORS
text: In the following, we show the deviation of left singular vectors, computed using Range-Net and Sketchy SVD, from the conventional SVD computed left singular vectors as a baseline. As shown before, random sketching introduces irreducible errors in randomized SVD methods resulting in this unchecked deviation. Fig. 25 shows a comparison of the cross-correlation against the common baseline for both Range-Net and Sketchy SVD.

heading: G SKETCHY SVD IMPLEMENTATION
text: A brief outline of the single-pass Sketch-SVD algorithm from Tropp et al. (2019). Note that for the numbers reported in terms of storage, we implemented this code with additional memory optimization and sparse matrices.
Algorithm 1 Sketchy SVD Input: X ∈ Rm×n, r : expected rank Output: X̃ ∈ Rm×k the approximated rank k-dim data 1: Initialize k = 4r + 1, s = 2k + 1 . Oversampling parameters 2: Projection maps: Υ ∈ Rk×m,Ω ∈ Rk×n,Φ ∈ Rs×m,Ψ ∈ Rs×n 3: Projection matrices: A ∈ Rk×n, B ∈ Rm×k, Z ∈ Rs×s as empty 4: for i = 1 : n do . Streaming Update 5: Form H ∈ Rm×n as a sparse empty matrix 6: H(i, :) = X(i, :) . Streamed columns 7: A← A+ ΥH . Update Co-Range 8: B ← B +HΩT . Update Range 9: Z ← Z + ΦHΨT . Update Core Sketch 10: Q ∈ Rm×k ← qr_econ(B) . Basis for Range 11: P ∈ Rn×k ← qr_econ(AT ) . Basis for Co-Range 12: C ∈ Rs×s ← ((ΦQ) \ Z)/(ΨP ) . Core Matrix 13: [U,Σ, V T ]← svd(C) . Full SVD of Core Matrix 14: Σ ∈ Rr×r ← Σ[1 : r, 1 : r] . Pick top r 15: U ∈ Rk×r ← U [:, 1 : r] . Pick top r 16: V T ∈ Rr×k ← V T [1 : r, :] . Pick top r 17: U ∈ Rm×r ← QU . Project to Row Space 18: V T ∈ Rr×n ← PV T . Project to Column Space 19: X̃ ∈ Rm×n ← UΣV T . Approximation

Review summaries:
The paper presents a multipass streaming algorithm for rank-r SVD.  Given an input matrix X in R^{mxn} the algorithm identifies two matrices V* in R^{nxr} and H in R^{rxr}. V* has orthonormal columns that span the top r right singular vectors.  H rotates V* so that V*H = Vr, where Vr in R^{nxr} is the matrix of the top-r right singular vectors.  V* and H are computed by minibatch gradient descent, trained until convergence with custom loss functions.  V* is computed in the first stage and then H in the second stage.   The algorithm is accurate and uses very little storage.

The paper presents a streaming method to compute an approximate Singular Value Decomposition, without requiring to load the entire data matrix into the main memory. The first stage of the method identifies a basis optimizing for the constraint of orthonormality, as well as minimizing the error between the input and the reconstruction of the input based on the learnt basis. The second stage corresponds to a rotation step, which aligns the low-rank approximation extracted with the SVD factors.

This paper proposes a neural network architecture for computing a rank-r approximation of a matrix. The first part consisted of learning the projection and the second part consisted of learning the rotation. An experimental comparison against the state-of-the-art Sketchy SVD algorithm is provided on three datasets. The authors claim that the proposed method is extremely memory efficient compared to previous approaches.