Subspace Approximation Research Articles

Approximating the large sparse matrix exponential using incomplete orthogonalization and Krylov subspaces of variable dimension

SummaryKrylov subspace approximations to the matrix exponential are popularly used with full orthogonalization instead of incomplete orthogonalization, even though the latter strategy is known to reduce the cost by truncating the recurrences of the modified Gram–Schmidt process. This study combines such a strategy with an adaptive step‐by‐step integration scheme that allows both the stepsize and the dimension of the Krylov subspace to vary. A convergence analysis is done. Numerical results on test problems drawn from systems biology and computer systems show a significant speedup over the standard implementation with full orthogonalization and fixed dimension.

An investigation of GPU-based stiff chemical kinetics integration methods

A fifth-order implicit Runge–Kutta method and two fourth-order exponential integration methods equipped with Krylov subspace approximations were implemented for the GPU and paired with the analytical chemical kinetic Jacobian software pyJac. The performance of each algorithm was evaluated by integrating thermochemical state data sampled from stochastic partially stirred reactor simulations and compared with the commonly used CPU-based implicit integrator CVODE. We estimated that the implicit Runge–Kutta method running on a single Tesla C2075 GPU is equivalent to CVODE running on 12–38 Intel Xeon E5-4640 v2 CPU cores for integration of a single global integration time step of 10−6s with hydrogen and methane kinetic models. In the stiffest case studied—the methane model with a global integration time step of 10−4s—thread divergence and higher memory traffic significantly decreased GPU performance to the equivalent of CVODE running on approximately three CPU cores. The exponential integration algorithms performed more slowly than the implicit integrators on both the CPU and GPU. Thread divergence and memory traffic were identified as the main limiters of GPU integrator performance, and techniques to mitigate these issues were discussed. Use of a finite-difference Jacobian on the GPU—in place of the analytical Jacobian provided by pyJac—greatly decreased integrator performance due to thread divergence, resulting in maximum slowdowns of 7.11−240.96×; in comparison, the corresponding slowdowns on the CPU were just 1.39−2.61×, underscoring the importance of use of an analytical Jacobian for efficient GPU integration. Finally, future research directions for working towards enabling realistic chemistry in reactive-flow simulations via GPU/SIMT accelerated stiff chemical kinetics integration were identified.

Rational approximation of affine coordinate subspaces of Euclidean space

We show that affine coordinate subspaces of dimension at least two in Euclidean space are of Khintchine type for divergence. For affine coordinate subspaces of dimension one, we prove a result which depends on the dual Diophantine type of the basepoint of the subspace. These results provide evidence for the conjecture that all affine subspaces of Euclidean space are of Khintchine type for divergence.

Rectifiable-Reifenberg and the regularity of stationary and minimizing harmonic maps

In this paper we study the regularity of stationary and minimizing harmonic maps $f:B_2(p)\subseteq M\to N$ between Riemannian manifolds. If $S^k(f)\equiv\{x\in M: \text{ no tangent map at $x$ is }k+1\text{-symmetric}\}$ is $k^{th}$-stratum of the singular set of $f$, then it is well known that $\dim S^k\leq k$, however little else about the structure of $S^k(f)$ is understood in any generality. Our first result is for a general stationary harmonic map, where we prove that $S^k(f)$ is $k$-rectifiable. In the case of minimizing harmonic maps we go further, and prove that the singular set $S(f)$, which is well known to satisfy $\dim S(f)\leq n-3$, is in fact $n-3$-rectifiable with uniformly {\it finite} $n-3$-measure. An effective version of this allows us to prove that $|\nabla f|$ has estimates in $L^3_{weak}$, an estimate which is sharp as $|\nabla f|$ may not live in $L^3$. The above results are in fact just applications of a new class of estimates we prove on the {\it quantitative} stratifications $S^k_{\epsilon,r}(f)$ and $S^k_{\epsilon}(f)\equiv S^k_{\epsilon,0}(f)$. Roughly, $S^k_{\epsilon}\subseteq M$ is the collection of points $x\in S^k_\epsilon$ for which no ball $B_r(x)$ is $\epsilon$-close to being $k+1$-symmetric. We show that $S^k_\epsilon$ is $k$-rectifiable and satisfies the Minkowski estimate $Vol(B_r\,S_\epsilon^k)\leq C r^{n-k}$. The proofs require a new $L^2$-subspace approximation theorem for stationary harmonic maps, as well as new $W^{1,p}$-Reifenberg and rectifiable-Reifenberg type theorems. These results are generalizations of the classical Reifenberg, and give checkable criteria to determine when a set is $k$-rectifiable with uniform measure estimates. The new Reifenberg type theorems may be of some independent interest.

An Efficient Reduced Basis Solver for Stochastic Galerkin Matrix Equations

Stochastic Galerkin finite element approximation of PDEs with random inputs leads to linear systems of equations with coefficient matrices that have a characteristic Kronecker product structure. By reformulating the systems as multiterm linear matrix equations, we develop an efficient solution algorithm which generalizes ideas from rational Krylov subspace approximation. Our working assumptions are that the number of random variables characterizing the random inputs is modest, in the order of a few tens, and that the dependence on these variables is linear, so that it is sufficient to seek only a reduction in the complexity associated with the spatial component of the approximation space. The new approach determines a low-rank approximation to the solution matrix by performing a projection onto a low-dimensional space and provides an efficient solution strategy whose convergence rate is independent of the spatial approximation. Moreover, it requires far less memory than the standard preconditioned conjuga...

A Hybrid Riemann Solver for Large Hyperbolic Systems of Conservation Laws

This paper presents a Riemann solver for nonlinear hyperbolic systems of conservation laws based on a Krylov subspace approximation of the upwinding dissipation vector. In the general case, the solver does not require any detailed information of the eigensystem, except an estimate of the global maximal propagation speed. It uses successive flux function evaluations to obtain a numerical flux which is almost equivalent to that of a Godunov scheme with complete upwinding. The new Krylov--Riemann solver is particularly efficient when used for large systems with many nonlinear equations such that typically no explicit expression for the eigensystem is available. Also, no numerical procedures are necessary to compute the eigensystem. Numerical examples demonstrate the excellent performance of the solver with respect to other solvers.

Hybrid and Iteratively Reweighted Regularization by Unbiased Predictive Risk and Weighted GCV for Projected Systems

Tikhonov regularization for projected solutions of large-scale ill-posed problems is considered. The Golub--Kahan iterative bidiagonalization is used to project the problem onto a subspace and regularization then applied to find a subspace approximation to the full problem. Determination of the regularization, parameter for the projected problem by unbiased predictive risk estimation, generalized cross validation, and discrepancy principle techniques is investigated. It is shown that the regularized parameter obtained by the unbiased predictive risk estimator can provide a good estimate which can be used for a full problem that is moderately to severely ill-posed. A similar analysis provides the weight parameter for the weighted generalized cross validation such that the approach is also useful in these cases, and also explains why the generalized cross validation without weighting is not always useful. All results are independent of whether systems are over-or underdetermined. Numerical simulations for standard one-dimensional test problems and two-dimensional data, for both image restoration and tomographic image reconstruction, support the analysis and validate the techniques. The size of the projected problem is found using an extension of a noise revealing function for the projected problem [I. Hnětynková, M. Plešinger, and Z. Strakoš, BIT Numer. Math., 49 (2009), pp. 669--696]. Furthermore, an iteratively reweighted regularization approach for edge preserving regularization is extended for projected systems, providing stabilization of the solutions of the projected systems and reducing dependence on the determination of the size of the projected subspace.

A Conversation on Data Mining Strategies in LC-MS Untargeted Metabolomics: Pre-Processing and Pre-Treatment Steps.

Untargeted metabolomic studies generate information-rich, high-dimensional, and complex datasets that remain challenging to handle and fully exploit. Despite the remarkable progress in the development of tools and algorithms, the “exhaustive” extraction of information from these metabolomic datasets is still a non-trivial undertaking. A conversation on data mining strategies for a maximal information extraction from metabolomic data is needed. Using a liquid chromatography-mass spectrometry (LC-MS)-based untargeted metabolomic dataset, this study explored the influence of collection parameters in the data pre-processing step, scaling and data transformation on the statistical models generated, and feature selection, thereafter. Data obtained in positive mode generated from a LC-MS-based untargeted metabolomic study (sorghum plants responding dynamically to infection by a fungal pathogen) were used. Raw data were pre-processed with MarkerLynxTM software (Waters Corporation, Manchester, UK). Here, two parameters were varied: the intensity threshold (50–100 counts) and the mass tolerance (0.005–0.01 Da). After the pre-processing, the datasets were imported into SIMCA (Umetrics, Umea, Sweden) for more data cleaning and statistical modeling. In addition, different scaling (unit variance, Pareto, etc.) and data transformation (log and power) methods were explored. The results showed that the pre-processing parameters (or algorithms) influence the output dataset with regard to the number of defined features. Furthermore, the study demonstrates that the pre-treatment of data prior to statistical modeling affects the subspace approximation outcome: e.g., the amount of variation in X-data that the model can explain and predict. The pre-processing and pre-treatment steps subsequently influence the number of statistically significant extracted/selected features (variables). Thus, as informed by the results, to maximize the value of untargeted metabolomic data, understanding of the data structures and exploration of different algorithms and methods (at different steps of the data analysis pipeline) might be the best trade-off, currently, and possibly an epistemological imperative.

Acceleration of contour integration techniques by rational Krylov subspace methods

We suggest a rational Krylov subspace approximation for products of matrix functions and a vector appearing in exponential integrators. We consider matrices with a field-of-values in a sector lying in the left complex half-plane. The choice of the poles for our method is suggested by a fixed rational approximation based on contour integration along a hyperbola around the sector. Compared to the fixed approximation, our rational Krylov subspace method exhibits an accelerated and more stable convergence of order O(e−Cn).

Sub-space approximations for MDO problems with disparate disciplinary variable dependence

An approach to solving multidisciplinary design optimisation problems using approximations built in sub-spaces of the design variable space is proposed. Each approximation is built in the sub-space significant to the corresponding discipline while the optimisation problem is solved in the full design variable space. Since the approximations are built in a space of reduced dimensionality, the computational budget associated with building them can be reduced without compromising their quality. The method requires the designer to make assumptions on which design variables are significant to each discipline. If such assumptions are deficient, the resulting approximations suffer from errors that are not possible to reduce by additional sampling. Therefore a recovery mechanism is proposed that updates the values of the insignificant variables at the end of each iteration to align with the current best point. The method is implemented within a trust region based optimisation framework and demonstrated on a multidisciplinary optimisation of a thin-walled beam section subject to stiffness and impact requirements.

Krylov Integration Factor Method on Sparse Grids for High Spatial Dimension Convection–Diffusion Equations

Krylov implicit integration factor (IIF) methods were developed in Chen and Zhang (J Comput Phys 230:4336---4352, 2011) for solving stiff reaction---diffusion equations on high dimensional unstructured meshes. The methods were further extended to solve stiff advection---diffusion---reaction equations in Jiang and Zhang (J Comput Phys 253:368---388, 2013). Recently we studied the computational power of Krylov subspace approximations on dealing with high dimensional problems. It was shown that the Krylov integration factor methods have linear computational complexity and are especially efficient for high dimensional convection---diffusion problems with anisotropic diffusions. In this paper, we combine the Krylov integration factor methods with sparse grid combination techniques and solve high spatial dimension convection---diffusion equations such as Fokker---Planck equations on sparse grids. Numerical examples are presented to show that significant computational times are saved by applying the Krylov integration factor methods on sparse grids.

A versatile sparse representation based post-processing method for improving image super-resolution

The objective of this work is single image super-resolution (SR), in which the input is specified by a low-resolution image and a consistent higher-resolution image should be returned. We propose a novel post-processing procedure named iterative fine-tuning and approximation (IFA) for mainstream SR methods. Internal image statistics are complemented by iteratively fine-tuning and performing linear subspace approximation on the outputs of existing external SR methods, helping to better reconstruct missing details and reduce unwanted artifacts. The primary concept of our method is that it first explores and enhances internal image information by grouping similar image patches and then finds their sparse or low-rank representations by iteratively learning the bases or primary components, thereby enhancing the primary structures and some details of the image. We evaluate the proposed IFA procedure over two standard benchmark datasets and demonstrate that IFA can yield substantial improvements for most existing methods via tweaking their outputs, achieving state-of-the-art performance.

Memory Efficient PCA Methods for Large Group ICA.

Principal component analysis (PCA) is widely used for data reduction in group independent component analysis (ICA) of fMRI data. Commonly, group-level PCA of temporally concatenated datasets is computed prior to ICA of the group principal components. This work focuses on reducing very high dimensional temporally concatenated datasets into its group PCA space. Existing randomized PCA methods can determine the PCA subspace with minimal memory requirements and, thus, are ideal for solving large PCA problems. Since the number of dataloads is not typically optimized, we extend one of these methods to compute PCA of very large datasets with a minimal number of dataloads. This method is coined multi power iteration (MPOWIT). The key idea behind MPOWIT is to estimate a subspace larger than the desired one, while checking for convergence of only the smaller subset of interest. The number of iterations is reduced considerably (as well as the number of dataloads), accelerating convergence without loss of accuracy. More importantly, in the proposed implementation of MPOWIT, the memory required for successful recovery of the group principal components becomes independent of the number of subjects analyzed. Highly efficient subsampled eigenvalue decomposition techniques are also introduced, furnishing excellent PCA subspace approximations that can be used for intelligent initialization of randomized methods such as MPOWIT. Together, these developments enable efficient estimation of accurate principal components, as we illustrate by solving a 1600-subject group-level PCA of fMRI with standard acquisition parameters, on a regular desktop computer with only 4 GB RAM, in just a few hours. MPOWIT is also highly scalable and could realistically solve group-level PCA of fMRI on thousands of subjects, or more, using standard hardware, limited only by time, not memory. Also, the MPOWIT algorithm is highly parallelizable, which would enable fast, distributed implementations ideal for big data analysis. Implications to other methods such as expectation maximization PCA (EM PCA) are also presented. Based on our results, general recommendations for efficient application of PCA methods are given according to problem size and available computational resources. MPOWIT and all other methods discussed here are implemented and readily available in the open source GIFT software.

Krylov single-step implicit integration factor WENO methods for advection–diffusion–reaction equations

Implicit integration factor (IIF) methods were developed in the literature for solving time-dependent stiff partial differential equations (PDEs). Recently, IIF methods were combined with weighted essentially non-oscillatory (WENO) schemes in Jiang and Zhang (2013) [19] to efficiently solve stiff nonlinear advection–diffusion–reaction equations. The methods can be designed for arbitrary order of accuracy. The stiffness of the system is resolved well and the methods are stable by using time step sizes which are just determined by the non-stiff hyperbolic part of the system. To efficiently calculate large matrix exponentials, Krylov subspace approximation is directly applied to the implicit integration factor (IIF) methods. So far, the IIF methods developed in the literature are multistep methods. In this paper, we develop Krylov single-step IIF-WENO methods for solving stiff advection–diffusion–reaction equations. The methods are designed carefully to avoid generating positive exponentials in the matrix exponentials, which is necessary for the stability of the schemes. We analyze the stability and truncation errors of the single-step IIF schemes. Numerical examples of both scalar equations and systems are shown to demonstrate the accuracy, efficiency and robustness of the new methods.

A Krylov semi-implicit discontinuous Galerkin method for the computation of ground and excited states in Bose–Einstein condensates

Numerical computations of ground and excited states of Bose–Einstein condensates (BECs) require high spatial resolution due to the presence of the multiscale structures and boundary layers and interior layers in the solution. In this paper, a new discontinuous Galerkin (DG) method is presented for the computation of ground and excited states of BECs. For the spatial discretization, the direct discontinuous Galerkin (DDG) method is applied by using the normalized gradient flow. For the time discretization, we integrate the ordinary differential equations (ODEs) which is obtained by linearizing the interaction term. To evaluate the matrix exponential operator efficiently, we apply the Krylov subspace approximations to the matrix exponential operator. Numerical examples with different potentials are reported to demonstrate the validity and effectiveness of the semi-implicit DG method.

Error bounds and estimates for Krylov subspace approximations of Stieltjes matrix functions

When using the Lanczos method to approximate $$f(A){\varvec{b}}$$ , the action of a matrix function on a vector, there is, in contrast to the solution of linear systems, no straightforward way to measure or estimate the error of the current iterate. Therefore, to be able to decide whether the desired accuracy has been reached, several different estimates and bounds for the error have been suggested, all of them specific to certain classes of functions. In this paper, we add to these results by developing a technique to compute error bounds for Stieltjes functions, using a recently suggested integral representation of the error, and we show how these bounds can be computed essentially for free, i.e., with cost independent of the iteration number and the dimension of the matrix A.

Robust Environmental Sound Recognition With Fast Noise Suppression for Home Automation

This paper proposes a robust environmental sound recognition system using a fast noise suppression approach for home automation applications. The system comprises a fast subspace-based noise suppression module and a sound classification module. For the noise suppression module, we propose a noise suppression method that applies fast subspace approximations in the wavelet domain. We show that this method offers a lower computational cost than conventional methods. In the sound classification module, we use a feature extraction method that is also based on the wavelet subspace, derived from seventeen critical bands in a signal's wavelet packet transform. Furthermore, we create a multiclass support vector machine by employing probability product kernels. The experimental results for ten classes of various environmental sounds show that the proposed system offers robust performance in environmental sound recognition tasks.

SubspaceEM: A fast maximum-a-posteriori algorithm for cryo-EM single particle reconstruction.

Single particle reconstruction methods based on the maximum-likelihood principle and the expectation-maximization (E-M) algorithm are popular because of their ability to produce high resolution structures. However, these algorithms are computationally very expensive, requiring a network of computational servers. To overcome this computational bottleneck, we propose a new mathematical framework for accelerating maximum-likelihood reconstructions. The speedup is by orders of magnitude and the proposed algorithm produces similar quality reconstructions compared to the standard maximum-likelihood formulation. Our approach uses subspace approximations of the cryo-electron microscopy (cryo-EM) data and projection images, greatly reducing the number of image transformations and comparisons that are computed. Experiments using simulated and actual cryo-EM data show that speedup in overall execution time compared to traditional maximum-likelihood reconstruction reaches factors of over 300.

Fast Subspace Approximation Via Greedy Least-Squares

In this note, we develop fast and deterministic dimensionality reduction techniques for a family of subspace approximation problems. We then utilize these dimensionality reduction techniques to help rapidly and accurately approximate the \(n\)-widths of point sets. Let \(P \subset \mathbb {R}^N\) be a given set of \(M\) points. The techniques developed herein find an \(O(n \log M)\)-dimensional subspace that is guaranteed to always contain a near-best fit \(n\)-dimensional hyperplane \(\mathcal {H}\) for \(P\) with respect to the cumulative projection error \(\left( \sum _{\mathbf{x} \in P} \Vert \mathbf{x} - \Pi _{\mathcal {H}} \mathbf{x} \Vert ^p_2 \right) ^{1/p}\), for any chosen \(p > 2\). The deterministic algorithm runs in \(\tilde{O} \left( MN^2 \right) \)-time, and can be randomized to run in only \(\tilde{O} \left( MNn \right) \)-time while maintaining its error guarantees with high probability. In the important \(p = \infty \) case, the dimensionality reduction techniques are then combined with efficient algorithms for computing the John ellipsoid of a data set in order to produce an \(n\)-dimensional subspace whose maximum \(\ell _2\)-distance to any point in the convex hull of \(P\) is minimized. The resulting algorithm remains in \(\tilde{O} \left( MNn \right) \)-time.

SubPatCNV: approximate subspace pattern mining for mapping copy-number variations.

BackgroundMany DNA copy-number variations (CNVs) are known to lead to phenotypic variations and pathogenesis. While CNVs are often only common in a small number of samples in the studied population or patient cohort, previous work has not focused on customized identification of CNV regions that only exhibit in subsets of samples with advanced data mining techniques to reliably answer questions such as “Which are all the chromosomal fragments showing nearly identical deletions or insertions in more than 30% of the individuals?”.ResultsWe introduce a tool for mining CNV subspace patterns, namely SubPatCNV, which is capable of identifying all aberrant CNV regions specific to arbitrary sample subsets larger than a support threshold. By design, SubPatCNV is the implementation of a variation of approximate association pattern mining algorithm under a spatial constraint on the positional CNV probe features. In benchmark test, SubPatCNV was applied to identify population specific germline CNVs from four populations of HapMap samples. In experiments on the TCGA ovarian cancer dataset, SubPatCNV discovered many large aberrant CNV events in patient subgroups, and reported regions enriched with cancer relevant genes. In both HapMap data and TCGA data, it was observed that SubPatCNV employs approximate pattern mining to more effectively identify CNV subspace patterns that are consistent within a subgroup from high-density array data.ConclusionsSubPatCNV available through http://sourceforge.net/projects/subpatcnv/is a unique scalable open-source software tool that provides the flexibility of identifying CNV regions specific to sample subgroups of different sizes from high-density CNV array data.Electronic supplementary materialThe online version of this article (doi:10.1186/s12859-014-0426-7) contains supplementary material, which is available to authorized users.