Discrete fourier transform computation using prime Ramanujan numbers

Ramanujan numbers were introduced in [2] to implement discrete fourier transform (DFT) without using any multiplication operation. Ramanujan numbers are related to π and integers which are powers of 2. If the transform sizeN, is a Ramanujan number, then the computational complexity of the algorithms used for computing isO(N ²) addition and shift operations, and no multiplications. In these algorithms, the transform can be computed sequentially with a single adder inO(N ²) addition times. Parallel implementation of the algorithm can be executed inO(N) addition times, withO(N) number of adders. Some of these Ramanujan numbers of order-2 are related to the Biblical and Babylonian values of π [1]. In this paper, we analytically obtain upper bounds on the degree of approximation in the computation of DFT if JV is a prime Ramanujan number.     

FFTs on the Rotation Group

We discuss an implementation of an efficient algorithm for the numerical computation of Fourier transforms of bandlimited functions defined on the rotation group SO(3). The implementation is freely available on the web. The algorithm described herein uses O(B ⁴) operations to compute the Fourier coefficients of a function whose Fourier expansion uses only (the O(B ³)) spherical harmonics of degree at most B. This compares very favorably with the direct O(B ⁶) algorithm derived from a basic quadrature rule on O(B ³) sample points. The efficient Fourier transform also makes possible the efficient calculation of convolution over SO(3) which has been used as the analytic engine for some new approaches to searching 3D databases (Funkhouser et al., ACM Trans. Graph. 83–105, [2003]; Kazhdan et al., Eurographics Symposium in Geometry Processing, pp. 167–175, [2003]). Our implementation is based on the “Separation of Variables” technique (see, e.g., Maslen and Rockmore, Proceedings of the DIMACS Workshop on Groups and Computation, pp. 183–237, [1997]). In conjunction with techniques developed for the efficient computation of orthogonal polynomial expansions (Driscoll et al., SIAM J. Comput. 26(4):1066–1099, [1997]), our fast SO(3) algorithm can be improved to give an algorithm of complexity O(B ³log ² B), but at a cost in numerical reliability. Numerical and empirical results are presented establishing the empirical stability of the basic algorithm. Examples of applications are presented as well. First author was supported by NSF ITR award; second author was supported by NSF Grant 0219717 and the Santa Fe Institute.     

Using Random Quasi-Monte-Carlo Within Particle Filters,With Application to Financial Time Series

This article presents a new particle filter algorithm which uses random quasi-Monte-Carlo to propagate particles. The filter can be used generally, but here it is shown that for one-dimensional state-space models, if the number of particles is N, then the rate of convergence of this algorithm is N^?1. This compares favorably with the N^?1/2 convergence rate of standard particle filters. The computational complexity of the new filter is quadratic in the number of particles, as opposed to the linear computational complexity of standard methods. I demonstrate the new filter on two important financial time series models, an ARCH model and a stochastic volatility model. Simulation studies show that for fixed CPU time, the new filter can be orders of magnitude more accurate than existing particle filters. The new filter is particularly efficient at estimating smooth functions of the states, where empirical rates of convergence are N^?3/2; and for performing smoothing, where both the new and existing filters have the same computational complexity.     

On Computing the Subset Graph of a Collection of Sets

Let a given collection of sets have size N measured by the sum of the cardinalities. Yellin and Jutla presented an algorithm which constructed the partial order induced by the subset relation (a “subset graph”) in a claimed O(N²/log N) operations over a dictionary ADT, and exhibited a collection whose subset graph had Θ(N²/log² N) edges. This paper first establishes a matching upper bound on the number of edges in a subset graph. It also presents a finer analysis of the algorithm, which confirms the claimed upper bound and shows it to be tight. A simple implementation requiring O(1) bit-parallel operations per ADT operation is presented, along with a variant of the algorithm with an implementation requiring O(N²/log N) RAM operations.     

A capacitance matrix method for Dirichlet problem on polygon region

Summary An efficient algorithm for the solution of linear equations arising in a finite element method for the Dirichlet problem is given. The cost of the algorithm is proportional toN ²log₂ N (N=1/h) where the cost of solving the capacitance matrix equations isNlog₂ N on regular grids andN ^3/2log₂ N on irregular ones.     

Massively parallel tabu search for the quadratic assignment problem

A new heuristic algorithm to perform tabu search on the Quadratic Assignment Problem (QAP) is developed. A massively parallel implementation of the algorithm on the Connection Machine CM-2 is provided. The implementation usesn ² processors, wheren is the size of the problem. The elements of the algorithm, calledPar_tabu, include dynamically changing tabu list sizes, aspiration criterion and long term memory. A new intensification strategy based on intermediate term memory is proposed and shown to be promising especially while solving large QAPs. The combination of all these elements gives a very efficient heuristic for the QAP: the best known or improved solutions are obtained in a significantly smaller number of iterations than in other comparative studies. Combined with the implementation on CM-2, this approach provides suboptimal solutions to QAPs of bigger dimensions in reasonable time.     

Remarks on a displacement-rank inversion method for Toeplitz systems

Comments are made regarding the implementation of a Toeplitz-matrix inversion algorithm described by Bitmead and Anderson in [1]. We show that although the algorithm is asymptotically efficient with O(N(logN)²) operations, it requires a 10⁶×10⁶ matrix to break even with the class of algorithms whose operation count is of the order of O(N²) (as found in [4]).     

A fast transform for spherical harmonics

Spherical harmonics arise on the sphere S² in the same way that the (Fourier) exponential functions {e^ik}k arise on the circle. Spherical harmonic series have many of the same wonderful properties as Fourier series, but have lacked one important thing: a numerically stable fast transform analogous to the Fast Fourier Transform (FFT). Without a fast transform, evaluating (or expanding in) spherical harmonic series on the computer is slow—for large computations probibitively slow. This paper provides a fast transform.For a grid ofO(N²) points on the sphere, a direct calculation has computational complexityO(N⁴), but a simple separation of variables and FFT reduce it toO(N³) time. Here we present algorithms with timesO(N^5/2 log N) andO(N²(log N)²). The problem quickly reduces to the fast application of matrices of associated Legendre functions of certain orders. The essential insight is that although these matrices are dense and oscillatory, locally they can be represented efficiently in trigonometric series.     

Numerical inversion of the Radon transform

Summary The problem of inverting the Radon transform, i.e. the reconstruction of a function inR ² from its line integrals arises e.g. in computerized tomography and in nondestructive testing. In the present paper the least squares method with piecewise constant trial functions is investigated. An error estimate is derived. An implementation using the fast Fourier transform is described and numerical results are reported.     

Spectral analysis of the finite Hankel transform and circular prolate spheroidal wave functions

In this paper, we develop two practical methods for the computation of the eigenvalues as well as the eigenfunctions of the finite Hankel transform operator. These different eigenfunctions are called circular prolate spheroidal wave functions (CPSWFs). This work is motivated by the potential applications of the CPSWFs as well as the development of practical methods for computing their values. Also, in this work, we should prove that the CPSWFs form an orthonormal basis of the space of Hankel band-limited functions, an orthogonal basis of L²([0,1]) and an orthonormal system of L²([0,+∞[). Our computation of the CPSWFs and their associated eigenvalues is done by the use of two different methods. The first method is based on a suitable matrix representation of the finite Hankel transform operator. The second method is based on the use of an efficient quadrature method based on a special family of orthogonal polynomials. Also, we give two Maple programs that implement the previous two methods. Finally, we present some numerical results that illustrate the results of this work.     

CoSaMP: Iterative signal recovery from incomplete and inaccurate samples

Compressive sampling offers a new paradigm for acquiring signals that are compressible with respect to an orthonormal basis. The major algorithmic challenge in compressive sampling is to approximate a compressible signal from noisy samples. This paper describes a new iterative recovery algorithm called CoSaMP that delivers the same guarantees as the best optimization-based approaches. Moreover, this algorithm offers rigorous bounds on computational cost and storage. It is likely to be extremely efficient for practical problems because it requires only matrix–vector multiplies with the sampling matrix. For compressible signals, the running time is just O(Nlog²N), where N is the length of the signal.     

Efficient dynamic programming implementations of Newton's method for unconstrained optimal control problems

Naive implementations of Newton's method for unconstrainedN-stage discrete-time optimal control problems with Bolza objective functions tend to increase in cost likeN ³ asN increases. However, if the inherent recursive structure of the Bolza problem is properly exploited, the cost of computing a Newton step will increase only linearly withN. The efficient Newton implementation scheme proposed here is similar to Mayne's DDP (differential dynamic programming) method but produces the Newton step exactly, even when the dynamical equations are nonlinear. The proposed scheme is also related to a Riccati treatment of the linear, two-point boundary-value problems that characterize optimal solutions. For discrete-time problems, the dynamic programming approach and the Riccati substitution differ in an interesting way; however, these differences essentially vanish in the continuous-time limit.This work was supported by the National Science Foundation, Grant No. DMS-85-03746.     

Vekua theory for the Helmholtz operator

Vekua operators map harmonic functions defined on domain in \mathbb R²{\mathbb R^{2}} to solutions of elliptic partial differential equations on the same domain and vice versa. In this paper, following the original work of I. Vekua (Ilja Vekua (1907–1977), Soviet-Georgian mathematician), we define Vekua operators in the case of the Helmholtz equation in a completely explicit fashion, in any space dimension N ≥ 2. We prove (i) that they actually transform harmonic functions and Helmholtz solutions into each other; (ii) that they are inverse to each other; and (iii) that they are continuous in any Sobolev norm in star-shaped Lipschitz domains. Finally, we define and compute the generalized harmonic polynomials as the Vekua transforms of harmonic polynomials. These results are instrumental in proving approximation estimates for solutions of the Helmholtz equation in spaces of circular, spherical, and plane waves.     

A quadrature formula for the Hankel transform

We present in this paper a quadrature formula for a certain Fourier-Bessel transform and, closely related to this, for the Hankel transform of order >–1. Such formulas originate in the context of a Galerkin-type projection of the weightedL ₂(–, ; ) space ( is the weight function mentioned below) used to get a discrete representation of a certain physical problem in Quantum Mechanics. The generalized Hermitee polynomialsH ₀ (x),H ₁ (x),..., with weight function (x), are used as the basis on which such a projection takes place. It is shown that theN-dimensional vectors representing certain projected functions as well as the entries of theN×N matrix representing the kernel of that Fourier-Bessel transform, approach the exact functional values at the zeros of theNth generalized Hermitee polynomial whenN.These properties lead to propose this matrix as a finite representation of the kernel of the Fourier-Bessel transform involved in this problem and theN zeros of the generalized Hermitee polynomialH _N (x) as abscissas to yield certain quadrature formulae for this integral and for the related Hankel transform. The error function produced by this algorithm is estimated at theN nodes and its is shown to be of a smaller order than 1/N. This error estimate is valid for piecewise continuous functions satisfying certain integral conditions involving their absolute values. The algorithm is presented with some numerical examples.     

About strongly polynomial time algorithms for quadratic optimization over submodular constraints

We present new strongly polynomial algorithms for special cases of convex separable quadratic minimization over submodular constraints. The main results are: an O(NM log(N ²/M)) algorithm for the problemNetwork defined on a network onM arcs andN nodes; an O(n logn) algorithm for thetree problem onn variables; an O(n logn) algorithm for theNested problem, and a linear time algorithm for theGeneralized Upper Bound problem. These algorithms are the best known so far for these problems. The status of the general problem and open questions are presented as well.This research has been supported in part by ONR grant N00014-91-J-1241.Corresponding author.     

Efficient computation of Fourier transforms on compact groups

This article genralizes the fast Fourier transform algorithm to the computation of Fourier transforms on compact Lie groups. The basic technique uses factorization of group elements and Gel'fand-Tsetlin bases to simplify the computations, and may be extended to treat the computation of Fourier transforms of finitely supported distributions on the group. Similar transforms may be defined on homogeneous spaces; in that case we show how special function properties of spherical functions lead to more efficient algorithms. These results may all be viewed as generalizations of the fast Fourier transform algorithms on the circle, and of recent results about Fourier transforms on finite groups. Acknowledgements and Notes. This paper was written while the author was supported by the Max-Planck-Institut für Mathematik, Bonn, Germany.     

On the k-error linear complexity of sequences with period 2pn over GF(q)

In this paper, we first optimize the structure of the Wei–Xiao–Chen algorithm for the linear complexity of sequences over GF(q) with period N =  2p ⁿ , where p and q are odd primes, and q is a primitive root modulo p ². The second, an union cost is proposed, so that an efficient algorithm for computing the k-error linear complexity of a sequence with period 2p ⁿ over GF(q) is derived, where p and q are odd primes, and q is a primitive root modulo p ². The third, we give a validity of the proposed algorithm, and also prove that there exists an error sequence e ^N , where the Hamming weight of e ^N is not greater than k, such that the linear complexity of (s + e) ^N reaches the k-error linear complexity c. We also present a numerical example to illustrate the algorithm. Finally, we present the minimum value k for which the k-error linear complexity is strictly less than the linear complexity.     

New fast divide‐and‐conquer algorithms for the symmetric tridiagonal eigenvalue problem

In this paper, two accelerated divide‐and‐conquer (ADC) algorithms are proposed for the symmetric tridiagonal eigenvalue problem, which cost O(N²r) flops in the worst case, where N is the dimension of the matrix and r is a modest number depending on the distribution of eigenvalues. Both of these algorithms use hierarchically semiseparable (HSS) matrices to approximate some intermediate eigenvector matrices, which are Cauchy‐like matrices and are off‐diagonally low‐rank. The difference of these two versions lies in using different HSS construction algorithms, one (denoted by ADC1) uses a structured low‐rank approximation method and the other (ADC2) uses a randomized HSS construction algorithm. For the ADC2 algorithm, a method is proposed to estimate the off‐diagonal rank. Numerous experiments have been carried out to show their stability and efficiency. These algorithms are implemented in parallel in a shared memory environment, and some parallel implementation details are included. Comparing the ADCs with highly optimized multithreaded libraries such as Intel MKL, we find that ADCs could be more than six times faster for some large matrices with few deflations. Copyright © 2016 John Wiley & Sons, Ltd.     

Implementing an efficient minimum capacity cut algorithm

In this paper, we present an efficient implementation of theO(mn + n ² logn) time algorithm originally proposed by Nagamochi and Ibaraki (1992) for computing the minimum capacity cut of an undirected network. To enhance computation, various ideas are added so that it can contract as many edges as possible in each iteration. To evaluate the performance of the resulting implementation, we conducted extensive computational experiments, and compared the results with that of Padberg and Rinaldi's algorithm (1990), which is currently known as one of the practically fastest programs for this problem. The results indicate that our program is considerably faster than Padberg and Rinaldi's program, and its running time is not significantly affected by the types of the networks being solved.Corresponding author.