18 results
Filters
Settings
Now showing1 - 10 of 18
- Some of the metrics are blocked by yourconsent settings
Item type:Publication, Access status: without files , Verified error bounds for sparse systems part II: inertia-based bounds, least squares and nonlinear systemsVerification methods provide mathematically correct error bounds for the solution of a numerical problem. That includes the proof of solvability of the problem and often uniqueness of the solution within the computed bounds. There are many verification methods for standard problems in numerical analysis, including linear and nonlinear systems of equations, matrix decompositions, eigenproblems, local and global optimization, and ordinary and partial differential equations. Many of those verification methods are included in INTLAB, the MATLAB/Octave toolbox for reliable computing. Despite many efforts, the solution of general sparse linear systems was an open problem. In Part I of this note, we presented an algorithm computing verified bounds for the solution of general real or complex sparse linear systems with condition number up to the limit [Formula: see text] in double precision. That algorithm splits into three subalgorithms for symmetric positive definite, symmetric indefinite, and general input matrix [Formula: see text]. It is based on a mathematically correct lower bound on the smallest singular value [Formula: see text]. In this Part II, we use a method published by the author in 1995 based on the inertia of a symmetric matrix. A key point is, as in Part I, a factorization [Formula: see text] such that [Formula: see text] and [Formula: see text] have identical sets of singular values with the smallest one close to [Formula: see text]. Numerical evidence suggests that the method is often a little slower than that in Part I, however, more robust. That means, for some of the few cases where the method in Part I could not compute verified bounds successfully, the method in this Part II succeeded. Furthermore, we show how to compute inclusions with almost maximal accuracy for all entries, i.e., all bounds differ by few bits, that is, based on a fast method to compute accurate approximations and bounds for extremely ill-conditioned dot products using a very efficient MATLAB implementation. We show how to treat a system with nearly symmetric input matrix in the sense that [Formula: see text] is small compared to [Formula: see text]. That is useful if symmetry is lost due to rounding and may improve the performance significantly. We show that Algorithm “verifySparseNearSym1” in Part I of the note is correct for unsymmetric input matrix without change. Moreover, algorithms are given to compute verified error bounds for a least squares problem and an underdetermined system of linear equations with sparse input matrix. Furthermore, we show how to compute verified error bounds for the solution of a system of nonlinear equations with sparse Jacobian. In all cases, the algorithms for sparse square linear systems of Part I and this Part II can be used.Publicationtype: Journal Article3 - Some of the metrics are blocked by yourconsent settings
Item type:Publication, Access status: without files , Verified error bounds for sparse systems part I: the splitting of a matrix into two factorsVerification methods provide mathematically correct error bounds for the solution of a numerical problem. That includes the proof of solvability of the problem and often uniqueness of the solution within the computed bounds. There are many verification methods for standard problems in numerical analysis, including linear and nonlinear systems of equations, matrix decompositions, eigenproblems, local and global optimization, and ordinary and partial differential equations. Many of those are included in INTLAB, the MATLAB/Octave toolbox for reliable computing. Despite several efforts, the verified solution of general sparse linear systems was an open problem. There are satisfactory algorithms for systems with symmetric positive definite input matrix. To that end error bounds for the solution of [Formula: see text] with general matrix [Formula: see text] could be computed using [Formula: see text], but that reduces the applicability in double precision to matrices with condition number up to [Formula: see text]. We give in this note an algorithm to compute entrywise error bounds for the solution of general real or complex sparse systems with condition number up to the limit [Formula: see text]. Our algorithm splits into three subalgorithms for symmetric positive definite, symmetric indefinite, and general input matrix [Formula: see text]. It is based on a mathematically correct lower bound on the smallest singular value [Formula: see text]. A key point is a factorization [Formula: see text] such that [Formula: see text] and [Formula: see text] have identical sets of singular values with the smallest one close to [Formula: see text]. A mathematically correct lower bound on [Formula: see text] is then computed using [Formula: see text]. Numerical evidence suggests that bounds for the solution of a linear system are computed for condition numbers up to [Formula: see text], and that often the bounds for all entries are close to maximally accurate, i.e., the bounds differ by few bits. We present a new a priori error bound based on Perron–Frobenius theory on the residual of a floating-point Cholesky decomposition which improves on existing ones by some two orders of magnitude. Our three subalgorithms benefit from this new bound. Based on the results in this Part I, an alternative approach will be presented in Part II of this note. Those methods are somewhat simpler, but often slower. However, they seem to be more robust, i.e., produce verified inclusions where the methods of Part I fail. Both approaches for square linear systems will be used in Part II of this note to compute verified error bounds for the solution of least squares problems and for underdetermined linear systems. Moreover, inclusions of the solution of general real or complex systems of nonlinear equations with sparse Jacobi matrix are computed by transforming the problem into a linear system with point matrix and interval right-hand side.Publicationtype: Journal Article19 - Some of the metrics are blocked by yourconsent settings
Item type:Publication, Access status: without files , A block householder-based algorithm for the QR decomposition of hierarchical matricesHierarchical matrices are dense but data-sparse matrices that use low-rank factorizations of suitable submatrices to reduce the storage and computational cost to linear-polylogarithmic complexity. In this paper, we propose a new approach to efficiently compute QR factorizations in the hierarchical matrix format based on block Householder transformations. To prevent unnecessarily high ranks in the resulting factors and to increase speed and accuracy, the algorithm meticulously tracks for which intermediate results low-rank factorizations are available. We also use a special storage scheme for the block Householder reflector to further reduce computational and storage costs. Numerical tests for two- and three-dimensional Laplacian boundary element matrices, different radial basis function kernel matrices, and matrices of typical hierarchical matrix structures but filled with random entries illustrate the performance of the new algorithm in comparison to some other QR algorithms for hierarchical matrices from the literature.Publicationtype: Journal Article64 - Some of the metrics are blocked by yourconsent settings
Item type:Publication, Access status: without files , Verified error bounds for matrix decompositions(2024); In this note we consider common matrix factorizations such as LU decomposition of a square and rectangular matrix, Cholesky and QR decomposition, singular value decomposition for square and rectangular matrices, and eigen-, Schur, and Takagi decomposition. We first note that well-conditioned factors tend to be sensitive to perturbations of the input matrix, while ill-conditioned factors tend to be insensitive. It seems that this behavior has not been recognized in numerical analysis. We develop a formula for the relation between the condition number of the factor and its sensitivity with respect to input perturbations, and give reasons for that. Our main focus is to describe verification methods for the factors of the mentioned decompositions. That means proving existence of the factorization together with rigorous entrywise error bounds for the factors. Our goal is to develop algorithms requiring \scrO(Pp2) operations for an m \times n matrix with P := max(m,n) and p := min(m,n). Moreover, bounds of high quality are aimed for, often not far from maximal accuracy. A main tool to achieve that is accurate dot products based on error-free transformations. Since preconditioning based on approximate inverses is used, our methods are restricted to full matrices.Publicationtype: Journal Article16 - Some of the metrics are blocked by yourconsent settings
Item type:Publication, Access status: without files , Verified error bounds for all eigenvalues and eigenvectors of a matrixA verification method is presented to compute error bounds for all eigenvectors and eigenvalues, including clustered and/or multiple ones of a general, real, or complex matrix. In case of a narrow cluster, error bounds for an invariant subspace are computed because computation of a single eigenvector may be ill-posed. Computer algebra and verification methods have in common that the computed results are correct with mathematical certainty. Unlike a computer algebra method, a verification method may fail in the sense that only partial or no inclusions at all are computed. That may happen for very ill conditioned problems being too sensitive for the arithmetical precision in use. That cannot happen for computer algebra methods which are “never-failing” because potentially infinite precision is used. In turn, however, that may slow down computer algebra methods significantly and may impose limitations on the problem size. In contrast, verification methods solely use floating-point operations so that their computing time and treatable problem size is of the order of that of purely numerical algorithms. For our problem it is proved that the union of the eigenvalue bounds contains the whole spectrum of the matrix, and bounds for corresponding invariant subspaces are computed. The computational complexity to compute inclusions of all eigenpairs of an n × nmatrix is O(n3).Publicationtype: Journal Article22 - Some of the metrics are blocked by yourconsent settings
Item type:Publication, Access status: without files , A linear relation approach to port-hamiltonian differential-algebraic equationsWe consider linear port-Hamiltonian differential-algebraic equations. Inspired by the geometric approach of van der Schaft and Maschke [System Control Lett., 121 (2018), pp. 31-37] and the linear algebraic approach of Mehl, Mehrmann, and Wojtylak [SIAM J. Matrix Anal. Appl., 39 (2018), pp. 1489-1519], we present another view by using the theory of linear relations. We show that this allows us to elaborate the differences and mutualities of the geometric and linear algebraic views, and we introduce a class of DAEs which comprises these two approaches. We further study the properties of matrix pencils arising from our approach via linear relations.Publicationtype: Journal Article60 - Some of the metrics are blocked by yourconsent settings
Item type:Publication, Access status: without files , Factorization, symmetrization, and truncated transformation of radial basis function-GA stabilized Gaussian radial basis functions(2019)Radial basis function (RBF) interpolation is a powerful, meshfree tool for function approximation but its direct implementation might suffer from severe ill-conditioning as the shape parameter decreases. We build upon RBF-GA, a stable approach for Gaussian RBFs, and derive its representation in terms of a matrix factorization which can then be generalized to a symmetrized version. This symmetrized version requires fewer function evaluations, yields new insight into the flat limit case, and, combined with diagonal scaling, allows a further reduction of the condition number. We also propose a truncated version of the basis transformation. We conclude with numerical tests to illustrate the performance of all introduced interpolation methods.Publicationtype: Journal Article214 - Some of the metrics are blocked by yourconsent settings
Item type:Publication, Access status: without files , The componentwise structured and unstructured backward errors can be arbitrarily far apartGiven a linear system Ax = b and some vector x, the backward error characterizes the smallest relative perturbation of A, b such that × is a solution of the perturbed system. If the input matrix has some structure such as being symmetric or Toeplitz, perturbations may be restricted to perturbations within the same class of structured matrices. For normwise perturbations, the symmetric and the general backward errors are equal, and the question about the relation between the symmetric and general componentwise backward errors arises. In this note we show for a number of common structures in numerical analysis that for componentwise perturbations the structured backward error can be equal to 1, whilst the unstructured backward error is arbitrarily small. Structures cover symmetric, persymmetric, skewsymmetric, Toeplitz, symmetric Toeplitz, Hankel, persymmetric Hankel, and circulant matrices. This is true although the normwise condition number ||A<sup>-1</sup>|| ||A|| is close to 1.Publicationtype: Journal Article60 - Some of the metrics are blocked by yourconsent settings
Item type:Publication, Access status: without files , Improved backward error bounds for LU and Cholesky factorizations(Society for Industrial & Applied Mathematics (SIAM), 2014-05-27); Assuming standard floating-point arithmetic (in base β, precision p) and barring underflow and overflow, classical rounding error analysis of the LU or Cholesky factorization of an n×n matrix A provides backward error bounds of the form |ΔA| < γn|L̂||Û| or |ΔA| < γn+1|R̂T||R̂|. Here, L̂, Û, and R̂ denote the computed factors, and γn is the usual fraction nu/(1-nu) = nu+O(u2) with u the unit roundoff. Similarly, when solving an n×n triangular system Tx = b by substitution, the computed solution x̂ satisfies (T + ΔT)x̂ = b with |ΔT| < γn|T|. All these error bounds contain quadratic terms in u and limit n to satisfy either nu < 1 or (n + 1)u < 1. We show in this paper that the constants γn and γn+1 can be replaced by nu and (n + 1)u, respectively, and that the restrictions on n can be removed. To get these new bounds the main ingredient is a general framework for bounding expressions of the form |ρ - s|, where s is the exact sum of a floating-point number and n-1 real numbers and where ρ is a real number approximating the computed sum ŝ. By instantiating this framework with suitable values of ρ, we obtain improved versions of the well-known Lemma 8.4 from [N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, Philadelphia, 2002] (used for analyzing triangular system solving and LU factorization) and of its Cholesky variant. All our results hold for rounding to nearest with any tie-breaking strategy and whatever the order of summation.Publicationtype: Journal Article55 - Some of the metrics are blocked by yourconsent settings
Item type:Publication, Access status: without files , Improved error bounds for inner products in floating-point arithmetic(Society for Industrial & Applied Mathematics (SIAM), 2013-04-04); Given two floating-point vectors x, y of dimension n and assuming rounding to nearest, we show that if no underflow or overflow occurs, any evaluation order for an inner product returns a floating-point number r̂ such that |r̂ - xT y| ≤ nu|x|T|y| with u the unit roundoff. This result, which holds for any radix and with no restriction on n, can be seen as a generalization of a similar bound given in [S. M. Rump, BIT, 52 (2012), pp. 201-220] for recursive summation in radix 2, namely, |r̂ - x T e| ≤ (n - 1)u|x|T e with e = [1, 1, ... , 1] T. As a direct consequence, the error bound for the floating-point approximation Ĉ of classical matrix multiplication with inner dimension n simplifies to |Ĉ - AB| ≤ nu|A||B|. © 2013 Society for Industrial and Applied Mathematics.Publicationtype: Journal Article50