37 research outputs found

    An overview of block Gram-Schmidt methods and their stability properties

    Full text link
    Block Gram-Schmidt algorithms serve as essential kernels in many scientific computing applications, but for many commonly used variants, a rigorous treatment of their stability properties remains open. This survey provides a comprehensive categorization of block Gram-Schmidt algorithms, particularly those used in Krylov subspace methods to build orthonormal bases one block vector at a time. All known stability results are assembled, and new results are summarized or conjectured for important communication-reducing variants. Additionally, new block versions of low-synchronization variants are derived, and their efficacy and stability are demonstrated for a wide range of challenging examples. Low-synchronization variants appear remarkably stable for s-step-like matrices built with Newton polynomials, pointing towards a new stable and efficient backbone for Krylov subspace methods. Numerical examples are computed with a versatile MATLAB package hosted at https://github.com/katlund/BlockStab, and scripts for reproducing all results in the paper are provided. Block Gram-Schmidt implementations in popular software packages are discussed, along with a number of open problems. An appendix containing all algorithms type-set in a uniform fashion is provided.Comment: 42 pages, 5 tables, 17 figures, 20 algorithm

    A new look at the Lanczos algorithm for solving symmetric systems of linear equations

    Get PDF
    AbstractSimple versions of the conjugate gradient algorithm and the Lanczos method are discussed, and some merits of the latter are described. A variant of Lanczos is proposed which maintains robust linear independence of the Lanczos vectors by keeping them in secondary storage and occasionally making use of them. The main applications are to problems in which (1) the cost of the matrix-vector product dominates other costs, (2) there is a sequence of right hand sides to be processed, and (3) the eigenvalue distribution of A is not too favorable

    Numerical simulation of a highly underexpanded carbon dioxide jet

    Get PDF
    The underexpanded jets are present in many processes such as rocket propulsion, mass spectrometry, fuel injection, as well as in the process called rapid expansion of supercritical solutions (RESS). In the RESS process a supercritical solution flows through a capillary nozzle until an expansion chamber where the strong changes in the thermodynamic properties of the solvent are used to encapsulate the solute in very fine particles. The research project was focused on the hydrodynamic modeling of an hypersonic carbon dioxide jet produced in the context of the RESS process. The mathematical modeling of the jet was developed using the set of the compressible Navier-Stokes equations along with the generalized Bender equation of state. This set of PDE was solved using an adaptive discontinuous Galerkin discretization for space and the exponential Rosenbrock-Euler method for the time integration. The numerical solver was implemented in C++ using several libraries such as deal.ii and Sacado-Trilinos

    Rank-adaptive structure-preserving reduced basis methods for Hamiltonian systems

    Get PDF
    This work proposes an adaptive structure-preserving model order reduction method for finite-dimensional parametrized Hamiltonian systems modeling non-dissipative phenomena. To overcome the slowly decaying Kolmogorov width typical of transport problems, the full model is approximated on local reduced spaces that are adapted in time using dynamical low-rank approximation techniques. The reduced dynamics is prescribed by approximating the symplectic projection of the Hamiltonian vector field in the tangent space to the local reduced space. This ensures that the canonical symplectic structure of the Hamiltonian dynamics is preserved during the reduction. In addition, accurate approximations with low-rank reduced solutions are obtained by allowing the dimension of the reduced space to change during the time evolution. Whenever the quality of the reduced solution, assessed via an error indicator, is not satisfactory, the reduced basis is augmented in the parameter direction that is worst approximated by the current basis. Extensive numerical tests involving wave interactions, nonlinear transport problems, and the Vlasov equation demonstrate the superior stability properties and considerable runtime speedups of the proposed method as compared to global and traditional reduced basis approaches

    Hierarchical Bayesian Regression with Application in Spatial Modeling and Outlier Detection

    Get PDF
    This dissertation makes two important contributions to the development of Bayesian hierarchical models. The first contribution is focused on spatial modeling. Spatial data observed on a group of areal units is common in scientific applications. The usual hierarchical approach for modeling this kind of dataset is to introduce a spatial random effect with an autoregressive prior. However, the usual Markov chain Monte Carlo scheme for this hierarchical framework requires the spatial effects to be sampled from their full conditional posteriors one-by-one resulting in poor mixing. More importantly, it makes the model computationally inefficient for datasets with large number of units. In this dissertation, we propose a Bayesian approach that uses the spectral structure of the adjacency to construct a low-rank expansion for modeling spatial dependence. We develop a computationally efficient estimation scheme that adaptively selects the functions most important to capture the variation in response. Through simulation studies, we validate the computational efficiency as well as predictive accuracy of our method. Finally, we present an important real-world application of the proposed methodology on a massive plant abundance dataset from Cape Floristic Region in South Africa. The second contribution of this dissertation is a heavy tailed hierarchical regression to detect outliers. We aim to build a linear model that can allow for small as well as large magnitudes of residuals through observation-specific error distribution. t-distribution is specifically suited for that purpose as we can parametrically control its degrees of freedom (df) to tune the heaviness of its tail - large df values represent observations in normal range and small ones represents potential outliers with high error magnitudes. In a hierarchical structure, we can write t-distribution as a scale mixture of a Gaussian distribution so that the standard MCMC algorithm for Gaussian setting can still be used. Post-MCMC, the posterior mean of degrees of freedom for any observation acts as a measure of outlyingness of that observation. We implemented this method on a real dataset consisting of biometric records

    Sparse Supernodal Solver Using Block Low-Rank Compression: design, performance and analysis

    Get PDF
    This paper presents two approaches using a Block Low-Rank (BLR) compressiontechnique to reduce the memory footprint and/or the time-to-solution of the sparse supernodalsolver PaStiX. This flat, non-hierarchical, compression method allows to take advantage of thelow-rank property of the blocks appearing during the factorization of sparse linear systems, whichcome from the discretization of partial differential equations. The first approach, called MinimalMemory, illustrates the maximum memory gain that can be obtained with the BLR compressionmethod, while the second approach, called Just-In-Time, mainly focuses on reducing the com-putational complexity and thus the time-to-solution. Singular Value Decomposition (SVD) andRank-Revealing QR (RRQR), as compression kernels, are both compared in terms of factorizationtime, memory consumption, as well as numerical properties. Experiments on a single node with24 threads and 128 GB of memory are performed to evaluate the potential of both strategies. Ona set of matrices from real-life problems, we demonstrate a memory footprint reduction of up to 4times using the Minimal Memory strategy and a computational time speedup of up to 3.5 timeswith the Just-In-Time strategy. Then, we study the impact of configuration parameters of theBLR solver that allowed us to solve a 3D laplacian of 36 million unknowns a single node, while thefull-rank solver stopped at 8 million due to memory limitation

    Controlling Sensitivity of Gaussian Bayes Predictions based on Eigenvalue Thresholding

    Get PDF
    Gaussian Bayes classifiers are widely used in machine learning for various purposes. Its special characteristic has provided a great capacity for estimating the likelihood and reliability of individual classification decision made, which has been used in many areas such as decision support assessments and risk analysis. However, Gaussian Bayes models tend to perform poorly when processing feature vectors of high dimensionality. This limitation is often resolved using dimension reduction techniques such as Principal Component Analysis. Conventional approaches on reducing dimensionalities usually rely on using a simple threshold based on accuracy measurements or sampling characteristics but rarely consider the sensitivity aspect of the prediction model created. In this paper, we have investigated the influence of eigenvalue selections on Gaussian Bayes classifiers in the context of sensitivity adjustment. Experiments based on real-life data have shown indicative and intriguing results

    Computational Intelligence and Complexity Measures for Chaotic Information Processing

    Get PDF
    This dissertation investigates the application of computational intelligence methods in the analysis of nonlinear chaotic systems in the framework of many known and newly designed complex systems. Parallel comparisons are made between these methods. This provides insight into the difficult challenges facing nonlinear systems characterization and aids in developing a generalized algorithm in computing algorithmic complexity measures, Lyapunov exponents, information dimension and topological entropy. These metrics are implemented to characterize the dynamic patterns of discrete and continuous systems. These metrics make it possible to distinguish order from disorder in these systems. Steps required for computing Lyapunov exponents with a reorthonormalization method and a group theory approach are formalized. Procedures for implementing computational algorithms are designed and numerical results for each system are presented. The advance-time sampling technique is designed to overcome the scarcity of phase space samples and the buffer overflow problem in algorithmic complexity measure estimation in slow dynamics feedback-controlled systems. It is proved analytically and tested numerically that for a quasiperiodic system like a Fibonacci map, complexity grows logarithmically with the evolutionary length of the data block. It is concluded that a normalized algorithmic complexity measure can be used as a system classifier. This quantity turns out to be one for random sequences and a non-zero value less than one for chaotic sequences. For periodic and quasi-periodic responses, as data strings grow their normalized complexity approaches zero, while a faster deceasing rate is observed for periodic responses. Algorithmic complexity analysis is performed on a class of certain rate convolutional encoders. The degree of diffusion in random-like patterns is measured. Simulation evidence indicates that algorithmic complexity associated with a particular class of 1/n-rate code increases with the increase of the encoder constraint length. This occurs in parallel with the increase of error correcting capacity of the decoder. Comparing groups of rate-1/n convolutional encoders, it is observed that as the encoder rate decreases from 1/2 to 1/7, the encoded data sequence manifests smaller algorithmic complexity with a larger free distance value
    corecore