Computer Science - Numerical Analysis Publications (50)


Computer Science - Numerical Analysis Publications

A new simple convergence acceleration method is proposed for a certain wide range class of convergent alternating series. The method has some common features with Smith's and Ford's modification of Levin's and Weniger's sequence transformations, but it is computationally less expensive. The similarities and differences between all three methods are analyzed and some common theoretical results are given. Read More

In an effort to increase the versatility of finite element codes, we explore the possibility of automatically creating the Jacobian matrix necessary for the gradient-based solution of nonlinear systems of equations. Particularly, we aim to assess the feasibility of employing the automatic differentiation tool TAPENADE for this purpose on a large Fortran codebase that is the result of many years of continuous development. As a starting point we will describe the special structure of finite element codes and the implications that this code design carries for an efficient calculation of the Jacobian matrix. Read More

Boolean matrix factorisation (BooMF) infers interpretable decompositions of a binary data matrix into a pair of low-rank, binary matrices: One containing meaningful patterns, the other quantifying how the observations can be expressed as a combination of these patterns. We introduce the OrMachine, a probabilistic generative model for BooMF and derive a Metropolised Gibbs sampler that facilitates very efficient parallel posterior inference. Our method outperforms all currently existing approaches for Boolean Matrix factorization and completion, as we show on simulated and real world data. Read More

Block Coordinate Update (BCU) methods enjoy low per-update computational complexity because every time only one or a few block variables would need to be updated among possibly a large number of blocks. They are also easily parallelized and thus have been particularly popular for solving problems involving large-scale dataset and/or variables. In this paper, we propose a primal-dual BCU method for solving linearly constrained convex program in multi-block variables. Read More

We address the statistical and optimization impacts of using classical or Hessian sketch to approximately solve the Matrix Ridge Regression (MRR) problem. Prior research has considered the effects of classical sketch on least squares regression (LSR), a strictly simpler problem. We establish that classical sketch has a similar effect upon the optimization properties of MRR as it does on those of LSR -- namely, it recovers nearly optimal solutions. Read More

For a wide variety of problems, creating detailed continuous models of (continuous) physical systems is, at the very least, impractical. Hybrid models can abstract away short transient behaviour (thus introducing discontinuities) in order to simplify the study of such systems. For example, when modelling a bouncing ball, the bounce can be abstracted as a discontinuous change of the velocity, instead of resorting to the physics of the ball (de-)compression to keep the velocity signal continuous. Read More

With the huge influx of various data nowadays, extracting knowledge from them has become an interesting but tedious task among data scientists, particularly when the data come in heterogeneous form and have missing information. Many data completion techniques had been introduced, especially in the advent of kernel methods. However, among the many data completion techniques available in the literature, studies about mutually completing several incomplete kernel matrices have not been given much attention yet. Read More

The emergent field of probabilistic numerics has thus far lacked rigorous statistical principals. This paper establishes Bayesian probabilistic numerical methods as those which can be cast as solutions to certain Bayesian inverse problems, albeit problems that are non-standard. This allows us to establish general conditions under which Bayesian probabilistic numerical methods are well-defined, encompassing both non-linear and non-Gaussian models. Read More

This paper presents a randomized algorithm for computing the near-optimal low-rank dynamic mode decomposition (DMD). Randomized algorithms are emerging techniques to compute low-rank matrix approximations. They are able to ease the computational challenges arising in the area of big data. Read More

Three complementary methods have been implemented in the code Denovo that accelerate neutral particle transport calculations with methods that use leadership-class computers fully and effectively: a multigroup block (MG) Krylov solver, a Rayleigh quotient iteration (RQI) eigenvalue solver, and a multigrid in energy preconditioner. The multigroup Krylov solver converges more quickly than Gauss Seidel and enables energy decomposition such that Denovo can scale to hundreds of thousands of cores. The new multigrid in energy preconditioner reduces iteration count for many problem types and takes advantage of the new energy decomposition such that it can scale efficiently. Read More

Lower bounds on the smallest eigenvalue of a symmetric positive definite matrices $A\in\mathbb{R}^{m\times m}$ play an important role in condition number estimation and in iterative methods for singular value computation. In particular, the bounds based on ${\rm Tr}(A^{-1})$ and ${\rm Tr}(A^{-2})$ attract attention recently because they can be computed in $O(m)$ work when $A$ is tridiagonal. In this paper, we focus on these bounds and investigate their properties in detail. Read More

This paper proposes an efficient method for computing selected generalized eigenpairs of a sparse Hermitian definite matrix pencil $(A,B)$. Based on Zolotarev's best rational function approximations of the signum function and conformal maps, we construct the best rational function approximation of a rectangular function supported on an arbitrary interval. This new best rational function approximation is applied to construct spectrum filters of $(A,B)$. Read More

This document contains the notes of a lecture I gave at the "Journ\'ees Nationales du Calcul Formel" (JNCF) on January 2017. The aim of the lecture was to discuss low-level algorithmics for p-adic numbers. It is divided into two main parts: first, we present various implementations of p-adic numbers and compare them and second, we introduce a general framework for studying precision issues and apply it in several concrete situations. Read More

The CANDECOMP/PARAFAC (CP) decomposition is a leading method for the analysis of multiway data. The standard alternating least squares algorithm for the CP decomposition (CP-ALS) involves a series of highly overdetermined linear least squares problems. We extend randomized least squares methods to tensors and show the workload of CP-ALS can be drastically reduced without a sacrifice in quality. Read More

Affiliations: 1Virginia Tech, 2Virginia Tech, 3Virginia Tech, 4University of California, Merced

Exponential integrators are special time discretization methods where the traditional linear system solves used by implicit schemes are replaced with computing the action of matrix exponential-like functions on a vector. A very general formulation of exponential integrators is offered by the Exponential Propagation Iterative methods of Runge-Kutta type (EPIRK) family of schemes. The use of Jacobian approximations is an important strategy to drastically reduce the overall computational costs of implicit schemes while maintaining the quality of their solutions. Read More

In this paper a new hp-adaptive strategy for elliptic problems based on refinement history is proposed, which chooses h-, p- or hp-refinement on individual elements according to a posteriori error estimate, as well as smoothness estimate of the solution obtained by comparing the actual and expected error reduction rate. Numerical experiments show that exponential convergence can be achieved with this strategy. Read More

Implicit schemes have been extensively used in building physics to compute the solution of moisture diffusion problems in porous materials for improving stability conditions. Nevertheless, these schemes require important sub-iterations when treating non-linear problems. To overcome this disadvantage, this paper explores the use of improved explicit schemes, such as Dufort-Frankel, Crank-Nicolson and hyperbolisation approaches. Read More

High order cumulants carry information about statistics of non--normally distributed multivariate data. Such cumulants are utilised in extreme events analysis, small target detection or outliers detection. In this work we present a new algorithm,for updating high order cumulant tensors of random multivariate data, if new package of data is recorded. Read More

In this paper we introduce a novel algorithm of calculating arbitrary order cumulants of multidimensional data. Since the n th order cumulant can be presented in the form of an n-dimensional tensor, the algorithm is presented using the tensor network notation. The presented algorithm exploits the super--symmetry of cumulant and moment tensors. Read More

This paper develops meshless methods for probabilistically describing discretisation error in the numerical solution of partial differential equations. This construction enables the solution of Bayesian inverse problems while accounting for the impact of the discretisation of the forward problem. In particular, this drives statistical inferences to be more conservative in the presence of significant solver error. Read More

On modern large-scale parallel computers, the performance of Krylov subspace iterative methods is limited by global synchronization. This has inspired the development of $s$-step Krylov subspace method variants, in which iterations are computed in blocks of $s$, which can reduce the number of global synchronizations per iteration by a factor of $O(s)$. Although the $s$-step variants are mathematically equivalent to their classical counterparts, they can behave quite differently in finite precision depending on the parameter $s$. Read More

This paper introduces multivariate input-output models to predict the errors and bases dimensions of local parametric Proper Orthogonal Decomposition reduced-order models. We refer to these multivariate mappings as the MP-LROM models. We employ Gaussian Processes and Artificial Neural Networks to construct approximations of these multivariate mappings. Read More

A fully (pseudo-)spectral solver for direct numerical simulations of large-scale turbulent channel flows is described. The solver utilizes the Chebyshev base functions suggested by J. Shen [SIAM J. Read More

In this work, we propose two-level space-time domain decomposition preconditioners for parabolic problems discretized using finite elements. They are motivated as an extension to space-time of balancing domain decomposition by constraints preconditioners. The key ingredients to be defined are the sub-assembled space and operator, the coarse degrees of freedom (DOFs) in which we want to enforce continuity among subdomains at the preconditioner level, and the transfer operator from the sub-assembled to the original finite element space. Read More

This work proposes a machine-learning-based framework for estimating the error introduced by surrogate models of parameterized dynamical systems. The framework applies high-dimensional regression techniques (e.g. Read More

Kernel matrices appear in machine learning and non-parametric statistics. Given $N$ points in $d$ dimensions and a kernel function that requires $\mathcal{O}(d)$ work to evaluate, we present an $\mathcal{O}(dN\log N)$-work algorithm for the approximate factorization of a regularized kernel matrix, a common computational bottleneck in the training phase of a learning task. With this factorization, solving a linear system with a kernel matrix can be done with $\mathcal{O}(N\log N)$ work. Read More

Differential (Ore) type polynomials with "approximate" polynomial coefficients are introduced. These provide an effective notion of approximate differential operators, with a strong algebraic structure. We introduce the approximate Greatest Common Right Divisor Problem (GCRD) of differential polynomials, as a non-commutative generalization of the well-studied approximate GCD problem. Read More

Design of filters for graph signal processing benefits from knowledge of the spectral decomposition of matrices that encode graphs, such as the adjacency matrix and the Laplacian matrix, used to define the shift operator. For shift matrices with real eigenvalues, which arise for symmetric graphs, the empirical spectral distribution captures the eigenvalue locations. Under realistic circumstances, stochastic influences often affect the network structure and, consequently, the shift matrix empirical spectral distribution. Read More

The paper derives and analyses the (semi-)discrete dispersion relation of the Parareal parallel-in-time integration method. It investigates Parareal's wave propagation characteristics with the aim to better understand what causes the well documented stability problems for hyperbolic equations. The analysis shows that the instability is caused by convergence of the amplification factor to the exact value from above for medium to high wave numbers. Read More

The paper covers a formulation of the inverse quadratic programming problem in terms of unconstrained optimization where it is required to find the unknown parameters (the matrix of the quadratic form and the vector of the quasi-linear part of the quadratic form) provided that approximate estimates of the optimal solution of the direct problem and those of the target function to be minimized in the form of pairs of values lying in the corresponding neighborhoods are only known. The formulation of the inverse problem and its solution are based on the least squares method. In the explicit form the inverse problem solution has been derived in the form a system of linear equations. Read More

When a measurement falls outside the quantization or measurable range, it becomes saturated and cannot be used in classical reconstruction methods. For example, in C-arm angiography systems, which provide projection radiography, fluoroscopy, digital subtraction angiography, and are widely used for medical diagnoses and interventions, the limited dynamic range of C-arm flat detectors leads to overexposure in some projections during an acquisition, such as imaging relatively thin body parts (e.g. Read More

We analyzed the performance of a biologically inspired algorithm called the Corrected Projections Algorithm (CPA) when a sparseness constraint is required to unambiguously reconstruct an observed signal using atoms from an overcomplete dictionary. By changing the geometry of the estimation problem, CPA gives an analytical expression for a binary variable that indicates the presence or absence of a dictionary atom using an L2 regularizer. The regularized solution can be implemented using an efficient real-time Kalman-filter type of algorithm. Read More

This thesis examines a modern concept for machine numbers based on interval arithmetic called 'Unums' and compares it to IEEE 754 floating-point arithmetic, evaluating possible uses of this format where floating-point numbers are inadequate. In the course of this examination, this thesis builds theoretical foundations for IEEE 754 floating-point numbers, interval arithmetic based on the projectively extended real numbers and Unums. Read More

This report describes the computation of gradients by algorithmic differentiation for statistically optimum beamforming operations. Especially the derivation of complex-valued functions is a key component of this approach. Therefore the real-valued algorithmic differentiation is extended via the complex-valued chain rule. Read More

Non-negative matrix factorization (NMF) is a prob- lem with many applications, ranging from facial recognition to document clustering. However, due to the variety of algorithms that solve NMF, the randomness involved in these algorithms, and the somewhat subjective nature of the problem, there is no clear "correct answer" to any particular NMF problem, and as a result, it can be hard to test new algorithms. This paper suggests some test cases for NMF algorithms derived from matrices with enumerable exact non-negative factorizations and perturbations of these matrices. Read More

This paper introduces optimally-blended quadrature rules for isogeometric analysis and analyzes the numerical dispersion of the resulting discretizations. To quantify the approximation errors when we modify the inner products, we generalize the Pythagorean eigenvalue theorem of Strang and Fix. The proposed blended quadrature rules have advantages over alternative integration rules for isogeometric analysis on uniform and non-uniform meshes as well as for different polynomial orders and continuity of the basis. Read More

As illustrated via numerical experiments with an implementation in Spark (the popular platform for distributed computation), randomized algorithms solve two ubiquitous problems: (1) calculating a full principal component analysis or singular value decomposition of a highly rectangular matrix, and (2) calculating a low-rank approximation in the form of a singular value decomposition to an arbitrary matrix. Several optimizations to recently introduced methods yield results that are uniformly superior to those of the stock implementations. Read More

This work is devoted to the design of interior penalty discontinuous Galerkin (dG) schemes that preserve maximum principles at the discrete level for the steady transport and convection-diffusion problems and the respective transient problems with implicit time integration. Monotonic schemes that combine explicit time stepping with dG space discretization are very common, but the design of such schemes for implicit time stepping is rare, and it had only been attained so far for 1D problems. The proposed scheme is based on an artificial diffusion that linearly depends on a shock detector that identifies the troublesome areas. Read More

Owing to their low-complexity iterations, Frank-Wolfe (FW) solvers are well suited for various large-scale learning tasks. When block-separable constraints are also present, randomized FW has been shown to further reduce complexity by updating only a fraction of coordinate blocks per iteration. In this context, the present work develops feasibility-ensuring step sizes, and provably convergent randomized block Frank-Wolfe (RB-FW) solvers that are flexible in selecting the number of blocks to update per iteration. Read More

This work develops a parallelized algorithm to compute the dynamic mode decomposition (DMD) on a graphics processing unit using the streaming method of snapshots singular value decomposition. This allows the algorithm to operate efficiently on streaming data by avoiding redundant inner-products as new data becomes available. In addition, it is possible to leverage the native compressed format of many data streams, such as HD video and computational physics codes that are represented sparsely in the Fourier domain, to massively reduce data transfer from CPU to GPU and to enable sparse matrix multiplications. Read More

The Kaczmarz method is an iterative algorithm for solving systems of linear equalities and inequalities, that iteratively projects onto these constraints. Recently, Strohmer and Vershynin [J. Fourier Anal. Read More

Laguerre and Laguerre-type polynomials are orthogonal polynomials on the interval $[0,\infty)$ with respect to a weight function of the form $w(x) = x^{\alpha} e^{-Q(x)}, Q(x) = \sum_{k=0}^m q_k x^k, \alpha > -1, q_m > 0$. The classical Laguerre polynomials correspond to $Q(x)=x$. The computation of higher-order terms of the asymptotic expansions of these polynomials for large degree becomes quite complicated, and a full description seems to be lacking in literature. Read More

In this paper, an efficient divide-and-conquer (DC) algorithm is proposed for the symmetric tridiagonal matrices based on ScaLAPACK and the hierarchically semiseparable (HSS) matrices. HSS is an important type of rank-structured matrices.Most time of the DC algorithm is cost by computing the eigenvectors via the matrix-matrix multiplications (MMM). Read More

We present a polynomial multigrid method for the nodal interior penalty formulation of the Poisson equation on three-dimensional Cartesian grids. Its key ingredient is a weighted overlapping Schwarz smoother operating on element-centered subdomains. The MG method reaches superior convergence rates corresponding to residual reductions of about two orders of magnitude within a single V(1,1) cycle. Read More

A stochastic iterative algorithm approximating second-order information using von Neumann series is discussed. We present an improved and simplified convergence analysis, in contrast to a similar algorithm and its analysis, LISSA. The algorithm is primarily suitable for training large scale linear models, where the number of data points is very large. Read More

The Fast Fourier Transform (FFT) is an algorithm of paramount importance in signal processing as it allows to apply the Fourier transform in O(n log n) instead of O(n 2) arithmetic operations. Graph Signal Processing (GSP) is a recent research domain that generalizes classical signal processing tools, such as the Fourier transform, to situations where the signal domain is given by any arbitrary graph instead of a regular grid. Today, there is no method to rapidly apply graph Fourier transforms. Read More

Functions of one or more variables are usually approximated with a basis; a complete, linearly independent set of functions that spans an appropriate function space. The topic of this paper is the numerical approximation of functions using the more general notion of frames; complete systems that are generally redundant but provide stable infinite representations. While frames are well-known in image and signal processing, coding theory and other areas of mathematics, their use in numerical analysis is far less widespread. Read More

Nanoindentation involves probing a hard diamond tip into a material, where the load and the displacement experienced by the tip is recorded continuously. This load-displacement data is a direct function of material's innate stress-strain behavior. Thus, theoretically it is possible to extract mechanical properties of a material through nanoindentation. Read More

We consider a new fictitious domain approach of higher order accuracy. To implement Dirichlet conditions we apply the classical Nitsche method combined with a facet-based stabilization (ghost penalty). Both techniques are combined with a higher order isoparametric finite element space which is based on a special mesh transformation. Read More