Julia中的快速低秩矩阵逼近_Julia_下载.zip

# LowRankApprox [](https://travis-ci.org/JuliaMatrices/LowRankApprox.jl) [](https://doi.org/10.5281/zenodo.1254147) This Julia package provides fast low-rank approximation algorithms for BLAS/LAPACK-compatible matrices based on some of the latest technology in adaptive randomized matrix sketching. Currently implemented algorithms include: - sketch methods: - random Gaussian - random subset - subsampled random Fourier transform - sparse random Gaussian - partial range finder - partial factorizations: - QR decomposition - interpolative decomposition - singular value decomposition - Hermitian eigendecomposition - CUR decomposition - spectral norm estimation By "partial", we mean essentially that these algorithms are early-terminating, i.e., they are not simply post-truncated versions of their standard counterparts. There is also support for "matrix-free" linear operators described only through their action on vectors. All methods accept a number of options specifying, e.g., the rank, estimated absolute precision, and estimated relative precision of approximation. Our implementation borrows heavily from the perspective espoused by [N. Halko, P.G. Martinsson, J.A. Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM Rev. 53 (2): 217-288, 2011.](http://dx.doi.org/10.1137/090771806), except that we choose the interpolative decomposition (ID) as our basic form of approximation instead of matrix range projection. The reason is that the latter requires expensive matrix-matrix multiplication to contract to the relevant subspace, while the former can sometimes be computed much faster, depending on the accelerated sketching strategy employed. This package has been developed with performance in mind, and early tests have shown large speedups over similar codes written in MATLAB and Python (and even some in Fortran and C). For example, computing an ID of a Hilbert matrix of order 1024 to relative precision ~1e-15 takes: - ~0.02 s using LowRankApprox in Julia - ~0.07 s using SciPy in Python (calling a Fortran backend; see [PyMatrixID](http://klho.github.io/PyMatrixID)) - ~0.3 s in MATLAB This difference can be attributed in part to both algorithmic improvements as well as to some low-level optimizations. ## Contents - [Installation](#installation) - [Getting Started](#getting-started) - [Low-Rank Factorizations](#low-rank-factorizations) - [QR Decomposition](#qr-decomposition) - [Interpolative Decomposition](#interpolative-decomposition-id) - [Singular Value Decomposition](#singular-value-decomposition-svd) - [Hermitian Eigendecomposition](#hermitian-eigendecomposition) - [CUR Decomposition](#cur-decomposition) - [Sketch Methods](#sketch-methods) - [Random Gaussian](#random-gaussian) - [Random Subset](#random-subset) - [Subsampled Random Fourier Transform](#subsampled-random-fourier-transform-srft) - [Sparse Random Gaussian](#sparse-random-gaussian) - [Other Capabilities](#other-capabilities) - [Partial Range](#partial-range) - [Spectral Norm Estimation](#spectral-norm-estimation) - [Core Algorithm](#core-algorithm) - [Options](#options) - [Accuracy Options](#accuracy-options) - [Sketching Options](#sketching-options) - [Other Options](#other-options) - [Computational Complexity](#computational-complexity) ## Installation To install LowRankApprox, simply type: ```julia Pkg.add("LowRankApprox") ``` at the Julia prompt. The package can then be imported as usual with `using` or `import`. ## Getting Started To illustrate the usage of this package, let's consider the computation of a partial QR decomposition of a Hilbert matrix, which is well known to have low rank. First, we load LowRankApprox via: ```julia using LowRankApprox ``` Then we construct a Hilbert matrix with: ```julia n = 1024 A = matrixlib(:hilb, n, n) ``` A partial QR decomposition can then be computed using: ```julia F = pqrfact(A) ``` This returns a `PartialQR` factorization with variables `Q`, `R`, and `p` denoting the unitary, trapezoidal, and permutation factors, respectively, constituting the decomposition. Alternatively, these can be extracted directly with: ```julia Q, R, p = pqr(A) ``` but the factorized form is often more convenient as it also implements various arithmetic operations. For example, the commands: ```julia x = rand(n) y = F *x z = F'*x ``` automatically invoke specialized multiplication routines to rapidly compute `y` and `z` using the low-rank structure of `F`. The rank of the factorization can be retrieved by `F[:k]`, which in this example usually gives 26 or 27. The reason for this variability is that the default interface uses randomized Gaussian sketching for acceleration. Likewise, the actual approximation error is also random but can be efficiently and robustly estimated: ```julia aerr = snormdiff(A, F) rerr = aerr / snorm(A) ``` This computes the absolute and relative errors in the spectral norm using power iteration. Here, the relative error achieved should be on the order of machine epsilon. (You may see a warning about exceeding the iteration limit, but this is harmless in this case.) The default interface requests ~1e-15 estimated relative precision. To request, say, only 1e-12 relative precision, use: ```julia F = pqrfact(A, rtol=1e-12) ``` which returns a factorization of rank ~22. We can also directly control the rank instead with, e.g.: ```julia F = pqrfact(A, rank=20) ``` Using both together as in: ```julia F = pqrfact(A, rank=20, rtol=1e-12) ``` sets two separate termination criteria: one on reaching rank 20 and the other on achieving estimated relative precision 1e-12---with the computation completing upon either of these being fulfilled. Many other options are available as well. All keyword arguments can also be encapsulated into an `LRAOptions` type for convenience. For example, we can equivalently write the above as: ```julia opts = LRAOptions(rank=20, rtol=1e-12) F = pqrfact(A, opts) ``` For further details, see the [Options](#options) section. All aforementioned considerations also apply when the input is a linear operator, i.e., when the matrix is described not by its entries but by its action on vectors. To demonstrate, we can convert `A` to type `LinearOperator` as follows: ```julia L = LinearOperator(A) ``` which inherits and stores methods for applying the matrix and its adjoint. (This command actually recognizes `A` as Hermitian and forms `L` as a `HermitianLinearOperator`.) A partial QR decomposition can then be computed with: ```julia F = pqrfact(L) ``` just as in the previous case. Of course, there is no real benefit to doing so in this particular example; the advantage comes when considering complicated matrix products that can be represented implicitly as a single `LinearOperator`. For instance, `A*A` can be represented as `L*L` without ever forming the resultant matrix explicitly, and we can even encapsulate entire factorizations as linear operators to exploit fast multiplication: ```julia L = LinearOperator(F) ``` Linear operators can be scaled, added, and composed together using the usual syntax. All methods in LowRankApprox transparently support both matrices and linear operators. ## Low-Rank Factorizations We now detail the various low-rank approximations implemented, which all nominally return compact `Factorization` types storing the matrix factors in structured form. All such factorizations provide optimized multiplication routines. Furthermore, the rank of any factorization `F` can be queried with `F[:k]` and the matrix approximant defined by `F` can be reconstructed as `Matrix(F)`. For concreteness of exposition, assume in the following that `A` has size `m` by `n` with factorization rank `F[:k] = k`. Note that certain matrix id