Fast generalized cross validation using Krylov subspace methods

The task of fitting smoothing spline surfaces to meteorological data such as temperature or rainfall observations is computationally intensive. The generalized cross validation (GCV) smoothing algorithm, if implemented using direct matrix techniques, is O(n ³) computationally, and memory requirements are O(n ²). Thus, for data sets larger than a few hundred observations, the algorithm is prohibitively slow. The core of the algorithm consists of solving series of shifted linear systems, and iterative techniques have been used to lower the computational complexity and facilitate implementation on a variety of supercomputer architectures. For large data sets though, the execution time is still quite high. In this paper we describe a Lanczos based approach that avoids explicitly solving the linear systems and dramatically reduces the amount of time required to fit surfaces to sets of data.      

A fast Fourier-collocation method for second boundary integral equations

In this paper we develop a fast collocation method for second boundary integral equations by the trigonometric polynomials. We propose a convenient way to compress the dense matrix representation of a compact integral operator with a smooth kernel under the Fourier basis and the corresponding collocation functionals. The compression leads to a sparse matrix with only O(nlog²n) number of nonzero entries, where 2n+1 denotes the order of the matrix. Thus we develop a fast Fourier-collocation method. We prove that the fast Fourier-collocation method gives the optimal convergence order up to a logarithmic factor. Moreover, we design a fast scheme for solving the corresponding truncated linear system. We establish that this algorithm preserves the quasi-optimal convergence of the approximate solution with requiring a number of O(nlog³n) multiplications.     

Multiple point evaluation on combined tensor product supports

We consider the multiple point evaluation problem for an n-dimensional space of functions [???1,1[ ^d ?? spanned by d-variate basis functions that are the restrictions of simple (say linear) functions to tensor product domains. For arbitrary evaluation points this task is faced in the context of (semi-)Lagrangian schemes using adaptive sparse tensor approximation spaces for boundary value problems in moderately high dimensions. We devise a fast algorithm for performing m?≥?n point evaluations of a function in this space with computational cost O(mlog ^d n). We resort to nested segment tree data structures built in a preprocessing stage with an asymptotic effort of O(nlog ^d???1 n).     

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.     

Multivariate interpolation with increasingly flat radial basis functions of finite smoothness

In this paper, we consider multivariate interpolation with radial basis functions of finite smoothness. In particular, we show that interpolants by radial basis functions in ℝ ^d with finite smoothness of even order converge to a polyharmonic spline interpolant as the scale parameter of the radial basis functions goes to zero, i.e., the radial basis functions become increasingly flat.     

An O(n log n) unidirectional distributed algorithm for extrema finding in a circle

In this paper we present algorithms, which given a circular arrangement of n uniquely numbered processes, determine the maximum number in a distributive manner. We begin with a simple unidirectional algorithm, in which the number of messages passed is bounded by 2 n log n + O(n). By making several improvements to the simple algorithm, we obtain a unidirectional algorithm in which the number of messages passed is bounded by 1.5nlogn + O(n). These algorithms disprove Hirschberg and Sinclair's conjecture that O(n²) is a lower bound on the number of messages passed in undirectional algorithms for this problem. At the end of the paper we indicate how our methods can be used to improve an algorithm due to Peterson, to obtain a unidirectional algorithm using at most 1.356nlogn + O(n) messages. This is the best bound so far on the number of messages passed in both the bidirectional and unidirectional cases.     

Competitive clustering algorithms based on ultrametric properties

We propose in this paper two new competitive unsupervised clustering algorithms: the first algorithm deals with ultrametric data, it has a computational cost of O(n). The second algorithm has two strong features: it is fast and flexible on the processed data type as well as in terms of precision. The second algorithm has a computational cost, in the worst case, of O(n²), and in the average case, of O(n). These complexities are due to exploitation of ultrametric distance properties. In the first method, we use the order induced by an ultrametric in a given space to demonstrate how we can explore quickly data proximity. In the second method, we create an ultrametric space from a sample data, chosen uniformly at random, in order to obtain a global view of proximities in the data set according to the similarity criterion. Then, we use this proximity profile to cluster the global set. We present an example of our algorithms and compare their results with those of a classic clustering method.     

Comparison of Radial Basis Functions

A survey of algorithms for approximation of multivariate functions with radial basis function (RBF) splines is presented. Algorithms of interpolating, smoothing, selecting the smoothing parameter, and regression with splines are described in detail. These algorithms are based on the feature of conditional positive definiteness of the spline radial basis function. Several families of radial basis functions generated by means of conditionally completely monotone functions are considered. Recommendations for the selection of the spline basis and preparation of initial data for approximation with the help of the RBF spline are given.     

On the complexity of calculating factorials

It is shown that n! can be evaluated with time complexity O(log log nM (n log n)), where M(n) is the complexity of multiplying two n-digit numbers together. This is effected, in part, by writing n! in terms of its prime factors. In conjunction with a fast multiplication this yields an O(n(log n log log n)²) complexity algorithm for n!. This might be compared to computing n! by multiplying 1 times 2 times 3, etc., which is ω(n² log n) and also to computing n! by binary splitting which is O(log nM(n log n)).     

Adaptive Itô-Taylor algorithm can optimally approximate the Itô integrals of singular functions

We deal with numerical approximation of stochastic Itô integrals of singular functions. We first consider the regular case of integrands belonging to the Hölder class with parameters r and ?. We show that in this case the classical Itô-Taylor algorithm has the optimal error Θ(n^−(r+?)). In the singular case, we consider a class of piecewise regular functions that have continuous derivatives, except for a finite number of unknown singular points. We show that any nonadaptive algorithm cannot efficiently handle such a problem, even in the case of a single singularity. The error of such algorithm is no less than n^{−min{1/2,r+?}}. Therefore, we must turn to adaptive algorithms. We construct the adaptive Itô-Taylor algorithm that, in the case of at most one singularity, has the optimal error O(n^−(r+?)). The best speed of convergence, known for regular functions, is thus preserved. For multiple singularities, we show that any adaptive algorithm has the error Ω(n^{−min{1/2,r+?}}), and this bound is sharp.     

Fast solution of toeplitz systems of equations and computation of Padé approximants

We present two new algorithms, ADT and MDT, for solving order-n Toeplitz systems of linear equations Tz = b in time O(n log²n) and space O(n). The fastest algorithms previously known, such as Trench's algorithm, require time $\text{Ω(n}{}^2\text{)}$ and require that all principal submatrices of T be nonsingular. Our algorithm ADT requires only that T be nonsingular. Both our algorithms for Toeplitz systems are derived from algorithms for computing entries in the Padé table for a given power series. We prove that entries in the Padé table can be computed by the Extended Euclidean Algorithm. We describe an algorithm EMGCD (Extended Middle Greatest Common Divisor) which is faster than the algorithm HGCD of Aho, Hopcroft and Ullman, although both require time O(n log²n), and we generalize EMGCD to produce PRSDC (Polynomial Remainder Sequence Divide and Conquer) which produces any iterate in the PRS, not just the middle term, in time O(n log²n). Applying PRSDC to the polynomials U₀(x) = x²ⁿ⁺¹ and U₁(x) = a₀ + a₁x + … + a_2nx²ⁿ gives algorithm AD (Anti-Diagonal), which computes any (m, p) entry along the antidiagonal m + p = 2n of the Padé table for U₁ in time O(n log²n). Our other algorithm, MD (Main-Diagonal), computes any diagonal entry (n, n) in the Padé table for a normal power series, also in time O(n log²n). MD is related to Schönhage's fast continued fraction algorithm. A Toeplitz matrix T is naturally associated with U₁, and the (n, n) Padé approximation to U₁ gives the first column of T^?1. We show how a formula due to Trench can be used to compute the solution z of Tz = b in time O(n log n) from the first row and column of T^?1. Thus, the Padé table algorithms AD and MD give O(n log²n) Toeplitz algorithms ADT and MDT. Trench's formula breaks down in certain degenerate cases, but in such cases a companion formula, the discrete analog of the Christoffel-Darboux formula, is valid and may be used to compute z in time O(n log²n) via the fast computation (by algorithm AD) of at most four Padé approximants. We also apply our results to obtain new complexity bounds for the solution of banded Toeplitz systems and for BCH decoding via Berlekamp's algorithm.     

An improved algorithm for the rectangle enclosure problem

Given a set of n rectangles in the plane, with sides parallel to the coordinate axes, the rectangle enclosure problem consists of finding all q pairs of rectangles such that one rectangle of the pair encloses the other. In this note we present an alternative algorithm to the one by Vaishnavi and Wood; while both techniques have worst-case running time O(nlog²n + q), ours uses optimal storage O(n) rather than O(nlog²n) of Vaishnavi and Wood. Our algorithm works entirely in-place and uses very conventional data structures.     

An exact algorithm for MAX-CUT in sparse graphs

We study exact algorithms for the MAX-CUT problem. Introducing a new technique, we present an algorithmic scheme that computes a maximum cut in graphs with bounded maximum degree. Our algorithm runs in time O^*(2^(1-(2/Δ))n). We also describe a MAX-CUT algorithm for general graphs. Its time complexity is O^*(2^mn/(m+n)). Both algorithms use polynomial space.     

On some algorithmic problems regarding the hairpin completion

We consider a few algorithmic problems regarding the hairpin completion. The first problem we consider is the membership problem of the hairpin and iterated hairpin completion of a language. We propose an O(nf(n)) and O(n²f(n)) time recognition algorithm for the hairpin completion and iterated hairpin completion, respectively, of a language recognizable in O(f(n)) time. We show that the n factor is not needed in the case of hairpin completion of regular and context-free languages. The n² factor is not needed in the case of iterated hairpin completion of context-free languages, but it is reduced to n in the case of iterated hairpin completion of regular languages. We then define the hairpin completion distance between two words and present a cubic time algorithm for computing this distance. A linear time algorithm for deciding whether or not the hairpin completion distance with respect to a given word is connected is also discussed. Finally, we give a short list of open problems which appear attractive to us.     

Fast enumeration of words generated by Dyck grammars

The problem of enumerating and denumerating words generated by Dyck grammars arises in the work of compilers for high-level programming languages and a number of other applications. The present paper proposes an algorithm for the fast enumeration and denumeration of words of Dyck languages; the complexity of this algorithm per one symbol of enumerated words is O(log³ n log log n) bit operations, provided that the Schönhage-Strassen multiplication and division algorithmis used. The well-knownmethods applied earlier possess complexityO(n) per one symbol of enumerated words. The construction of the proposed algorithm is based on the Ryabko method.     

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.     

An energy-minimization framework for monotonic cubic spline interpolation

This paper describes the use of cubic splines for interpolating monotonic data sets. Interpolating cubic splines are popular for fitting data because they use low-order polynomials and have C² continuity, a property that permits them to satisfy a desirable smoothness constraint. Unfortunately, that same constraint often violates another desirable property: monotonicity. It is possible for a set of monotonically increasing (or decreasing) data points to yield a curve that is not monotonic, i.e., the spline may oscillate. In such cases, it is necessary to sacrifice some smoothness in order to preserve monotonicity.The goal of this work is to determine the smoothest possible curve that passes through its control points while simultaneously satisfying the monotonicity constraint. We first describe a set of conditions that form the basis of the monotonic cubic spline interpolation algorithm presented in this paper. The conditions are simplified and consolidated to yield a fast method for determining monotonicity. This result is applied within an energy minimization framework to yield linear and nonlinear optimization-based methods. We consider various energy measures for the optimization objective functions. Comparisons among the different techniques are given, and superior monotonic C² cubic spline interpolation results are presented. Extensions to shape preserving splines and data smoothing are described.     

Reduction Algorithm for the NPMLE for the Distribution Function of Bivariate Interval-Censored Data

This article considers computational aspects of the nonparametric maximum likelihood estimator (NPMLE) for the distribution function of bivariate interval-censored data. The computation of the NPMLE consists of a parameter reduction step and an optimization step. This article focuses on the reduction step and introduces two new reduction algorithms: the Tree algorithm and the HeightMap algorithm. The Tree algorithm is mentioned only briefly. The HeightMap algorithm is discussed in detail and also given in pseudo code. It is a fast and simple algorithm of time complexityO(n²). This is an order faster than the best known algorithm thus far by Bogaerts and Lesaffre. We compare the new algorithms to earlier algorithms in a simulation study, and demonstrate that the new algorithms are significantly faster. Finally, we discuss how the HeightMap algorithm can be generalized to d-dimensional data with d > 2. Such a multivariate version of the HeightMap algorithm has time complexity O(n^d).