Bounded Condition Number Research Articles

Numerical algorithm for estimating a conditioned symmetric positive definite matrix under constraints

SummaryWe present RCO (regularized Cholesky optimization): a numerical algorithm for finding a symmetric positive definite (PD) matrix with a bounded condition number that minimizes an objective function. This task arises when estimating a covariance matrix from noisy data or due to model constraints, which can cause spurious small negative eigenvalues. A special case is the problem of finding the nearest well‐conditioned PD matrix to a given matrix. RCO explicitly optimizes the entries of the Cholesky factor. This requires solving a regularized non‐linear, non‐convex optimization problem, for which we apply Newton‐CG and exploit the Hessian's sparsity. The regularization parameter is determined via numerical continuation with an accuracy‐conditioning trade‐off criterion. We apply RCO to our motivating educational measurement application of estimating the covariance matrix of an empirical best linear prediction (EBLP) of school growth scores. We present numerical results for two empirical datasets, state and urban. RCO outperforms general‐purpose near‐PD algorithms, obtaining ‐smaller matrix reconstruction bias and smaller EBLP estimator mean‐squared error. It is in fact the only algorithm that solves the right minimization problem, which strikes a balance between the objective function and the condition number. RCO can be similarly applied to the stable estimation of other posterior means. For the task of finding the nearest PD matrix, RCO yields similar condition numbers to near‐PD methods, but provides a better overall near‐null space.

Quasi-local and frequency-robust preconditioners for the Helmholtz first-kind integral equations on the disk

We propose preconditioners for the Helmholtz scattering problems by a planar, disk-shaped screen in ℝ3. Those preconditioners are approximations of the square-roots of some partial differential operators acting on the screen. Their matrix-vector products involve only a few sparse system resolu- tions and can thus be evaluated cheaply in the context of iterative methods. For the Laplace equation (i.e. for the wavenumber k = 0) with Dirichlet condition on the disk and on regular meshes, we prove that the preconditioned linear system has a bounded condition number uniformly in the mesh size. We further provide numerical evidence indicating that the preconditioners also perform well for large values of k and on locally refined meshes.

On the stability of unevenly spaced samples for interpolation and quadrature

Unevenly spaced samples from a periodic function are common in signal processing and can often be viewed as a perturbed equally spaced grid. In this paper, the question of how the uneven distribution of the samples impacts the quality of interpolation and quadrature is analyzed. Starting with equally spaced nodes on $$[-\pi ,\pi )$$ with grid spacing h, suppose the unevenly spaced nodes are obtained by perturbing each uniform node by an arbitrary amount $$\le \alpha h$$ , where $$0 \le \alpha < 1/2$$ is a fixed constant. A discrete version of the Kadec-1/4 theorem is proved, which states that the nonuniform discrete Fourier transform associated with perturbed nodes has a bounded condition number independent of h, for any $$\alpha <1/4$$ . Then, it is shown that unevenly spaced quadrature rules converge for all continuous functions and interpolants converge uniformly for all differentiable functions whose derivative has bounded variation when $$0\le \alpha <1/4$$ . Though, quadrature rules at perturbed nodes can have negative weights for any $$\alpha >0$$ , a bound on the absolute sum of the quadrature weights is provided, which shows that perturbed equally spaced grids with small $$\alpha $$ can be used without numerical woes. While the proof techniques work primarily when $$0 \le \alpha < 1/4$$ , it is shown that a small amount of oversampling extends our results to the case when $$1/4\le \alpha <1/2$$ .

Global Linear Convergence of Evolution Strategies on More than Smooth Strongly Convex Functions

Evolution strategies (ESs) are zeroth-order stochastic black-box optimization heuristics invariant to monotonic transformations of the objective function. They evolve a multivariate normal distribution, from which candidate solutions are generated. Among different variants, CMA-ES is nowadays recognized as one of the state-of-the-art zeroth-order optimizers for difficult problems. Despite ample empirical evidence that ESs with a step-size control mechanism converge linearly, theoretical guarantees of linear convergence of ESs have been established only on limited classes of functions. In particular, theoretical results on convex functions are missing, where zeroth-order and also first-order optimization methods are often analyzed. In this paper, we establish almost sure linear convergence and a bound on the expected hitting time of an ES family, namely, the $(1+1)_\kappa$-ES, which includes the (1+1)-ES with (generalized) one-fifth success rule and an abstract covariance matrix adaptation with bounded condition number, on a broad class of functions. The analysis holds for monotonic transformations of positively homogeneous functions and of quadratically bounded functions, the latter of which particularly includes monotonic transformation of strongly convex functions with Lipschitz continuous gradient. As far as the authors know, this is the first work that proves linear convergence of ES on such a broad class of functions.

Exact Solution of New Magnetic Current Based Surface-Volume-Surface EFIE and Analysis of Its Spectral Properties

A novel magnetic current based Surface-Volume-Surface Electric Field Integral Equation (SVS-EFIE-M) is presented for the problem of scattering on homogeneous non-magnetic dielectric objects. The exact Galerkin Method of Moments (MoM) utilizing both the rotational and irrotational vector spherical harmonics as orthogonal basis and test functions according to the Helmholtz decomposition is implemented to solve SVS-EFIE-M analytically for the case of dielectric sphere excited by an electric dipole. The field throughout the sphere is evaluated and compared against the exact classical Mie series solution. The two are shown to agree to 12 digits of accuracy upon a sufficient number of basis/test functions taken in the MoM solution and the Mie series expansion. This exact solution validates the rigorous nature of the new SVS-EFIE-M formulation. It also reveals the spectral properties of its individual operators, their products and their linear combination. The spectrum of the MoM impedance matrix is also obtained. It is shown that upon choosing basis and test functions in <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"><tex-math notation="LaTeX">$L^{2}(S)$</tex-math></inline-formula> space and evaluating testing inner products in the same space, the MoM impedance matrix features bounded condition number with increasing order of discretization and/or at low frequencies. This makes the proposed SVS-EFIE-M formulation free of oversampling and low-frequency breakdowns giving it advantage both over its SVS-EFIE-J predecessor and classical double-source integral equations such as PMCHWT, Muller, and others suffering from this type of numerical instabilities inherent to their inferior spectral properties.

Analytic Solution of Surface–Volume–Surface Electric Field Integral Equation on Dielectric Sphere and Analysis of Its Spectral Properties

The exact solution of the surface–volume–surface electric field integral equation (SVS-EFIE) is presented for the problem of radiation in the vicinity of the homogeneous dielectric sphere. The solution is obtained via the Galerkin method of moments (MoM) utilizing rotational and irrotational complete sets of orthogonal vector spherical harmonics as basis and test functions according to the Helmholtz decomposition. In the case of radial or tangential electric dipole radiation, the electric field throughout the sphere evaluated via the analytic MoM solution of the SVS-EFIE is compared against the exact classical Mie series solution. The two are shown to agree to 12 digits of accuracy upon a sufficient number of basis/test functions taken in the MoM solution and the Mie series expansion. The exact solution confirms and validates the rigorous nature of the SVS-EFIE formulation. It also reveals the spectral properties of its individual operators, their products, and their linear combination, as well as the spectrum of the MoM impedance matrix. It is shown that, upon choosing basis and test functions in the Sobolev space <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$H_{div}^{1/2}(S)$ </tex-math></inline-formula> and performing testing inner products evaluation in the space <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$H_{div}^{-1/2}(S)$ </tex-math></inline-formula> , with <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$S$ </tex-math></inline-formula> being the surface of the sphere, the MoM impedance matrix features bounded condition number with increasing order of discretization similar with analogous exact MoM solution of the surface EFIE (S-EFIE) on the perfectly electrically conducting (PEC) sphere.

Linear systems arising in interior methods for convex optimization: a symmetric formulation with bounded condition number

We provide eigenvalues bounds for a new formulation of the step equations in interior methods for convex quadratic optimization. The matrix of our formulation, named , has bounded condition number, converges to a well-defined limit under strict complementarity, and has the same size as the traditional, ill-conditioned, saddle-point formulation. We evaluate the performance in the context of a Matlab object-oriented implementation of PDCO, an interior-point solver for minimizing a smooth convex function subject to linear constraints. The main benefit of our implementation, named PDCOO, is to separate the logic of the interior-point method from the formulation of the system used to compute a step at each iteration and the method used to solve the system. Thus, PDCOO allows easy addition of a new system formulation and/or solution method for experimentation. Our numerical experiments indicate that the formulation has the same storage requirements as the traditional ill-conditioned saddle-point formulation, and its condition is often more favourable than the unsymmetric block formulation.

Multilevel decompositions and norms for negative order Sobolev spaces

We consider multilevel decompositions of piecewise constants on simplicial meshes that are stable in H − s H^{-s} for s ∈ ( 0 , 1 ) s\in (0,1) . Proofs are given in the case of uniformly and locally refined meshes. Our findings can be applied to define local multilevel diagonal preconditioners that lead to bounded condition numbers (independent of the mesh-sizes and levels) and have optimal computational complexity. Furthermore, we discuss multilevel norms based on local (quasi-)projection operators that allow the efficient evaluation of negative order Sobolev norms. Numerical examples and a discussion on several extensions and applications conclude this article.

New Magnetic Current Based Surface-Volume-Surface Electric Field Integral Equation For Solution of Scattering Problems on Dielectric Objects

New Surface-Volume-Surface Electric Field Integral Equation (SVS-EFIE) based on magnetic current single source field representation (SVS-EFIE-M) is proposed for the problems of scattering on homogeneous non-magnetic dielectric objects. It's analytic MoM solution for a dielectric sphere is compared against the Mie series expansion to validate the exact nature of this new magnetic current based SVS-EFIE formulation. This analytic solution also allows to rigorously study spectral properties of the SVS-EFIE- M operators and to determine proper functional spaces in which MoM solutions of it should be sought in order to maintain bounded condition number of the pertinent matrix equation.

A parallel-in-time two-sided preconditioning for all-at-once system from a non-local evolutionary equation with weakly singular kernel

In this paper, we study a parallel-in-time (PinT) algorithm for all-at-once system from a non-local evolutionary equation with weakly singular kernel where the temporal term involves a non-local convolution with a weakly singular kernel and the spatial term is the usual Laplacian operator with variable coefficients. Such a problem has been intensively studied in recent years thanks to the numerous real world applications. However, due to the non-local property of the time evolution, solving the equation in PinT manner is difficult. We propose to use a two-sided preconditioning technique for the all-at-once discretization of the equation. Our preconditioner is constructed by replacing the variable diffusion coefficients with a constant coefficient to obtain a constant-coefficient all-at-once matrix. We split a square root of the constant Laplacian operator out of the constant-coefficient all-at-once matrix as a right preconditioner and take the remaining part as a left preconditioner, which constitutes our two-sided preconditioning. Exploiting the diagonalizability of the constant-Laplacian matrix and the triangular Toeplitz structure of the temporal discretization matrix, we obtain efficient representations of inverses of the right and the left preconditioners, because of which the iterative solution can be fast updated in a PinT manner. Theoretically, the condition number of the two-sided preconditioned matrix is proven to be uniformly bounded by a constant independent of the matrix size. To the best of our knowledge, for the non-local evolutionary equation with variable coefficients, this is the first attempt to develop a PinT preconditioning technique that has fast and exact implementation and that the corresponding preconditioned system has a uniformly bounded condition number. Numerical results are reported to confirm the efficiency of the proposed two-sided preconditioning technique.

An inverse problem formulation of the immersed‐boundary method

AbstractWe formulate the immersed‐boundary method (IBM) as an inverse problem. A control variable is introduced on the boundary of a larger domain that encompasses the target domain. The optimal control is the one that minimizes the mismatch between the state and the desired boundary value along the immersed target‐domain boundary. We begin by investigating a naïve problem formulation that we show is ill‐posed: in the case of the Laplace equation, we prove that the solution is unique, but it fails to depend continuously on the data; for the linear advection equation, even solution uniqueness fails to hold. These issues are addressed by two complimentary strategies. The first strategy is to ensure that the enclosing domain tends to the true domain, as the mesh is refined. The second strategy is to include a specialized parameter‐free regularization that is based on penalizing the difference between the control and the state on the boundary. The proposed inverse IBM is applied to the diffusion, advection, and advection‐diffusion equations using a high‐order discontinuous Galerkin discretization. The numerical experiments demonstrate that the regularized scheme achieves optimal rates of convergence and that the reduced Hessian of the optimization problem has a bounded condition number, as the mesh is refined.

Point Cloud Denoising via Feature Graph Laplacian Regularization.

Point cloud is a collection of 3D coordinates that are discrete geometric samples of an object's 2D surfaces. Imperfection in the acquisition process means that point clouds are often corrupted with noise. Building on recent advances in graph signal processing, we design local algorithms for 3D point cloud denoising. Specifically, we design a signal-dependent feature graph Laplacian regularizer (SDFGLR) that assumes surface normals computed from point coordinates are piecewise smooth with respect to a signal-dependent graph Laplacian matrix. Using SDFGLR as a signal prior, we formulate an optimization problem with a general 'p-norm fidelity term that can explicitly remove only two types of additive noise: small but non-sparse noise like Gaussian (using '2 fidelity term) and large but sparser noise like Laplacian (using '1 fidelity term). To establish a linear relationship between normals and 3D point coordinates, we first perform bipartite graph approximation to divide the point cloud into two disjoint node sets (red and blue). We then optimize the red and blue nodes' coordinates alternately. For '2-norm fidelity term, we iteratively solve an unconstrained quadratic programming (QP) problem, efficiently computed using conjugate gradient with a bounded condition number to ensure numerical stability. For '1-norm fidelity term, we iteratively minimize an '1-'2 cost function using accelerated proximal gradient (APG), where a good step size is chosen via Lipschitz continuity analysis. Finally, we propose simple mean and median filters for flat patches of a given point cloud to estimate the noise variance given the noise type, which in turn is used to compute a weight parameter trading off the fidelity term and signal prior in the problem formulation. Extensive experiments show state-of-the-art denoising performance among local methods using our proposed algorithms.

Positivity preserving finite difference methods for Poisson–Nernst–Planck equations with steric interactions: Application to slit-shaped nanopore conductance

To study ion transport in electrolyte solutions, we propose numerical methods for a modified Poisson–Nernst–Planck model with ionic steric effects (SPNP). Positivity preserving schemes based on harmonic-mean approximations are proposed on a nonuniform mesh for the spatial discretization of the SPNP equations. Both explicit and semi-implicit discretization are considered in time. Numerical analysis shows that explicit forward Euler and semi-implicit trapezoidal discretization lead to schemes that maintain fully discrete solution positivity under a constraint on a mesh ratio, while the semi-implicit backward Euler discretization maintain fully discrete solution positivity unconditionally. We further study the condition numbers of the matrix associated with the semi-implicit backward Euler discretization, and establish an upper bound on condition numbers, indicating that the developed discretization based on harmonic-mean approximations effectively solves the issue of the presence of large condition numbers when using the Slotboom-type variables. Further numerical simulations confirm the analysis results on accuracy, positivity, and bounded condition numbers. The proposed schemes are also applied to study practical applications, such as the impact of self and cross steric interactions on ion distribution and rectifying behavior in a slit-shaped nanopore with surface charges. Possible extensions of the numerical methods to other modified PNP models with ion correlations and steric effects are also discussed.

Galerkin method with new quadratic spline wavelets for integral and integro-differential equations

The paper is concerned with the wavelet-Galerkin method for the numerical solution of Fredholm linear integral equations and second-order integro-differential equations. We propose a construction of a quadratic spline-wavelet basis on the unit interval, such that the wavelets have three vanishing moments and the shortest support among such wavelets. We prove that this basis is a Riesz basis in the space L 2 0 , 1 . We adapt the basis to homogeneous Dirichlet boundary conditions, and using a tensor product we construct a wavelet basis on the hyperrectangle. We use the wavelet-Galerkin method with the constructed bases for solving integral and integro-differential equations, and we show that the matrices arising from discretization have uniformly bounded condition numbers and that they can be approximated by sparse matrices. We present numerical examples and compare the results with the Galerkin method using other quadratic spline wavelet bases and other methods.

Frames and Numerical Approximation

Functions of one or more variables are usually approximated with a basis: a complete, linearly-independent system of functions that spans a suitable function space. The topic of this paper is the numerical approximation of functions using the more general notion of frames: that is, complete systems that are generally redundant but provide infinite representations with bounded coefficients. While frames are well-known in image and signal processing, coding theory and other areas of applied mathematics, their use in numerical analysis is far less widespread. Yet, as we show via a series of examples, frames are more flexible than bases, and can be constructed easily in a range of problems where finding orthonormal bases with desirable properties (rapid convergence, high resolution power, etc.) is difficult or impossible. A key concern when using frames is that computing a best approximation requires solving an ill-conditioned linear system. Nonetheless, we construct a frame approximation via regularization with bounded condition number (with respect to perturbations in the data), and which approximates any function up to an error of order $\sqrt{\epsilon}$, or even of order $\epsilon$ with suitable modifications. Here $\epsilon$ is a threshold value that can be chosen by the user. Crucially, rate of decay of the error down to this level is determined by the existence of approximate representations of $f$ in the frame possessing small-norm coefficients. We demonstrate the existence of such representations in all of our examples. Overall, our analysis suggests that frames are a natural generalization of bases in which to develop numerical approximation. In particular, even in the presence of severely ill-conditioned linear systems, the frame condition imposes sufficient mathematical structure in order to give rise to accurate, well-conditioned approximations.

Efficient preconditioner of one-sided space fractional diffusion equation

In this paper, we propose an efficient preconditioner for the linear systems arising from the one-sided space fractional diffusion equation with variable coefficients. The shifted Gr $$\ddot{\mathrm{u}}$$ nwald formula is employed to discretize the one-sided Riemann–Liouville fractional derivative. The matrix structure of resulting linear systems is Toeplitz-like, which is a summation of an identity matrix and a diagonal-times-nonsymmetric-Toeplitz matrix. A diagonal-times-nonsymmetric-Toeplitz preconditioner is proposed to reduce the condition number of the Toeplitz-like matrix, where the diagonal part comes from the variable coefficients and the nonsymmetric Toeplitz part comes from the Riemann–Liouville derivative. Theoretically, we show that the condition number of the preconditioned matrix is uniformly bounded by a constant independent of discretization step-sizes under certain assumptions on the coefficient function. Due to the uniformly bounded condition number, the Krylov subspace method for the preconditioned linear systems converges linearly and independently on discretization step-sizes. Numerical results are reported to show the efficiency of the proposed preconditioner and to demonstrate its superiority over other tested preconditioners.

Quadratic spline wavelets with short support satisfying homogeneous boundary conditions

In this paper, we construct a new quadratic spline-wavelet basis on the interval and on the unit square satisfying homogeneous Dirichlet boundary conditions of the first order. The wavelets have one vanishing moment and the shortest support among quadratic spline wavelets with at least one vanishing moment adapted to the same type of boundary conditions. The stiffness matrices arising from the discretization of the second-order elliptic problems using the constructed wavelet basis have uniformly bounded condition numbers, and the condition numbers are small. We present some quantitative properties of the constructed basis. We provide numerical examples to show that the Galerkin method and the adaptive wavelet method using our wavelet basis require fewer iterations than methods with other quadratic spline wavelet bases. Moreover, due to the small support of the wavelets, when using these methods with the new wavelet basis, the system matrix is sparser, and thus one iteration requires a smaller number of floating point operations than for other quadratic spline wavelet bases.

Spectrally Accurate and Well-Posed Jacobi Galerkin Method for Multiple Delay Pantograph Integral Equations

Operational Jacobi Galerkin method is developed and analyzed to solve multiple delay pantograph integral equations. This algorithm is an application of the matrix vector product approach in Galerkin formulation of the problem when the classical Jacobi polynomials are employed as natural basis functions. Our scheme has two main advantages. First, it produces a well-posed approximate method in the sense that the algebraic form of the Galerkin discretization of the problem has an upper triangular structure with a bounded condition number which can be solved directly by the forward substitution method and second, it has a potential convergence property such that its numerical error decay with the familiar spectral rate of convergence. Finally, some numerical examples are used to illustrate the efficiency and accuracy of the method.

An Inversion-Free Estimating Equations Approach for Gaussian Process Models

ABSTRACTOne of the scalability bottlenecks for the large-scale usage of Gaussian processes is the computation of the maximum likelihood estimates of the parameters of the covariance matrix. The classical approach requires a Cholesky factorization of the dense covariance matrix for each optimization iteration. In this work, we present an estimating equations approach for the parameters of zero-mean Gaussian processes. The distinguishing feature of this approach is that no linear system needs to be solved with the covariance matrix. Our approach requires solving an optimization problem for which the main computational expense for the calculation of its objective and gradient is the evaluation of traces of products of the covariance matrix with itself and with its derivatives. For many problems, this is an O(nlog n) effort, and it is always no larger than O(n2). We prove that when the covariance matrix has a bounded condition number, our approach has the same convergence rate as does maximum likelihood in that the Godambe information matrix of the resulting estimator is at least as large as a fixed fraction of the Fisher information matrix. We demonstrate the effectiveness of the proposed approach on two synthetic examples, one of which involves more than 1 million data points.

Large-scale bounded distortion mappings

We propose an efficient algorithm for computing large-scale bounded distortion maps of triangular and tetrahedral meshes. Specifically, given an initial map, we compute a similar map whose differentials are orientation preserving and have bounded condition number. Inspired by alternating optimization and Gauss-Newton approaches, we devise a first order method which combines the advantages of both. On the one hand, its iterations are as computationally efficient as those of alternating optimization. On the other hand, it enjoys preferable convergence properties, associated with Gauss-Newton like approaches. We demonstrate the utility of the proposed approach in efficiently solving geometry processing problems, focusing on challenging large-scale problems.