The increasing number of applications requiring the solution of large scale singular value problems have rekindled interest in iterative methods for the SVD. Some promising recent ad- vances in large scale iterative methods are still plagued by slow convergence and accuracy limitations for computing smallest singular triplets. Furthermore, their current implementations in MATLAB cannot address the required large problems. Recently, we presented a preconditioned, two-stage method to effectively and accurately compute a small number of extreme singular triplets. In this research, we present a high-performance software, PRIMME SVDS, that implements our hybrid method based on the state-of-the-art eigensolver package PRIMME for both largest and smallest singular values. PRIMME SVDS fills a gap in production level software for computing the partial SVD, especially with preconditioning. The numerical experiments demonstrate its superior performance compared to other state-of-the-art software and its good parallel performance under strong and weak scaling.
The first approach computes eigenpairs of the normal equations matrix @math (or @math if @math ). The eigenpairs of @math correspond to the squares of the singular values and the right singular vectors of @math , @math . If @math , the corresponding left singular vectors are obtained as @math . Although extreme Hermitian eigenvalue problems can be solved very efficiently, @math can be ill-conditioned. Assume @math is a computed singular triplet corresponding to the exact triplet @math with residual satisfying @math . Then the relative errors satisfy the relations @cite_51 , where @math is the closest singular value to @math , @math is a constant that depends on the dimensions of @math , and @math is the machine precision. Since all singular values (except the largest @math ) will lose accuracy, this approach is typically followed by a second stage of iterative refinement. The singular triplets can be refined one by one @cite_18 @cite_20 @cite_28 or collectively as an initial subspace for an iterative eigensolver @cite_12 . This is especially needed when seeking smallest singular values.
@cite_18 @cite_28 @cite_20 @cite_51 @cite_12
2167193748 1918870341 2052926498 1747539222
The harmonic Lanczos bidiagonalization method can be used to compute the smallest singular triplets of a large matrix @math . We prove that for good enough projection subspaces harmonic Ritz values converge if the columns of @math are strongly linearly independent. On the other hand, harmonic Ritz values may miss some desired singular values when the columns of @math are almost linearly dependent. Furthermore, harmonic Ritz vectors may converge irregularly and may even fail to converge. Based on the refined projection principle for large matrix eigenproblems due to the first author, we propose a refined harmonic Lanczos bidiagonalization method that takes the Rayleigh quotients of the harmonic Ritz vectors as approximate singular values and extracts the best approximate singular vectors, called the refined harmonic Ritz approximations, from the given subspaces in the sense of residual minimizations. The refined approximations are shown to converge to the desired singular vectors once the subspaces are sufficiently good and the Rayleigh quotients converge. An implicitly restarted refined harmonic Lanczos bidiagonalization algorithm (IRRHLB) is developed. We study how to select the best possible shifts, and suggest refined harmonic shifts that are theoretically better than the harmonic shifts used within the implicitly restarted harmonic Lanczos bidiagonalization algorithm (IRHLB). We propose a novel procedure that can numerically compute the refined harmonic shifts efficiently and accurately. Numerical experiments are reported that compare IRRHLB with five other algorithms based on the Lanczos bidiagonalization process. It appears that IRRHLB is at least competitive with them and can be considerably more efficient when computing the smallest singular triplets. We consider the behavior of the GMRES method for solving a linear system @math when @math is singular or nearly so, i.e., ill conditioned. The (near) singularity of @math may or may not affect the performance of GMRES, depending on the nature of the system and the initial approximate solution. For singular @math , we give conditions under which the GMRES iterates converge safely to a least-squares solution or to the pseudoinverse solution. These results also apply to any residual minimizing Krylov subspace method that is mathematically equivalent to GMRES. A practical procedure is outlined for efficiently and reliably detecting singularity or ill conditioning when it becomes a threat to the performance of GMRES. We describe a parallel iterative least squares solver named LSRN that is based on random normal projection. LSRN computes the min-length solution to @math , where @math with @math or @math , and where @math may be rank-deficient. Tikhonov regularization may also be included. Since @math is only involved in matrix-matrix and matrix-vector multiplications, it can be a dense or sparse matrix or a linear operator, and LSRN automatically speeds up when @math is sparse or a fast linear operator. The preconditioning phase consists of a random normal projection, which is embarrassingly parallel, and a singular value decomposition of size @math , where @math is moderately larger than 1, e.g., @math . We prove that the preconditioned system is well-conditioned, with a strong concentration result on the extreme singular values, and hence that the number of iterations is fully predictable when we apply LSQR or the Chebyshev semi-iterative method. As we demonstrate, the Chebyshev method is particularly efficient for solving large problems on clusters with high communication cost. Numerical results demonstrate that on a shared-memory machine, LSRN outperforms LAPACK's DGELSD on large dense problems, and MATLAB's backslash (SuiteSparseQR) on sparse problems. Further experiments demonstrate that LSRN scales well on an Amazon Elastic Compute Cloud cluster. Let A be a real m×n matrix with m≧n. It is well known (cf. [4]) that @math (1) where @math The matrix U consists of n orthonormalized eigenvectors associated with the n largest eigenvalues of AA T , and the matrix V consists of the orthonormalized eigenvectors of A T A. The diagonal elements of ∑ are the non-negative square roots of the eigenvalues of A T A; they are called singular values. We shall assume that @math Thus if rank(A)=r, σ r+1 = σ r+2=⋯=σ n = 0. The decomposition (1) is called the singular value decomposition (SVD). We describe the design and implementation of a new algorithm for computing the singular value decomposition (SVD) of a real bidiagonal matrix. This algorithm uses ideas developed by Grosser and Lang that extend Parlett’s and Dhillon’s multiple relatively robust representations (MRRR) algorithm for the tridiagonal symmetric eigenproblem. One key feature of our new implementation is that @math singular triplets can be computed using only @math storage units and floating point operations, where @math is the dimension of the matrix. The algorithm will be made available as routine xBDSCR in the upcoming new release of the LAPACK library. Let @math denote a random symmetric @math matrix whose upper diagonal entries are independent and identically distributed Bernoulli random variables (which take values @math and @math with probability @math each). It is widely conjectured that @math is singular with probability at most @math . On the other hand, the best known upper bound on the singularity probability of @math , due to Vershynin (2011), is @math , for some unspecified small constant @math . This improves on a polynomial singularity bound due to Costello, Tao, and Vu (2005), and a bound of Nguyen (2011) showing that the singularity probability decays faster than any polynomial. In this paper, improving on all previous results, we show that the probability of singularity of @math is at most @math for all sufficiently large @math . The proof utilizes and extends a novel combinatorial approach to discrete random matrix theory, which has been recently introduced by the authors together with Luh and Samotij. An approximate Spielman-Teng theorem for the least singular value @math of a random @math square matrix @math is a statement of the following form: there exist constants @math such that for all @math , @math . The goal of this paper is to develop a simple and novel framework for proving such results for discrete random matrices. As an application, we prove an approximate Spielman-Teng theorem for @math -valued matrices, each of whose rows is an independent vector with exactly @math zero components. This improves on previous work of Nguyen and Vu, and is the first such result in a truly combinatorial' setting. In statistics and machine learning, people are often interested in the eigenvectors (or singular vectors) of certain matrices (e.g. covariance matrices, data matrices, etc). However, those matrices are usually perturbed by noises or statistical errors, either from random sampling or structural patterns. One usually employs Davis-Kahan @math theorem to bound the difference between the eigenvectors of a matrix @math and those of a perturbed matrix @math , in terms of @math norm. In this paper, we prove that when @math is a low-rank and incoherent matrix, the @math norm perturbation bound of singular vectors (or eigenvectors in the symmetric case) is smaller by a factor of @math or @math for left and right vectors, where @math and @math are the matrix dimensions. The power of this new perturbation result is shown in robust covariance estimation, particularly when random variables have heavy tails. There, we propose new robust covariance estimators and establish their asymptotic properties using the newly developed perturbation bound. Our theoretical results are verified through extensive numerical experiments. Super-resolution refers to the process of recovering the locations and amplitudes of a collection of point sources, represented as a discrete measure, given @math of its noisy low-frequency Fourier coefficients. The recovery process is highly sensitive to noise whenever the distance @math between the two closest point sources is less than @math . This paper studies the fundamental difficulty of super-resolution and the performance guarantees of a subspace method called MUSIC in the regime that @math . The most important quantity in our theory is the minimum singular value of the Vandermonde matrix whose nodes are specified by the source locations. Under the assumption that the nodes are closely spaced within several well-separated clumps, we derive a sharp and non-asymptotic lower bound for this quantity. Our estimate is given as a weighted @math sum, where each term only depends on the configuration of each individual clump. This implies that, as the noise increases, the super-resolution capability of MUSIC degrades according to a power law where the exponent depends on the cardinality of the largest clump. Numerical experiments validate our theoretical bounds for the minimum singular value and the resolution limit of MUSIC. When there are @math point sources located on a grid with spacing @math , the fundamental difficulty of super-resolution can be quantitatively characterized by a min-max error, which is the reconstruction error incurred by the best possible algorithm in the worst-case scenario. We show that the min-max error is closely related to the minimum singular value of Vandermonde matrices, and we provide a non-asymptotic and sharp estimate for the min-max error, where the dominant term is @math . Many data analysis applications deal with large matrices and involve approximating the matrix using a small number of “components.” Typically, these components are linear combinations of the rows and columns of the matrix, and are thus difficult to interpret in terms of the original features of the input data. In this paper, we propose and study matrix approximations that are explicitly expressed in terms of a small number of columns and or rows of the data matrix, and thereby more amenable to interpretation in terms of the original data. Our main algorithmic results are two randomized algorithms which take as input an @math matrix @math and a rank parameter @math . In our first algorithm, @math is chosen, and we let @math , where @math is the Moore-Penrose generalized inverse of @math . In our second algorithm @math , @math , @math are chosen, and we let @math . ( @math and @math are matrices that consist of actual columns and rows, respectively, of @math , and @math is a generalized inverse of their intersection.) For each algorithm, we show that with probability at least @math , @math , where @math is the “best” rank- @math approximation provided by truncating the SVD of @math , and where @math is the Frobenius norm of the matrix @math . The number of columns of @math and rows of @math is a low-degree polynomial in @math , @math , and @math . Both the Numerical Linear Algebra community and the Theoretical Computer Science community have studied variants of these matrix decompositions over the last ten years. However, our two algorithms are the first polynomial time algorithms for such low-rank matrix approximations that come with relative-error guarantees; previously, in some cases, it was not even known whether such matrix decompositions exist. Both of our algorithms are simple and they take time of the order needed to approximately compute the top @math singular vectors of @math . The technical crux of our analysis is a novel, intuitive sampling method we introduce in this paper called “subspace sampling.” In subspace sampling, the sampling probabilities depend on the Euclidean norms of the rows of the top singular vectors. This allows us to obtain provable relative-error guarantees by deconvoluting “subspace” information and “size-of- @math ” information in the input matrix. This technique is likely to be useful for other matrix approximation and data analysis problems. Consider a system of diagonal equations satisfying the property that the (fixed) integral coefficient matrix @math contains no singular @math submatrix. A recent paper of the authors [ 3 ] establishes that whenever @math and @math , then the expected asymptotic formula holds for the number @math of integral solutions @math of ( @math ) with @math @math . This paper proves that an "old dog", namely -- the classical Johnson-Lindenstrauss transform, "performs new tricks" -- it gives a novel way of preserving differential privacy. We show that if we take two databases, @math and @math , such that (i) @math is a rank-1 matrix of bounded norm and (ii) all singular values of @math and @math are sufficiently large, then multiplying either @math or @math with a vector of iid normal Gaussians yields two statistically close distributions in the sense of differential privacy. Furthermore, a small, deterministic and alteration of the input is enough to assert that all singular values of @math are large. We apply the Johnson-Lindenstrauss transform to the task of approximating cut-queries: the number of edges crossing a @math -cut in a graph. We show that the JL transform allows us to that preserves edge differential privacy (where two graphs are neighbors if they differ on a single edge) while adding only @math random noise to any given query (w.h.p). Comparing the additive noise of our algorithm to existing algorithms for answering cut-queries in a differentially private manner, we outperform all others on small cuts ( @math ). We also apply our technique to the task of estimating the variance of a given matrix in any given direction. The JL transform allows us to that preserves differential privacy w.r.t bounded changes (each row in the matrix can change by at most a norm-1 vector) while adding random noise of magnitude independent of the size of the matrix (w.h.p). In contrast, existing algorithms introduce an error which depends on the matrix dimensions. This paper makes two contributions to the areas of anti-concentration and non-asymptotic random matrix theory. First, we study the counting problem in inverse Littlewood-Offord theory for general random variables: for random variables @math which are i.i.d. copies of a random variable @math (satisfying some mild hypotheses), how many integer vectors @math in a prescribed box have large @math ? Building on recent work of Ferber, Jain, Luh, and Samotij (who treated the case when @math is a Rademacher random variable), we provide significantly better bounds for this problem than those obtained using the inverse Littlewood-Offord theorems of Tao and Vu, and Nguyen and Vu. Next, we study the non-asymptotic behavior of the least singular value @math of a random @math square matrix @math with i.i.d. entries, which are only assumed to be centered with variance @math . As an application of our counting theorem, and utilizing and developing a recent work by the author, we show that for all @math , @math , where @math is an absolute constant, thereby providing a considerably simpler proof of an approximate version of an essentially optimal result due to Rebrova and Tikhomirov. New methods are proposed for the numerical evaluation of @math or @math , where @math is a function such as @math or @math with singularities in @math and @math is a matrix with eigenvalues on or near @math . The methods are based on combining contour integrals evaluated by the periodic trapezoid rule with conformal maps involving Jacobi elliptic functions. The convergence is geometric, so that the computation of @math is typically reduced to one or two dozen linear system solves, which can be carried out in parallel. Given an @math matrix @math , its Schatten- @math norm, @math , is defined as @math , where @math is the @math -th largest singular value of @math . These norms have been studied in functional analysis in the context of non-commutative @math -spaces, and recently in data stream and linear sketching models of computation. Basic questions on the relations between these norms, such as their embeddability, are still open. Specifically, given a set of matrices @math , suppose we want to construct a linear map @math such that @math for each @math , where @math and @math , and further, @math for a given approximation factor @math and real number @math . Then how large do @math and @math need to be as a function of @math ? We nearly resolve this question for every @math , for the case where @math can be expressed as @math , where @math and @math are arbitrary matrices that are allowed to depend on @math , that is, @math can be implemented by left and right matrix multiplication. Namely, for every @math , we provide nearly matching upper and lower bounds on the size of @math and @math as a function of @math . Importantly, our upper bounds are oblivious , meaning that @math and @math do not depend on the @math , while our lower bounds hold even if @math and @math depend on the @math . As an application of our upper bounds, we answer a recent open question of about space-approximation trade-offs for the Schatten @math -norm, showing in a data stream it is possible to estimate the Schatten- @math norm up to a factor of @math using @math space. A numerically stable and fairly fast scheme is described to compute the unitary matrices U and V which transform a given matrix A into a diagonal form @math , thus exhibiting A’s singular values on @math ’s diagonal. The scheme first transforms A to a bidiagonal matrix J, then diagonalizes J. The scheme described here is complicated but does not suffer from the computational difficulties which occasionally afflict some previously known methods. Some applications are mentioned, in particular the use of the pseudo-inverse @math to solve least squares problems in a way which dampens spurious oscillation and cancellation. The @math -means clustering algorithm is a ubiquitous tool in data mining and machine learning that shows promising performance. However, its high computational cost has hindered its applications in broad domains. Researchers have successfully addressed these obstacles with dimensionality reduction methods. Recently, [1] develop a state-of-the-art random projection (RP) method for faster @math -means clustering. Their method delivers many improvements over other dimensionality reduction methods. For example, compared to the advanced singular value decomposition based feature extraction approach, [1] reduce the running time by a factor of @math for data matrix @math with @math data points and @math features, while losing only a factor of one in approximation accuracy. Unfortunately, they still require @math for matrix multiplication and this cost will be prohibitive for large values of @math and @math . To break this bottleneck, we carefully build a sparse embedded @math -means clustering algorithm which requires @math ( @math denotes the number of non-zeros in @math ) for fast matrix multiplication. Moreover, our proposed algorithm improves on [1]'s results for approximation accuracy by a factor of one. Our empirical studies corroborate our theoretical findings, and demonstrate that our approach is able to significantly accelerate @math -means clustering, while achieving satisfactory clustering performance. We define the hierarchical singular value decomposition (SVD) for tensors of order @math . This hierarchical SVD has properties like the matrix SVD (and collapses to the SVD in @math ), and we prove these. In particular, one can find low rank (almost) best approximations in a hierarchical format ( @math -Tucker) which requires only @math parameters, where @math is the order of the tensor, @math the size of the modes, and @math the (hierarchical) rank. The @math -Tucker format is a specialization of the Tucker format and it contains as a special case all (canonical) rank @math tensors. Based on this new concept of a hierarchical SVD we present algorithms for hierarchical tensor calculations allowing for a rigorous error analysis. The complexity of the truncation (finding lower rank approximations to hierarchical rank @math tensors) is in @math and the attainable accuracy is just 2-3 digits less than machine precision. We consider the problem of robustly recovering a @math -sparse coefficient vector from the Fourier series that it generates, restricted to the interval @math . The difficulty of this problem is linked to the superresolution factor SRF, equal to the ratio of the Rayleigh length (inverse of @math ) by the spacing of the grid supporting the sparse vector. In the presence of additive deterministic noise of norm @math , we show upper and lower bounds on the minimax error rate that both scale like @math , providing a partial answer to a question posed by Donoho in 1992. The scaling arises from comparing the noise level to a restricted isometry constant at sparsity @math , or equivalently from comparing @math to the so-called @math -spark of the Fourier system. The proof involves new bounds on the singular values of restricted Fourier matrices, obtained in part from old techniques in complex analysis. We consider the functions @math defined as the @math th partial derivative of Lebesgue's singular function @math with respect to @math at @math . This sequence includes a multiple of the Takagi function as the case @math . We show that @math is continuous but nowhere differentiable for each @math , and determine the Holder order of @math . From this, we derive that the Hausdorff dimension of the graph of @math is one. Using a formula of Lomnicki and Ulam, we obtain an arithmetic expression for @math using the binary expansion of @math , and use this to find the sets of points where @math and @math take on their absolute maximum and minimum values. We show that these sets are topological Cantor sets. In addition, we characterize the sets of local maximum and minimum points of @math and @math .