Parallel defect control

How can small-scale parallelism best be exploited in the solution of nonstiff initial value problems? It is generally accepted that only modest gains inefficiency are possible, and it is often the case that “fast” parallel algorithms have quite crude error control and stepsize selection components. In this paper we consider the possibility of using parallelism to improvereliability andfunctionality rather than efficiency. We present an algorithm that can be used with any explicit Runge-Kutta formula. The basic idea is to take several smaller substeps in parallel with the main step. The substeps provide an interpolation facility that is essentially free, and the error control strategy can then be based on a defect (residual) sample. If the number of processors exceeds (p ? 1)/2, wherep is the order of the Runge-Kutta formula, then the interpolant and the error control scheme satisfy very strong reliability conditions. Further, for a given orderp, the asymptotically optimal values for the substep lengths are independent of the problem and formula and hence can be computed a priori. Theoretical comparisons between the parallel algorithm and optimal sequential algorithms at various orders are given. We also report on numerical tests of the reliability and efficiency of the new algorithm, and give some parallel timing statistics from a 4-processor machine.     

Factoring matrices with a tree-structured sparsity pattern

Let A be a matrix whose sparsity pattern is a tree with maximal degree d_max. We show that if the columns of A are ordered using minimum degree on |A|+|A|^∗, then factoring A using a sparse LU with partial pivoting algorithm generates only O(d_maxn) fill, requires only O(d_maxn) operations, and is much more stable than LU with partial pivoting on a general matrix. We also propose an even more efficient and just-as-stable algorithm called sibling-dominant pivoting. This algorithm is a strict partial pivoting algorithm that modifies the column preordering locally to minimize fill and work. It leads to only O(n) work and fill. More conventional column pre-ordering methods that are based (usually implicitly) on the sparsity pattern of |A|^∗|A| are not as efficient as the approaches that we propose in this paper.     

Multi-dimensional option pricing using radial basis functions and the generalized Fourier transform

We show that the generalized Fourier transform can be used for reducing the computational cost and memory requirements of radial basis function methods for multi-dimensional option pricing. We derive a general algorithm, including a transformation of the Black–Scholes equation into the heat equation, that can be used in any number of dimensions. Numerical experiments in two and three dimensions show that the gain is substantial even for small problem sizes. Furthermore, the gain increases with the number of dimensions.     

A communication-less parallel algorithm for tridiagonal Toeplitz systems

Diagonally dominant tridiagonal Toeplitz systems of linear equations arise in many application areas and have been well studied in the past. Modern interest in numerical linear algebra is often focusing on solving classic problems in parallel. In McNally [Fast parallel algorithms for tri-diagonal symmetric Toeplitz systems, MCS Thesis, University of New Brunswick, Saint John, 1999], an m processor Split & Correct algorithm was presented for approximating the solution to a symmetric tridiagonal Toeplitz linear system of equations. Nemani [Perturbation methods for circulant-banded systems and their parallel implementation, Ph.D. Thesis, University of New Brunswick, Saint John, 2001] and McNally (2003) adapted the works of Rojo [A new method for solving symmetric circulant tri-diagonal system of linear equations, Comput. Math. Appl. 20 (1990) 61–67], Yan and Chung [A fast algorithm for solving special tri-diagonal systems, Computing 52 (1994) 203–211] and McNally et al. [A split-correct parallel algorithm for solving tri-diagonal symmetric Toeplitz systems, Internat. J. Comput. Math. 75 (2000) 303–313] to the non-symmetric case. In this paper we present relevant background from these methods and then introduce an m processor scalable communication-less approximation algorithm for solving a diagonally dominant tridiagonal Toeplitz system of linear equations.     

Optimal adaptive solution of initial-value problems with unknown singularities

The optimal solution of initial-value problems in ODEs is well studied for smooth right-hand side functions. Much less is known about the optimality of algorithms for singular problems. In this paper, we study the (worst case) solution of scalar problems with a right-hand side function having r   continuous bounded derivatives in R$\mathbf{R}$, except for an unknown singular point. We establish the minimal worst case error for such problems (which depends on r similarly as in the smooth case), and define optimal adaptive algorithms. The crucial point is locating an unknown singularity of the solution by properly adapting the grid. We also study lower bounds on the error of an algorithm for classes of singular problems. In the case of a single singularity with nonadaptive information, or in the case of two or more singularities, the error of any algorithm is shown to be independent of r.     

Point estimation of simultaneous methods for solving polynomial equations: A survey (II)

The construction of computationally verifiable initial conditions which provide both the guaranteed and fast convergence of the numerical root-finding algorithm is one of the most important problems in solving nonlinear equations. Smale's “point estimation theory” from 1981 was a great advance in this topic; it treats convergence conditions and the domain of convergence in solving an equation f(z)=0$f(z)=0$ using only the information of f   at the initial point _z0$z_0$. The study of a general problem of the construction of initial conditions of practical interest providing guaranteed convergence is very difficult, even in the case of algebraic polynomials. In the light of Smale's point estimation theory, an efficient approach based on some results concerning localization of polynomial zeros and convergent sequences is applied in this paper to iterative methods for the simultaneous determination of simple zeros of polynomials. We state new, improved initial conditions which provide the guaranteed convergence of frequently used simultaneous methods for solving algebraic equations: Ehrlich–Aberth's method, Ehrlich–Aberth's method with Newton's correction, Börsch-Supan's method with Weierstrass’ correction and Halley-like (or Wang–Zheng) method. The introduced concept offers not only a clear insight into the convergence analysis of sequences generated by the considered methods, but also explicitly gives their order of convergence. The stated initial conditions are of significant practical importance since they are computationally verifiable; they depend only on the coefficients of a given polynomial, its degree n and initial approximations to polynomial zeros.     

Computing iterative roots of polygonal functions

Based on the iterative root theory for monotone functions, an algorithm for computing polygonal iterative roots of increasing polygonal functions was given by J. Kobza. In this paper we not only give an algorithm for roots of decreasing polygonal functions but also generalize Kobza's results to the general n. Furthermore, we extend our algorithms for polygonal PM functions, a class of non-monotonic functions.     

B-spline collocation for solution of two-point boundary value problems

A numerical method based on B-spline is developed to solve the general nonlinear two-point boundary value problems up to order 6. The standard formulation of sextic spline for the solution of boundary value problems leads to non-optimal approximations. In order to derive higher orders of accuracy, high order perturbations of the problem are generated and applied to construct the numerical algorithm. The error analysis and convergence properties of the method are studied via Green’s function approach. O(h⁶) global error estimates are obtained for numerical solution of these classes of problems. Numerical results are given to illustrate the efficiency of the proposed method. Results of numerical experiments verify the theoretical behavior of the orders of convergence.     

Lattices on simplicial partitions which are not simply connected

In this paper, (d+1)-pencil lattices on simplicial partitions in R^d, which are not simply connected, are studied. It is shown, how the fact that a partition is not simply connected can be used to increase the flexibility of a lattice. A local modification algorithm is developed also to deal with slight partition topology changes that may appear afterwards a lattice has already been constructed.     

A generalized unitary Hessenberg matrix

This paper describes the use of a generalized isometric Arnoldi algorithm to reduce a unitary matrix, via unitary similarity, to a product of elementary reflectors and permutations. The computation is analogous to the reduction of a unitary matrix to a unitary Hessenberg matrix using the isometric Arnoldi algorithm. In the case in which A is a shift matrix, the reduction provides a novel recurrence for the factor R in the QR factorization of a Toeplitz-like matrix.     

Particular rules for the Θ-algorithm

The -algorithm is an extrapolation algorithm which can be very useful in accelerating some slowly convergent sequences. Like the other acceleration algorithms, the -algorithm is quite sensitive to the propagation of rounding errors due to cancellation in the difference between two almost equal quantities.In order to (partially) avoid this drawback, particular rules are given. They have to be used, instead of the usual rules of the algorithm, when two adjacent quantities in a column are nearly equal. Numerical examples show that these rules can improve the numerical stability of the algorithm in some cases while, in other cases, the improvement is non-existent.     

B-splines and discretization in an inverse problem for Poisson processes

Recent results on quasi-maximum likelihood histogram sieve estimators in inverse problems for Poisson processes are generalized to B-spline sieves. The impact of discretization effects on strong L² consistency and convergence rates are studied in detail. In particular, a “rates saturation effect”, caused by discretization, is demonstrated. Finite-sample implementation is proposed and tested in a Monte Carlo experiment with the Wicksell problem, which shows a superior performance of the new approach, when compared to other methods commonly used in that context. The proposed algorithm can also be used in cases with only approximately known folding kernel.     

A full-Newton step non-interior continuation algorithm for a class of complementarity problems

In this paper, we investigate a class of nonlinear complementarity problems arising from the discretization of the free boundary problem, which was recently studied by Sun and Zeng [Z. Sun, J. Zeng, A monotone semismooth Newton type method for a class of complementarity problems, J. Comput. Appl. Math. 235 (5) (2011) 1261–1274]. We propose a new non-interior continuation algorithm for solving this class of problems, where the full-Newton step is used in each iteration. We show that the algorithm is globally convergent, where the iteration sequence of the variable converges monotonically. We also prove that the algorithm is globally linearly and locally superlinearly convergent without any additional assumption, and locally quadratically convergent under suitable assumptions. The preliminary numerical results demonstrate the effectiveness of the proposed algorithm.     

D-optimal designs embedded in Hadamard matrices and their effect on the pivot patterns

In this paper we develop a new approach for detecting if specific D-optimal designs exist embedded in Sylvester-Hadamard matrices. Specifically, we investigate the existence of the D-optimal designs of orders 5, 6, 7 and 8. The problem is motivated to explaining why specific values appear as pivot elements when Gaussian elimination with complete pivoting is applied to Hadamard matrices. Using this method and a complete search algorithm we explain, for the first time, the appearance of concrete pivot values for equivalence classes of Hadamard matrices of orders n = 12, 16 and 20.     

Fast Fourier-Galerkin methods for solving singular boundary integral equations: Numerical integration and precondition

We develop a fast fully discrete Fourier-Galerkin method for solving a class of singular boundary integral equations. We prove that the number of multiplications used in generating the compressed matrix is O(nlog³n), and the solution of the proposed method preserves the optimal convergence order O(n^−t), where n is the order of the Fourier basis functions used in the method and t denotes the degree of regularity of the exact solution. Moreover, we propose a preconditioning which ensures the numerical stability when solving the preconditioned linear system. Numerical examples are presented to confirm the theoretical estimates and to demonstrate the approximation accuracy and computational efficiency of the proposed algorithm.     

A reliable modification of the cross rule for rational Hermite interpolation

Claessens' cross rule [8] enables simple computation of the values of the rational interpolation table if the table is normal, i.e. if the denominators in the cross rule are non-zero. In the exceptional case of a vanishing denominator a singular block is detected having certain structural properties so that some values are known without further computations. Nevertheless there remain entries which cannot be determined using only the cross rule.In this note we introduce a simple recursive algorithm for computation of the values of neighbours of the singular block. This allows to compute entries in the rational interpolation table along antidiagonals even in the presence of singular blocks. Moreover, in the case of non-square singular blocks, we discuss a facility to monitor the stability.Dedicated to Professor G. Mühlbach on the occasion of his 50th birthday     

A variable step implicit block multistep method for solving first-order ODEs

A new four-point implicit block multistep method is developed for solving systems of first-order ordinary differential equations with variable step size. The method computes the numerical solution at four equally spaced points simultaneously. The stability of the proposed method is investigated. The Gauss-Seidel approach is used for the implementation of the proposed method in the PE(CE)^m mode. The method is presented in a simple form of Adams type and all coefficients are stored in the code in order to avoid the calculation of divided difference and integration coefficients. Numerical examples are given to illustrate the efficiency of the proposed method.     

On Khachiyan's algorithm for the computation of minimum-volume enclosing ellipsoids

Given A?{a¹,…,a^m}⊂R^d whose affine hull is R^d, we study the problems of computing an approximate rounding of the convex hull of A and an approximation to the minimum-volume enclosing ellipsoid of A. In the case of centrally symmetric sets, we first establish that Khachiyan's barycentric coordinate descent (BCD) method is exactly the polar of the deepest cut ellipsoid method using two-sided symmetric cuts. This observation gives further insight into the efficient implementation of the BCD method. We then propose a variant algorithm which computes an approximate rounding of the convex hull of A, and which can also be used to compute an approximation to the minimum-volume enclosing ellipsoid of A. Our algorithm is a modification of the algorithm of Kumar and Y?ld?r?m, which combines Khachiyan's BCD method with a simple initialization scheme to achieve a slightly improved polynomial complexity result, and which returns a small “core set.” We establish that our algorithm computes an approximate solution to the dual optimization formulation of the minimum-volume enclosing ellipsoid problem that satisfies a more complete set of approximate optimality conditions than either of the two previous algorithms. Furthermore, this added benefit is achieved without any increase in the improved asymptotic complexity bound of the algorithm of Kumar and Y?ld?r?m or any increase in the bound on the size of the computed core set. In addition, the “dropping idea” used in our algorithm has the potential of computing smaller core sets in practice. We also discuss several possible variants of this dropping technique.     

A direct method to solve block banded block Toeplitz systems with non-banded Toeplitz blocks

A fast solution algorithm is proposed for solving block banded block Toeplitz systems with non-banded Toeplitz blocks. The algorithm constructs the circulant transformation of a given Toeplitz system and then by means of the Sherman-Morrison-Woodbury formula transforms its inverse to an inverse of the original matrix. The block circulant matrix with Toeplitz blocks is converted to a block diagonal matrix with Toeplitz blocks, and the resulting Toeplitz systems are solved by means of a fast Toeplitz solver.The computational complexity in the case one uses fast Toeplitz solvers is equal to ξ(m,n,k)=O(mn³)+O(k³n³) flops, there are m block rows and m block columns in the matrix, n is the order of blocks, 2k+1 is the bandwidth. The validity of the approach is illustrated by numerical experiments.