Exact Arithmetic Research Articles

Language and Math: What If We Have Two Separate Naming Systems?

The role of language in numerical processing has traditionally been restricted to counting and exact arithmetic. Nevertheless, the impact that each of a bilinguals’ languages may have in core numerical representations has not been questioned until recently. What if the language in which math has been first acquired (LLmath) had a bigger impact in our math processing? Based on previous studies on language switching we hypothesize that balanced bilinguals would behave like unbalanced bilinguals when switching between the two codes for math. In order to address this question, we measured the brain activity with magneto encephalography (MEG) and source estimation analyses of 12 balanced Basque-Spanish speakers performing a task in which participants were unconscious of the switches between the two codes. The results show an asymmetric switch cost between the two codes for math, and that the brain areas responsible for these switches are similar to those thought to belong to a general task switching mechanism. This implies that the dominances for math and language could run separately from the general language dominance.

Sound and robust solid modeling via exact real arithmetic and continuity

Algorithms for solid modeling, i.e., Computer-Aided Design (CAD) and computer graphics, are often specified on real numbers and then implemented with finite-precision arithmetic, such as floating-point. The result is that these implementations do not soundly compute the results that are expected from their specifications. We present a new library, StoneWorks, that provides sound and robust solid modeling primitives. We implement StoneWorks in MarshallB, a pure functional programming language for exact real arithmetic in which types denote topological spaces and functions denote continuous maps, ensuring that all programs are sound and robust. We developed MarshallB as an extension of the Marshall language. We also define a new shape representation, compact representation ( K-rep ), that enables constructions such as Minkowski sum and analyses such as Hausdorff distance that are not possible with traditional representations. K-rep is a nondeterminism monad for describing all the points in a shape. With our library, language, and representation together, we show that short StoneWorks programs can specify and execute sound and robust solid modeling algorithms and tasks.

A New Projected Variant of the Deflated Block Conjugate Gradient Method

The deflated block conjugate gradient (D-BCG) method is an attractive approach for the solution of symmetric positive definite linear systems with multiple right-hand sides. However, the orthogonality between the block residual vectors and the deflation subspace is gradually lost along with the process of the underlying algorithm implementation, which usually causes the algorithm to be unstable or possibly have delayed convergence. In order to maintain such orthogonality to keep certain level, full reorthogonalization could be employed as a remedy, but the expense required is quite costly. In this paper, we present a new projected variant of the deflated block conjugate gradient (PD-BCG) method to mitigate the loss of this orthogonality, which is helpful to deal with the delay of convergence and thus further achieve the theoretically faster convergence rate of D-BCG. Meanwhile, the proposed PD-BCG method is shown to scarcely have any extra computational cost, while having the same theoretical properties as D-BCG in exact arithmetic. Additionally, an automated reorthogonalization strategy is introduced as an alternative choice for the PD-BCG method. Numerical experiments demonstrate that PD-BCG is more efficient than its counterparts especially when solving ill-conditioned linear systems or linear systems suffering from rank deficiency.

Exact Arithmetic of Zero and Infinity Based on the Conventional Division by Zero 0/0 = 1

The chief object of this work is to create an exact and consistent arithmetic of zero, denoted 0, and infinity (zero divisor), written as 1/0 and denoted ∞ , based on the conventional division by zero $$ \dfrac{{0}}{{0}}=1. $$ Manifold and undeniable applications of this arithmetic are given in this work in order to show its usefulness.

An extension of Mobius--Lie geometry with conformal ensembles of cycles and its implementation in a GiNaC library

We propose to consider ensembles of cycles (quadrics), which are interconnected through conformal-invariant geometric relations (e.g. ``to be orthogonal'', ``to be tangent'', etc.), as new objects in an extended M\"obius--Lie geometry. It was recently demonstrated in several related papers, that such ensembles of cycles naturally parameterize many other conformally-invariant families of objects, e.g. loxodromes or continued fractions. The paper describes a method, which reduces a collection of conformally in\-vari\-ant geometric relations to a system of linear equations, which may be accompanied by one fixed quadratic relation. To show its usefulness, the method is implemented as a {\CPP} library. It operates with numeric and symbolic data of cycles in spaces of arbitrary dimensionality and metrics with any signatures. Numeric calculations can be done in exact or approximate arithmetic. In the two- and three-dimensional cases illustrations and animations can be produced. An interactive {\Python} wrapper of the library is provided as well.

Block Modified Gram--Schmidt Algorithms and Their Analysis

New block modified Gram--Schmidt (BMGS) methods for the Q-R factorization of a full column rank matrix $X \in {{R}}^{m \times n}$, $m \geq n$, are considered. Such methods factor $X$ into $Q \in {{R}}^{m \times n}$ and an upper triangular $R \in {{R}}^{m \times n}$ such that $X = QR,$ where, in exact arithmetic, $Q$ is left orthogonal (i.e., $Q^T\;{Q}=I_n$). Gram--Schmidt-based algorithms play an important role in the implementation of Krylov space methods, such as GMRES, Arnoldi, and Lanczos. For these applications, the left orthogonal factor $Q$ is needed and the matrix is produced either one column at a time or one block of columns at a time. For block implementations of Krylov methods, a block of columns of $X$ is introduced at each step, and a new block of columns of $Q$ must then be produced. That is a task for which BMGS methods are ideally suited. However, for these Krylov methods to converge properly, the BMGS algorithms need to have numerical behavior that is similar to that of modified Gram--Schmidt. To design such BMGS algorithms, we build upon the block Householder representation of Schreiber and Van Loan [SIAM J. Sci. Stat. Comput., 10 (1989), pp. 53--57] and an observation by Charles Sheffield analyzed by Paige [SIAM J. Matrix Anal. Appl., 31 (2009), pp. 565--583] about the relationship between modified Gram--Schmidt and Householder Q-R factorization. Our new BMGS algorithms exploit the Sheffield framework so that they share a similar relationship to Householder Q-R and thus have error analysis properties similar to modified Gram--Schmidt. The last BMGS algorithm developed is based entirely upon matrix multiplications and the “tall, skinny” Q-R (TSQR) factorization---two operations that have been studied extensively for cache-based architectures and distributed architectures. It is shown that if the TSQR part of the BMGS algorithm satisfies error analysis properties connected to the Sheffield structure, then so does the entire BMGS algorithm. Thus new criteria for when a BMGS algorithm has error analysis properties similar to those of MGS are proposed.

Möbius--Lie Geometry and Its Extension

These lectures review the classical Moebius-Lie geometry and recent work on its extension. The latter considers ensembles of cycles (quadrics), which are interconnected through conformal-invariant geometric relations (e.g. "to be orthogonal", "to be tangent", etc.), as new objects in an extended Moebius--Lie geometry. It is shown on examples, that such ensembles of cycles naturally parameterise many other conformally-invariant families of objects, two examples---the Poincare extension and continued fractions are considered in detail. Further examples, e.g. loxodromes, wave fronts and integrable systems, are published elsewhere. The extended Moebius--Lie geometry is efficient due to a method, which reduces a collection of conformally invariant geometric relations to a system of linear equations, which may be accompanied by one fixed quadratic relation. The algorithmic nature of the method allows to implement it as a C++ library, which operates with numeric and symbolic data of cycles in spaces of arbitrary dimensionality and metrics with any signatures. Numeric calculations can be done in exact or approximate arithmetic. In the two- and three-dimensional cases illustrations and animations can be produced. An interactive Python wrapper of the library is provided as well.

Euclidean-Norm Error Bounds for SYMMLQ and CG

For positive definite and semidefinite consistent $Ax_\star=b$, we use the Gauss--Radau approach of Golub and Meurant (1997) to obtain an upper bound on the error $\|x_\star-x_k^L\|_2$ for SYMMLQ iterates, assuming exact arithmetic. Such a bound, computable in constant time per iteration, was not previously available. We show that the CG error $\|x_\star-x_k^C\|_2$ is always smaller and can also be bounded in constant time per iteration. Our approach is computationally cheaper than other bounds or estimates of the CG error in the literature. As with other approaches using Gauss--Radau quadrature, we require a positive lower bound on the smallest nonzero eigenvalue of $A$. For indefinite $A$, we obtain an estimate of $\|x_\star-x_k^L\|_2$. Numerical experiments demonstrate that our bounds are remarkably tight for SYMMLQ on positive definite systems and therefore provide reliable bounds for CG.

When integration sparsification fails: Banded Galerkin discretizations for Hermite functions, rational Chebyshev functions and sinh-mapped Fourier functions on an infinite domain, and Chebyshev methods for solutions with C[formula omitted] endpoint singularities

Chebyshev polynomial spectral methods are very accurate, but are plagued by the cost and ill-conditioning of dense discretization matrices. Modified schemes, collectively known as “integration sparsification”, have mollified these problems by discretizing the highest derivative as a diagonal matrix. Here, we examine five case studies where the highest derivative diagonalization fails. Nevertheless, we show that Galerkin discretizations do yield banded matrices that retain most of the advantages of “integration sparsification”. Symbolic computer algebra greatly extends the reach of spectral methods. When spectral methods are implemented using exact rational arithmetic, as is possible for small truncation N in Maple, Mathematica and their ilk, roundoff error is irrelevant, and sparsification failure is not worrisome. When the discretization contains a parameter L, symbolic algebra spectral methods return, as answer to an eigenproblem, not discrete numbers but rather a plane algebraic curve defined as the zero set of a bivariate polynomial P(λ,L); the optimal approximations to the eigenvalues λj are in the middle of the straight portions of the zero contours of P(λ;L) where the isolines are parallel to the L axis.

Fuzzy linear least squares for the identification of possibilistic regression models

The stochastic approach to regression, e.g. the least-squares method, is a commonly used technique in a wide range of disciplines as it offers an intuitive and well-founded solution to the problem of fitting a set of parameters to existing data. Specifically, in the case where the parameters of meta models are identified, the linear least-squares formulation receives wide appreciation among engineers due to its simple yet powerful formulation and solution. Usually, this approach provides a reasonable fit, depending largely on the quality of the model that is employed. However, it has two main drawbacks. From a statistical point of view, the premises justifying its use cannot always be assumed in good faith, and most importantly, the resulting model will almost certainly fail at correctly predicting the output without falling victim to the effect of over-fitting.The presented fuzzy regression algorithm originates from a set of different axioms. For the solution of the linear least-squares problem, it is able to provide the worst-case optimal parameter variation necessary to cover all of the crisp data points while encoding meaningful information in the membership function. This is accomplished by means of a sensible fuzzification procedure of the crisp data points and the application of an exact inverse fuzzy arithmetic for estimating the fuzzy-valued parameters of the model in question.

Differentiating through conjugate gradient

We show that, although the conjugate gradient (CG) algorithm has a singularity at the solution, it is possible to differentiate forward through the algorithm automatically by re-declaring all the variables as truncated Taylor series, the type of active variable widely used in automatic differentiation (AD) tools such as ADOL-C. If exact arithmetic is used, this approach gives a complete sequence of correct directional derivatives of the solution, to arbitrary order, in a single cycle of at most n iterations, where n is the number of dimensions. In the inexact case, the approach emphasizes the need for a means by which the programmer can communicate certain conditions involving derivative values directly to an AD tool.

The $QR$ Steps with Perfect Shifts

In this paper we revisit the problem of performing a $QR$ step on an unreduced Hessenberg matrix $H$ when we know an “exact” eigenvalue $\lambda_0$ of $H$. In exact arithmetic, this eigenvalue will...

Computing the Wave-Kernel Matrix Functions

We derive an algorithm for computing the wave-kernel functions $\cosh\sqrt{A}$ and $\mathrm{sinhc}\sqrt{A}$ for an arbitrary square matrix $A$, where $\mathrm{sinhc}z=\sinh(z)/z$. The algorithm is based on Pade approximation and the use of double angle formulas. We show that the backward error of any approximation to $\cosh\sqrt{A}$ can be explicitly expressed in terms of a hypergeometric function. To bound the backward error we derive and exploit a new bound for $\|A^k\|^{1/k}$ that is sharper than one previously obtained by Al-Mohy and Higham [SIAM J. Matrix Anal. Appl., 31 (2009), pp. 970--989]. The amount of scaling and the degree of the Pade approximant are chosen to minimize the computational cost subject to achieving backward stability for $\cosh\sqrt{A}$ in exact arithmetic. Numerical experiments show that the algorithm behaves in a forward stable manner in floating-point arithmetic and is superior in this respect to the general purpose Schur--Parlett algorithm applied to these functions.

Solving rank-constrained semidefinite programs in exact arithmetic

We consider the problem of minimizing a linear function over an affine section of the cone of positive semidefinite matrices, with the additional constraint that the feasible matrix has prescribed rank. When the rank constraint is active, this is a non-convex optimization problem, otherwise it is a semidefinite program. Both find numerous applications especially in systems control theory and combinatorial optimization, but even in more general contexts such as polynomial optimization or real algebra. While numerical algorithms exist for solving this problem, such as interior-point or Newton-like algorithms, in this paper we propose an approach based on symbolic computation. We design an exact algorithm for solving rank-constrained semidefinite programs, whose complexity is essentially quadratic on natural degree bounds associated to the given optimization problem: for subfamilies of the problem where the size of the feasible matrix, or the dimension of the affine section, is fixed, the algorithm is polynomial time. The algorithm works under assumptions on the input data: we prove that these assumptions are generically satisfied. We implement it in Maple and discuss practical experiments.

SPECTRA – a Maple library for solving linear matrix inequalities in exact arithmetic

This document describes our freely distributed Maple library spectra, for Semidefinite Programming solved Exactly with Computational Tools of Real Algebra. It solves linear matrix inequalities with symbolic computation in exact arithmetic and it is targeted to small-size, possibly degenerate problems for which symbolic infeasibility or feasibility certificates are required.

Exact duals and short certificates of infeasibility and weak infeasibility in conic linear programming

In conic linear programming -- in contrast to linear programming -- the Lagrange dual is not an exact dual: it may not attain its optimal value, or there may be a positive duality gap. The corresponding Farkas' lemma is also not exact (it does not always prove infeasibility). We describe exact duals, and certificates of infeasibility and weak infeasibility for conic LPs which are nearly as simple as the Lagrange dual, but do not rely on any constraint qualification. Some of our exact duals generalize the SDP duals of Ramana, and Klep and Schweighofer to the context of general conic LPs. Some of our infeasibility certificates generalize the row echelon form of a linear system of equations: they consist of a small, trivially infeasible subsystem obtained by elementary row operations. We prove analogous results for weakly infeasible systems. We obtain some fundamental geometric corollaries: an exact characterization of when the linear image of a closed convex cone is closed, and an exact characterization of nice cones. Our infeasibility certificates provide algorithms to generate {\em all} infeasible conic LPs over several important classes of cones; and {\em all} weakly infeasible SDPs in a natural class. Using these algorithms we generate a public domain library of infeasible and weakly infeasible SDPs. The status of our instances can be verified by inspection in exact arithmetic, but they turn out to be challenging for commercial and research codes.

Certifying solutions to overdetermined and singular polynomial systems over [formula omitted

This paper is concerned with certifying that a given point is near an exact root of an overdetermined or singular polynomial system with rational coefficients. The difficulty lies in the fact that consistency of overdetermined systems is not a continuous property. Our certification is based on hybrid symbolic-numeric methods to compute an exact rational univariate representation (RUR) of a component of the input system from approximate roots. For overdetermined polynomial systems with simple roots, we compute an initial RUR from approximate roots. The accuracy of the RUR is increased via Newton iterations until the exact RUR is found, which we certify using exact arithmetic. Since the RUR is well-constrained, we can use it to certify the given approximate roots using α-theory. To certify isolated singular roots, we use a determinantal form of the isosingular deflation, which adds new polynomials to the original system without introducing new variables. The resulting polynomial system is overdetermined, but the roots are now simple, thereby reducing the problem to the overdetermined case. We prove that our algorithms have complexity that are polynomial in the input plus the output size upon successful convergence, and we use worst case upper bounds for termination when our iteration does not converge to an exact RUR. Examples are included to demonstrate the approach.

A cooperative conjugate gradient method for linear systems permitting efficient multi-thread implementation

This paper revisits, in a multi-thread context, the so-called mul-tiparameter or block conjugate gradient (BCG) methods, first proposed as sequential algorithms by O'Leary and Brezinski, for the solution of the linear system Ax = b, for an n-dimensional symmetric positive definite matrix A. Instead of the scalar parameters of the classical CG algorithm, which minimizes a scalar functional at each iteration, multiple descent and conjugate directions are updated simultaneously. Implementation involves the use of multiple threads and the algorithm is referred to as cooperative CG (CCG) in order to emphasize that each thread now uses information that comes from the other threads. It is shown that for a sufficiently large matrix dimension n, the use of an optimal number of threads results in a worst case flop count of O(n 7/3) in exact arithmetic. Numerical experiments on a multicore, multi-thread computer , for synthetic and real matrices, illustrate the theoretical results.

Anisotropic geometry-conforming d-simplicial meshing via isometric embeddings

We develop a dimension-independent, Delaunay-based anisotropic mesh generation algorithm suitable for integration with adaptive numerical solvers. As such, the mesh produced by our algorithm conforms to an anisotropic metric prescribed by the solver as well as the domain geometry, given as a piecewise smooth complex. Motivated by the work of Lévy and Dassi [10-12,20], we use a discrete manifold embedding algorithm to transform the anisotropic problem to a uniform one. This work differs from previous approaches in several ways. First, the embedding algorithm is driven by a Riemannian metric field instead of the Gauss map, lending itself to general anisotropic mesh generation problems. Second we describe our method for computing restricted Voronoi diagrams in a dimension-independent manner which is used to compute constrained centroidal Voronoi tessellations. In particular, we compute restricted Voronoi simplices using exact arithmetic and use data structures based on convex polytope theory. Finally, since adaptive solvers require geometry-conforming meshes, we offer a Steiner vertex insertion algorithm for ensuring the extracted dual Delaunay triangulation is homeomorphic to the input geometries.The two major contributions of this paper are: a method for isometrically embedding arbitrary mesh-metric pairs in higher dimensional Euclidean spaces and a dimension-independent vertex insertion algorithm for producing geometry-conforming Delaunay meshes. The former is demonstrated on a two-dimensional anisotropic problem whereas the latter is demonstrated on both 3d and 4d problems.

Speeding up Exact Real Arithmetic on Fast Binary Cauchy Sequences by Using Memoization Based on Quantized Precision

Speeding up Exact Real Arithmetic on Fast Binary Cauchy Sequences by Using Memoization Based on Quantized Precision