Search

An earlier paper proved the convergence of a variable stepsize Bregman operator splitting algorithm (BOSVS) for minimizing ø(Bu) + H(u), where H and ø are convex functions, and ø is possibly nonsmooth. The algorithm was shown to be relatively efficient when applied to partially parallel magnetic resonance image reconstruction problems. In this paper, the convergence rate of BOSVS is analyzed. When H(u) = ‖Au - f ‖2, where A is a matrix, it is shown that for an ergodic approximation uk obtained by averaging k BOSVS iterates, the error in the objective value ø(Buk)+H(uk) is O(1/k). When the optimization problem has a unique solution u*, we obtain the estimate ‖uk-u*‖ = O(1/ √ k). The theoretical analysis is compared to observed convergence rates for partially parallel magnetic resonance image reconstruction problems where A is a large dense ill-conditioned matrix. [ABSTRACT FROM AUTHOR]

A framework, PhaseLift, was recently proposed to solve the phase retrieval problem. In this framework, the problem is solved by optimizing a cost function over the set of complex Hermitian positive semidefinite matrices. This approach to phase retrieval motivates a more general consideration of optimizing cost functions on semidefinite Hermitian matrices where the desired minimizers are known to have low rank. This paper considers an approach based on an alternative cost function defined on a union of appropriate manifolds. It is related to the original cost function in a manner that preserves the ability to find a global minimizer and is significantly more efficient computationally. A rank-based optimality condition for stationary points is given and optimization algorithms based on state-of-the-art Riemannian optimization and dynamically reducing rank are proposed. Empirical evaluations are performed using the PhaseLift problem. The new approach is shown to be an effective method of phase retrieval with computational efficiency increased substantially compared to a state-of-the-art algorithm, the Wirtinger flow algorithm. A preliminary version of this paper can be found in [W. Huang, K. A. Gallivan, and X. Zhang, Procedia Comput. Sci., 80 (2016), pp. 1125-1134]. [ABSTRACT FROM AUTHOR]

A systematic mathematical framework for the study of numerical algorithms would allow comparisons, facilitate conjugacy arguments, and enable the discovery of improved, accelerated, data-driven algorithms. Over the course of the past century, the Koopman operator has provided a mathematical framework for the study of dynamical systems which facilitates conjugacy arguments and can provide efficient reduced descriptions. More recently, numerical approximations of the operator have enabled the analysis of a large number of deterministic and stochastic dynamical systems in a completely data-driven, essentially equation-free pipeline. Discrete- or continuous-time numerical algorithms (integrators, nonlinear equation solvers, optimization algorithms are themselves dynamical systems. In this paper, we use this insight to leverage the Koopman operator framework in the data-driven study of such algorithms and discuss benefits for analysis and acceleration of numerical computation. For algorithms acting on high-dimensional spaces by quickly contracting them toward low-dimensional manifolds, we demonstrate how basis functions adapted to the data help to construct efficient reduced representations of the operator. Our illustrative examples include the gradient descent and Nesterov optimization algorithms as well as the Newton-Raphson algorithm. [ABSTRACT FROM AUTHOR]

Let F be a compact set of a Banach space X. This paper analyzes the "generalized empirical interpolation method," which, given a function f ∊ F, builds an interpolant Jn[f] in an n dimensional subspace Xn ⊂ X with the knowledge of n outputs (σi(f))ni=1, where σi ∊ X' and X' is the dual space of X. The space Xn is built with a greedy algorithm that is adapted to F in the sense that it is generated by elements of F itself. The algorithm also selects the linear functionals (σi)ni=1 from a dictionary Σ ⊂ X'. In this paper, we study the interpolation error maxf∊F ‖f - Jn[f]‖X by comparing it with the best possible performance on an n dimensional space, i.e., the Kolmogorov n-width of F in X, dn(F,X). For polynomial or exponential decay rates of dn(F,X), we prove that the interpolation error has the same behavior modulo the norm of the interpolation operator. Sharper results are obtained in the case where X is a Hilbert space. [ABSTRACT FROM AUTHOR]

The generalized null space decomposition (GNSD) is a unitary reduction of a general matrix A of order n to a block upper triangular form that reveals the structure of the Jordan blocks of A corresponding to a zero eigenvalue. The reduction was introduced by Kublanovskaya. It was extended first by Ruhe and then by Golub and Wilkinson, who based the reduction on the singular value decomposition. Unfortunately, if A has large Jordan blocks, the complexity of these algorithms can approach the order of n4. This paper presents an alternative algorithm, based on repeated updates of a QR decomposition of A, that is guaranteed to be of order n³. Numerical experiments confirm the stability of this algorithm, which turns out to produce essentially the same form as that of Golub and Wilkinson. The effect of errors in A on the ability to recover the original structure is investigated empirically. Several applications are discussed, including the computation of the Drazin inverse. [ABSTRACT FROM AUTHOR]

We introduce a new algorithm of proper generalized decomposition (PGD) for parametric symmetric elliptic partial differential equations. For any given dimension, we prove the existence of an optimal subspace of at most that dimension which realizes the best approximation---in the mean parametric norm associated to the elliptic operator---of the error between the exact solution and the Galerkin solution calculated on the subspace. This is analogous to the best approximation property of the proper orthogonal decomposition (POD) subspaces, except that in our case the norm is parameter-dependent. We apply a deflation technique to build a series of approximating solutions on finite-dimensional optimal subspaces, directly in the online step, and we prove that the partial sums converge to the continuous solution in the mean parametric elliptic norm. We show that the standard PGD for the considered parametric problem is strongly related to the deflation algorithm introduced in this paper. This opens the possibility of computing the PGD expansion by directly solving the optimization problems that yield the optimal subspaces. [ABSTRACT FROM AUTHOR]

The Euler genus of a graph is a fundamental and well-studied parameter in graph theory and topology. Computing it has been shown to be NP-hard by Thomassen [J. Algorithms, 10 (1989), pp. 568{576; J. Combin. Theory, Ser. B, 57 (1993), pp. 196{206], and it is known to be fixed-parameter tractable. However, the approximability of the Euler genus is wide open. While the existence of an O(1)-approximation is not ruled out, only an O( p n)-approximation [J. Chen, S. P. Kanchi, and A. Kanevsky, Inform. Process. Lett., 61 (1997), pp. 317{322] is known even in bounded-degree graphs. In this paper we give a polynomial-time algorithm which, given a boundeddegree graph of Euler genus g, computes a drawing in a surface of Euler genus gO(1) · logO(1) n. Combined with the upper bound from [J. Chen, S. P. Kanchi, and A. Kanevsky, Inform. Process. Lett., 61 (1997), pp. 317{322], our result also implies a O(n1/2-α)-approximation for some constant α > 0. Using our algorithm for approximating the Euler genus as a subroutine, we obtain, in a uniform fashion, algorithms with approximation ratios of the form OPTO(1) · logO(1) n for several related problems on bounded-degree graphs. These include the problems of orientable genus, crossing number, and planar edge and vertex deletion. Our algorithm and proof of correctness for the crossing number problem are simpler compared to the long and difficult proof in the recent breakthrough by Chuzhoy [Proceedings of the ACM Symposium on Theory of Computing, 2011, pp. 303{312], while essentially obtaining a qualitatively similar result. For planar edge and vertex deletion problems our results are the first to obtain a bound of the form poly(OPT; log n). We also highlight some further applications of our results in the design of algorithms for graphs with small genus. Many such algorithms require that a drawing of the graph is given as part of the input. Our results imply that in several interesting cases, we can implement such algorithms even when the drawing is unknown. [ABSTRACT FROM AUTHOR]

Threshold dynamics and its extensions have proven useful in computing interfacial motions in applications as diverse as materials science and machine learning. Certain desirable properties of the algorithm, such as unconditional monotonicity in two-phase ows and gradient stability more generally, hinge on positivity properties of the convolution kernel and its Fourier transform. Recent developments in the analysis of this class of algorithms indicate that sometimes, as in the case of certain anisotropic curvature ows arising in materials science, these properties of the convolution kernel cannot be expected. Other applications, such as machine learning, would benefit from as great a level of exibility in choosing the convolution kernel as possible. In this paper, we establish certain desirable properties of threshold dynamics, such as gamma convergence of its associated energy, for a substantially wider class of kernels than has been hitherto possible. We also present variants of the algorithm that extend some of these properties to even wider classes of convolution kernels. [ABSTRACT FROM AUTHOR]

The roundoff-error-free (REF) LU and Cholesky factorizations, combined with the REF substitution algorithms, allow rational systems of linear equations to be solved exactly and efficiently by working entirely in integer arithmetic. These REF computational tools share two key properties: their constituent divisions are exact, and their matrix entries' bit-lengths (i.e., required number of storage bits) are bounded polynomially. However, as this paper explains, updating the REF factorizations via the common delete-insert-reduce approach turns out to be inefficient in terms of operand bit-length growth and increased computational effort. In fact, we prove that this inefficiency can eliminate the computational savings expected of the factorization update process. Hence, the current work develops REF update algorithms that differ significantly from their traditional counterparts. To achieve this, it introduces the frame matrix methodology for representing and applying the REF factorizations, and then it devises new REF operations and shortcuts for transitioning between "nearby" frame matrices. The presented push-and-swap update approach utilizes these transitions to preserve the special structure of the REF factorizations and to avoid additional growth in their matrix entries' bit-lengths; this update process requires the expected O(n²) operations. The featured updates are column addition, deletion, and replacement with respect to the REF LU factorization. We also prove that analogous row updates can be performed via the column updates and discuss special considerations for updating the REF Cholesky factorization. An accompanying set of computational tests attests that for fully dense random instances, the push-and-swap column replacement update for the REF LU factorization is on average an order of magnitude faster than the corresponding delete-insert-reduce update for the exact rational arithmetic LU factorization. Altogether, the REF update algorithms represent a promising preliminary step toward enhancing the basic solution validation subroutines employed by state-of-the-art mixed-precision methodologies for exact linear and mixed-integer programming. [ABSTRACT FROM AUTHOR]

The development of strong linearizations preserving whatever structure a matrix polynomial might possess has been a very active area of research in the last years, since such linearizations are the starting point of numerical algorithms for computing eigenvalues of structured matrix polynomials with the properties imposed by the considered structure. In this context, Hermitian matrix polynomials are one of the most important classes of matrix polynomials arising in applications and their real eigenvalues are of great interest. The sign characteristic is a set of signs attached to these real eigenvalues which is crucial for determining the behavior of systems described by Hermitian matrix polynomials and, therefore, it is desirable to develop linearizations that preserve the sign characteristic of these polynomials, but, at present, only one such linearization is known. In this paper, we present a complete characterization of all the Hermitian strong linearizations that preserve the sign characteristic of a given Hermitian matrix polynomial and identify several families of such linearizations that can be constructed very easily from the coefficients of the polynomial. [ABSTRACT FROM AUTHOR]

Solving feasibility problems is a central task in mathematics and the applied sciences. One particularly successful method is the Douglas--Rachford algorithm. In this paper, we provide many new conditions sufficient for finite convergence. Numerou s examples illustrate our results. [ABSTRACT FROM AUTHOR]

Fast algorithms for matrix multiplication, namely those that perform asymptotically fewer scalar operations than the classical algorithm, have been considered primarily of theoretical interest. Apart from Strassen's original algorithm, few fast algorithms have been efficiently im- plemented or used in practical applications. However, there exist many practical alternatives to Strassen's algorithm with varying performance and numerical properties. Fast algorithms are known to be numerically stable, but because their error bounds are slightly weaker than the classical al- gorithm, they are not used even in cases where they provide a performance benefit. We argue in this paper that the numerical sacrifice of fast algorithms, particularly for the typical use cases of practical algorithms, is not prohibitive, and we explore ways to improve the accuracy both theoret- ically and empirically. The numerical accuracy of fast matrix multiplication depends on properties of the algorithm and of the input matrices, and we consider both contributions independently. We generalize and tighten previous error analyses of fast algorithms and compare their properties. We discuss algorithmic techniques for improving the error guarantees from two perspectives: manipulat- ing the algorithms, and reducing input anomalies by various forms of diagonal scaling. Finally, we benchmark performance and demonstrate our improved numerical accuracy. [ABSTRACT FROM AUTHOR]

In this paper, we introduce the uncertain optimal control problem of determining a control that minimizes the expectation of an objective functional for a system with parameter uncertainty in both dynamics and objective. We present a computational framework for the numerical solution of this problem, wherein an independently drawn random sample is taken from the space of uncertain parameters, and the expectation in the objective functional is approximated by a sample average. The result is a sequence of approximating standard optimal control problems that can be solved using existing techniques. To analyze the performance of this computational framework, we develop necessary conditions for both the original and approximate problems and show that the approximation based on sample averages is consistent in the sense of Polak [Optimization: Algorithms and Consistent Approximations, Springer, New York, 1997]. This property guarantees that accumulation points of a sequence of global minimizers (stationary points) of the approximate problem are global minimizers (stationary points) of the original problem. We show that the uncertain optimal control problem can further be approximated in a consistent manner by a sequence of nonlinear programs under mild regularity assumptions. In numerical examples, we demonstrate that the framework enables the solution of optimal search and optimal ensemble control problems. [ABSTRACT FROM AUTHOR]

Max-plus methods have previously been used to solve deterministic control problems. The methods are based on max-plus (or min-plus) expansions and can yield curse-of-dimensionalityfree numerical methods. In this paper, we explore min-plus methods for continuous-time stochastic control on a finite-time horizon. We first approximate the original value function via timediscretization. By generalizing the min-plus distributive property to continuum spaces, we obtain an algorithm for recursive computation of the time-discretized values, which we refer to as the idempotent distributed dynamic programming principle (IDDPP). Under the IDDPP, the value function at each step can be represented as an infimum of functions in a certain class. This is a min-plus expansion for the value function. For the specific class of problems considered here, we see that the class can be taken as that consisting of the quadratic functions. A means for reducing the numbers of constituent quadratic functions is discussed. [ABSTRACT FROM AUTHOR]

We present a fast and approximate multifrontal solver for large-scale sparse linear systems arising from finite-difference, finite-volume, or finite-element discretization of high-frequency wave equations. The proposed solver leverages the butterfly algorithm and its hierarchical matrix extension for compressing and factorizing large frontal matrices via graph distance guided entry evaluation or randomized matrix-vector multiplication-based schemes. Complexity analysis and numerical experiments demonstrate O (N log² N) computation and O (N) memory complexity when applied to an N X N sparse system arising from 3D high-frequency Helmholtz and Maxwell problems. [ABSTRACT FROM AUTHOR]

Although mixed precision arithmetic has recently garnered interest for training dense neural networks, many other applications could benefit from the speedups and lower storage cost if applied appropriately. The growing interest in employing mixed precision computations motivates the need for rounding error analysis that properly handles behavior from mixed precision arithmetic. We develop mixed precision variants of existing Householder QR algorithms and show error analyses supported by numerical experiments. [ABSTRACT FROM AUTHOR]

The hp-adaptive approximation is formulated as an approximation problem on a full binary tree T, where for each of the leaves Δ an order p(Δ ) ≥ 1 is assigned in such a way that the sum of all orders p(Δ ) does not exceed N, which is called the complexity of the approximation. The leaves Δ correspond to the cells of the partition, while p(Δ ) is the dimension of the polynomial space used for the local approximation on Δ . Devising an incremental algorithm for near-best adaptive approximation for the problem of finding the best possible tree T and assignments p(Δ ) leads to building a construction that attaches a ghost tree with p(Δ ) leaves to each leaf Δ of T with p(Δ ) > 1. The resulting full binary tree τ has at most N leaves and can be used as a proxy of T for assembling hp-adaptive procedures. Under the standard assumptions about the local errors, we prove that the error of our approximation of complexity N is bounded by 2N-1 N-n+1σ n, where σ n, n ≤ N, is the error of the best possible approximation of complexity n. [ABSTRACT FROM AUTHOR]

Graphs with bounded highway dimension were introduced by Abraham et al. [Proceedings of SODA 2010, pp. 782{793] as a model of transportation networks. We show that any such graph can be embedded into a distribution over bounded treewidth graphs with arbitrarily small distortion. More concretely, given a weighted graph G = (V;E) of constant highway dimension, we show how to randomly compute a weighted graph H = (V;E0) that distorts shortest path distances of G by at most a 1 + ε factor in expectation, and whose treewidth is polylogarithmic in the aspect ratio of G. Our probabilistic embedding implies quasi-polynomial time approximation schemes for a number of optimization problems that naturally arise in transportation networks, including Travelling Salesman, Steiner Tree, and Facility Location. To construct our embedding for low highway dimension graphs we extend Talwar's [Proceedings of STOC 2004, pp. 281-290] embedding of low doubling dimension metrics into bounded treewidth graphs, which generalizes known results for Euclidean metrics. We add several nontrivial ingredients to Talwar's techniques, and in particular thoroughly analyze the structure of low highway dimension graphs. Thus we demonstrate that the geometric toolkit used for Euclidean metrics extends beyond the class of low doubling metrics. [ABSTRACT FROM AUTHOR]

We consider the a posteriori error estimation for an atomistic-to-continuum coupling scheme for a generic one-dimensional many-body next-nearest-neighbor interaction model in one dimension. We derive and rigorously prove the efficiency of the residual type estimator. We prove the equivalence between the residual type and the gradient recovery-type estimator in the continuum region and propose a (novel) hybrid a posteriori error estimator by combining the two types of estimators. Our numerical experiments illustrate the optimal convergence rate of the adaptive algorithms using these estimators, whose effectivity indices are also presented. [ABSTRACT FROM AUTHOR]

Consider the scattering of a time-harmonic point source from an unbounded rough surface. The surface is assumed to be different from a plane surface over a finite interval, which leads to the situation that the model problem can be reduced to an equivalent boundary value problem with compactly supported boundary data. Well-posedness is thus proved by employing the variational method with classical Fredholm theory. By constructing a modified near-field equation, we then propose a novel sampling-type method for reconstructing the shape and location of the local rough surface, yielding a fast imaging algorithm. Numerical experiments are presented to illustrate the effectiveness of the method. [ABSTRACT FROM AUTHOR]

In this work, we develop an O(N) implicit real space method in 1D and 2D for the Cahn-Hilliard (CH) and vector Cahn-Hilliard (VCH) equations, based on the method of lines transpose (MOLT) formulation. This formulation results in a semidiscrete time stepping algorithm, which we prove is gradient stable in the H-1 norm. The spatial discretization follows from dimensional splitting and an O(N) matrix-free solver, which applies fast convolution to the modified Helmholtz equation. We propose a novel factorization technique, in which fourth-order spatial derivatives are incorporated into the solver. The splitting error is included in the nonlinear fixed point iteration, resulting in a high-order, logically Cartesian (line-by-line) update. Our method is fast but not restricted to periodic boundaries like the fast Fourier transform (FFT). The basic solver is implemented using the backward Euler formulation, and we extend this to both backward difference formula (BDF) stencils, singly diagonal implicit Runge-Kutta (SDIRK), and spectral deferred correction (SDC) frameworks to achieve high orders of temporal accuracy. We demonstrate with numerical results that the CH and VCH equations maintain gradient stability in one and two spatial dimensions. We also explore time-adaptivity, so that meta-stable states and ripening events can be simulated both quickly and efficiently. [ABSTRACT FROM AUTHOR]

In tree-based adaptive mesh refinement, elements are partitioned between processes using a space-filling curve. The curve establishes an ordering between all elements that derive from the same root element, the tree. When representing complex geometries by connecting several trees, the roots of these trees form an unstructured coarse mesh. We present an algorithm to partition the elements of the coarse mesh such that (a) the fine mesh can be load-balanced to equal element counts per process regardless of the element-to-tree map, and (b) each process that holds fine mesh elements has access to the meta data of all relevant trees. As an additional feature, the algorithm partitions the meta data of relevant ghost (halo) trees as well. We develop in detail how each process computes the communication pattern for the partition routine without handshaking and with minimal data movement. We demonstrate the scalability of this approach on up to 917e3 MPI ranks and 371e9 coarse mesh elements, measuring run times of one second or less. [ABSTRACT FROM AUTHOR]

We present a multi-index quasi-Monte Carlo method for the solution of elliptic partial differential equations with random coefficients. By combining the multi-index sampling idea with randomly shifted rank-1 lattice rules, the algorithm constructs an estimator for the expected value of some functional of the solution. The efficiency of this new method is illustrated on a three-dimensional subsurface flow problem with lognormal diffusion coefficient with underlying Matérn covariance function. This example is particularly challenging because of the small correlation length considered, and thus the large number of uncertainties that must be included. We show numerical evidence that it is possible to achieve a cost inversely proportional to the requested tolerance on the root-mean-square error, for problems with a smoothly varying random field. [ABSTRACT FROM AUTHOR]

The integration of functions over measurable sets is a fundamental problem in computational science. When the measurable sets belong to high-dimensional spaces or the function is computationally complex, it may only be practical to estimate integrals based on weighted sums of function values from a finite collection of samples. Monte Carlo, quasi-Monte Carlo, and other (pseudo-)random schemes are common choices for determining a set of samples. These schemes are appealing for their conceptual ease and ability to circumvent, with various degrees of success, the so-called curse of dimensionality. However, convergence is often slow and described in terms of probability. We consider a general measure-theoretic interpretation of any sample based algorithm for numerically approximating an integral. A priori error bounds are proven that provide insight into defining adaptive sampling algorithms solving error optimization problems. We use these bounds to improve integral approximations for both forward and inverse problems. [ABSTRACT FROM AUTHOR]

We develop an a posteriori analysis of C0 interior penalty methods for the displacement obstacle problem of clamped Kirchhoff plates. We show t hat a residual based error estimator originally designed for C0 interior penalty methods for the bound ary value problem of clamped Kirchhoff plates can also be used for the obstacle problem. We obtain reliability and efficiency estimates for the error estimator and introduce an adaptive algorithm based on this error estimator. Numerical results indicate that the performance of the adaptive algorithm is optimal for both quadratic and cubic C0 interior penalty methods. [ABSTRACT FROM AUTHOR]

Starting with far field data of time-harmonic acoustic or electromagnetic waves radiated by a collection of compactly supported sources in two-dimensional free space, we develop criteria and algorithms for the recovery of the far field components radiated by each of the individual sources, and the simultaneous restoration of missing data segments. Although both parts of this inverse problem are severely ill-conditioned in general, we give precise conditions relating the wavelength, the diameters of the supports of the individual source components a nd the distances between them, and the size of the missing data segments, which guarantee that stab le recovery in the presence of noise is possible. The only additional requirement is that a priori information on the approximate location of the individual sources is available. We give analytic and numerical examples to confirm the sharpness of our results and to illustrate the performance of corresponding reconstruction algorithms, and we discuss consequences for stability and resolution in inverse source and inverse scattering problems. [ABSTRACT FROM AUTHOR]

We develop a fast phase retrieval method which can utilize a large class of local phaseless correlationbased measurements in order to recover a given signal x ∊ Cd (up to an unknown global phase) in near-linear O(d log4 d)-time. Accompanying theoretical analysis proves that the proposed algorithm is guaranteed to deterministically recover all signals x satisfying a natural atness (i.e., nonsparsity) condition for a particular choice of deterministic correlation-based measurements. A randomized version of these same measurements is then shown to provide nonuniform probabilistic recovery guarantees for arbitrary signals x ∊ Cd. Numerical experiments demonstrate the method's speed, accuracy, and robustness in practice|all code is made publicly available. In its simplest form, our proposed phase retrieval method employs a modified lifting scheme akin to the one utilized by the well-known PhaseLift algorithm. In particular, it interprets quadratic magnitude measurements of x as linear measurements of a restricted set of lifted variables, xixj, for |j - i| < δ≪ d. This leads to a linear system involving a total of (2δ - 1)d unknown lifted variables, all of which can then be solved for using only O(δd) measurements. Once these lifted variables, xixj for |j - i| < δ ≪ d, have been recovered, a fast angular synchronization method can then be used to propagate the local phase difference information they provide across the entire vector in order to estimate the (relative) phases of every entry of x. In addition, the lifted variables corresponding to x4xj = |xj|² automatically provide magnitude estimates for each entry, xj, of x. The proposed phase retrieval method then approximates x by carefully combining these entrywise phase and magnitude estimates. Finally, we conclude by developing an extension of the proposed method to the sparse phase retrieval problem; specifically, we demonstrate a sublinear-time compressive phase retrieval algorithm which is guaranteed to recover a given s-sparse vector x ∊ Cd with high probability in just O(s log5 s.log d)- time using only O(s log4 s.log d) magnitude measurements. In doing so we demonstrate the existence of compressive phase retrieval algorithms with near-optimal linear-in-sparsity runtime complexities. [ABSTRACT FROM AUTHOR]

This manuscript describes a technique for computing partial rank-revealing factorizations, such as a partial QR factorization or a partial singular value decomposition. The method takes as input a tolerance ε and an m×n matrix A and returns an approximate low-rank factorization of A that is accurate to within precision ε in the Frobenius norm (or some other easily computed norm). The rank k of the computed factorization (which is an output of the algorithm) is in all examples we examined very close to the theoretically optimal ε-rank. The proposed method is inspired by the Gram-Schmidt algorithm and has the same O(mnk) asymptotic flop count. However, the method relies on randomized sampling to avoid column pivoting, which allows it to be blocked, and hence accelerates practical computations by reducing communication. Numerical experiments demonstrate that the accuracy of the scheme is for every matrix that was tried at least as good as column-pivoted QR and is sometimes much better. Computational speed is also improved substantially, in particular on GPU architectures. [ABSTRACT FROM AUTHOR]

The Laplace Transform is frequently encountered in mathematics, physics, engineering, and other fields. However, the spectral properties of the Laplace Transform tend to complicate its numerical treatment; therefore, the closely related "Truncated" Laplace Transforms are often used in applications. We have constructed efficient algorithms for the evaluation of the left singular functions and singular values of the Truncated Laplace Transform. Together with the previously introduced algorithms for the evaluation of the right singular functions, these algorithms provide the singular value decomposition of the Truncated Laplace Transform. The resulting algorithms are applicable to all environments likely to be encountered in applications, including the evaluation of singular functions corresponding to extremely small singular values (e.g., 10-1000). [ABSTRACT FROM AUTHOR]

Hidden-variable resultant methods are a class of algorithms for solving multidimensional polynomial rootfinding problems. In two dimensions, when significant care is taken, they are competitive practical rootfinders. However, in higher dimensions they are known to miss zeros, calculate roots to low precision, and introduce spurious solutions. We show that the hidden variable resultant method based on the Cayley (Dixon or Bézout) matrix is inherently and spectacularly numerically unstable by a factor that grows exponentially with the dimension. We also show that the Sylvester matrix for solving bivariate polynomial systems can square the condition number of the problem. In other words, two popular hidden variable resultant methods are numerically unstable, and this mathematically explains the difficulties that are frequently reported by practitioners. Regardless of how the constructed polynomial eigenvalue problem is solved, severe numerical difficulties will be present. Along the way, we prove that the Cayley resultant is a generalization of Cramer's rule for solving linear systems and generalize Clenshaw's algorithm to an evaluation scheme for polynomials expressed in a degree-graded polynomial basis. [ABSTRACT FROM AUTHOR]

In this article we introduce a novel numerical method to solve the problem of optimal transport and the related elliptic Monge-Ampère equation. It is one of the few numerical algorithms capable of solving this problem efficiently with the proper transport boundary condition. The computation time scales well with the grid size and has the additional advantage that the target domain may be nonconvex. We present the method and several numerical experiments. [ABSTRACT FROM AUTHOR]

Compressive sensing has been widely applied to problems in signal and imaging processing. In this work, we present an algorithm for predicting optimal real-time compression rates for video. The video data we consider is spatially compressed during the acquisition process, unlike in many of the standard methods. Rather than temporally compressing the frames at a fixed rate, our algorithm adaptively predicts the compression rate given the behavior of a few previous compressed frames. The algorithm uses polynomial fitting and simple filters, making it computationally feasible and easy to implement in hardware. Based on numerical simulations of real videos, the algorithm is able to capture object motion and approximate dynamics within the compressed frames. The adaptive video compression improves the quality of the reconstructed video (as compared to an equivalent fixed rate compression scheme) by several dB of peak signal-to-noise ratio without increasing the amount of information stored, as seen in numerical simulations presented here. [ABSTRACT FROM AUTHOR]

This article details a novel numerical scheme to approximate gradient flows for optimal transport (i.e., Wasserstein) metrics. These flows have proved useful to tackle theoretically and numerically nonlinear diffusion equations that model, for instance, porous media or crowd evolutions. These gradient flows define a suitable notion of weak solutions for these evolutions and they can be approximated in a stable way using discrete flows. These discrete flows are implicit Euler time stepping according to the Wasserstein metric. A bottleneck of these approaches is the high computational load induced by the resolution of each step. Indeed, this corresponds to the resolution of a convex optimization problem involving a Wasserstein distance to the previous iterate. Following several recent works on the approximation of Wasserstein distances, we consider a discrete flow induced by an entropic regularization of the transportation coupling. This entropic regularization allows one to trade the initial Wasserstein fidelity term for a Kullback-Leibler divergence, which is easier to deal with numerically. We show how Kullback-Leibler first order proximal schemes and, in particular, Dykstra's algorithm, can be used to compute each step of the regularized flow. The resulting algorithm is both fast, parallelizable, and versatile, because it only requires multiplications by the Gibbs kernel e -c/γ, where c is the ground cost and γ > 0 the regularization strength. On Euclidean domains discretized on a uniform grid, this corresponds to a linear filtering (for instance, a Gaussian filtering when c is the squared Euclidean distance) which can be computed in nearly linear time. On more general domains, such as (possibly nonconvex) shapes or on manifolds discretized by a triangular mesh, following a recently proposed numerical scheme for optimal transport, this Gibbs kernel multiplication is approximated by a short-time heat diffusion. We show numerical illustrations of this method to approximate crowd motions with congestion on complicated domains as well as extensions to take into account interactions between multiple densities. [ABSTRACT FROM AUTHOR]

The infinitesimal generator matrix underlying a multidimensional Markov chain can be represented compactly by using sums of Kronecker products of small rectangular matrices. For such compact representations, analysis methods based on vector-Kronecker product multiplication need to be employed. When the factors in the Kronecker product terms are relatively dense, vector-Kronecker product multiplication can be performed efficiently by the shuffle algorithm. When the factors are relatively sparse, it may be more efficient to obtain nonzero elements of the generator matrix in Kronecker form on the fly and multiply them with corresponding elements of the vector. This work proposes a modification to the shuffle algorithm that multiplies relevant elements of the vector with submatrices of factors in which zero rows and columns are omitted. This approach avoids unnecessary floating-point operations that evaluate to zero during the course of the multiplication and possibly reduces the amount of memory used. Numerical experiments on a large number of models indicate that in many cases the modified shuffle algorithm performs a smaller number of floating-point operations than the shuffle algorithm and the algorithm that generates nonzeros on the fly, sometimes with a minimum number of floating-point operations and as little of memory possible. [ABSTRACT FROM AUTHOR]

We consider the solution of large linear systems with multiple right-hand sides using a block GMRES approach. We introduce a new algorithm that effectively handles the situation of almost rank deficient block generated by the block Arnoldi procedure and that enables the recycling of spectral information at restart. The first feature is inherited from an algorithm introduced by Robb'e and Sadkane [Linear Algebra Appl., 419 (2006), pp. 265-285], while the second one is obtained by extending the deflated restarting strategy proposed by Morgan [Appl. Numer. Math., 54 (2005), pp. 222-236]. Through numerical experiments, we show that the new algorithm combines efficiently the attractive numerical features of its two parents and outperforms them. [ABSTRACT FROM AUTHOR]