Optimized look-ahead recurrences for adjacent rows in the Padé table

A new look-ahead algorithm for recursively computing Padé approximants is introduced. It generates a subsequence of the Padé approximants on two adjacent rows (defined by fixed numerator degree) of the Padé table. Its two basic versions reduce to the classical Levinson and Schur algorithms if no look-ahead is required. The new algorithm can be viewed as a combination of the look-ahead sawtooth and the look-ahead Levinson and Schur algorithms that we proposed before, but now the look-ahead step size is minimal (as in the sawtooth version) and the computational costs are as low as in the least expensive competing algorithms (including our look-ahead Levinson and Schur algorithms). The underlying recurrences link well-conditioned basic pairs,i.e., pairs of sufficiently different neighboring Padé forms.The algorithm can be used to solve Toeplitz systems of equationsTx = b. In this application it comes in several versions: anO(N ²) Levinson-type form, anO(N ²) Schur-type form, and a superfastO(N log² N) Schur-type version. As an option of the first two versions, the corresponding block LDU decompositions ofT ^–1 orT, respectively, can be found.     

Fast Algorithms with Preprocessing for Matrix-Vector Multiplication Problems

In this paper the problem of complexity of multiplication of a matrix with a vector is studied for Toeplitz, Hankel, Vandermonde, and Cauchy matrices and for matrices connected with them (i.e., for transpose, inverse, and transpose to inverse matrices). The proposed algorithms have complexities of at most O(n log²n) flops and in a number of cases they improve the known estimates. In these algorithms, in a separate preprocessing phase, are singled out all the actions on the preparation of a given matrix which aimed at the reduction of the complexity of the second stage of computations directly connected with multiplication by an arbitrary vector. Effective algorithms for computing the Vandermonde determinant and the determination of a Cauchy matrix are given.     

Fast numerical solution of fredholm integral equations with stationary kernels

A fast recursive matrix method for the numerical solution of Fredholm integral equations with stationary kernels is derived. IfN denotes the number of nodal points, the complexity of the algorithm isO(N ²), which should be compared toO(N ³) for conventional algorithms for solving such problems. The method is related to fast algorithms for inverting Toeplitz matrices.Applications to equations of the first and second kind as well as miscellaneous problems are discussed and illustrated with numerical examples. These show that the theoretical improvement in efficiency is indeed obtained, and that no problems with numerical stability or accuracy are encountered.     

Displacement ranks of matrices and linear equations

It takes of the order of N³ operations to solve a set of N linear equations in N unknowns or to invert the corresponding coefficient matrix. When the underlying physical problem has some time- or shift-invariance properties, the coefficient matrix is of Toeplitz (or difference or convolution) type and it is known that it can be inverted with O(N²) operations. However non-Toeplitz matrices often arise even in problems with some underlying time-invariance, e.g., as inverses or products or sums of products of possibly rectangular Toeplitz matrices. These non-Toeplitz matrices should be invertible with a complexity between O(N²) and O(N³). In this paper we provide some content for this feeling by introducing the concept of displacement ranks, which serve as a measure of how ‘close’ to Toeplitz a given matrix is.     

Block band Toeplitz preconditioners derived from generating function approximations: analysis and applications

We are concerned with the study and the design of optimal preconditioners for ill-conditioned Toeplitz systems that arise from a priori known real-valued nonnegative generating functions f(x,y) having roots of even multiplicities. Our preconditioned matrix is constructed by using a trigonometric polynomial θ(x,y) obtained from Fourier/kernel approximations or from the use of a proper interpolation scheme. Both of the above techniques produce a trigonometric polynomial θ(x,y) which approximates the generating function f(x,y), and hence the preconditioned matrix is forced to have clustered spectrum. As θ(x,y) is chosen to be a trigonometric polynomial, the preconditioner is a block band Toeplitz matrix with Toeplitz blocks, and therefore its inversion does not increase the total complexity of the PCG method. Preconditioning by block Toeplitz matrices has been treated in the literature in several papers. We compare our method with their results and we show the efficiency of our proposal through various numerical experiments.This research was co-funded by the European Union in the framework of the program “Pythagoras I” of the “Operational Program for Education and Initial Vocational Training” of the 3rd Community Support Framework of the Hellenic Ministry of Education, funded by national sources (25%) and by the European Social Fund - ESF (75%). The work of the second and of the third author was partially supported by MIUR (Italian Ministry of University and Research), grant number 2004015437.     

Band plus algebra preconditioners for two-level Toeplitz systems

In this paper we are interested in the fast and efficient solution of nm×nm symmetric positive definite ill-conditioned Block Toeplitz with Toeplitz Blocks (BTTB) systems of the form T _nm (f)x=b, where the generating function f is a priori known. The preconditioner that we propose and analyze is an extension of the one proposed in (D. Noutsos and P. Vassalos, Comput. Math. Appl., 56 (2008), pp. 1255–1270) and it arises as a product of a Block band Toeplitz matrix and matrices that may belong to any trigonometric matrix algebra. The underlying idea of the proposed scheme is to embody the well known advantages characterizing each component of the product when used alone. As a result we obtain spectral equivalence and a weak clustering of the eigenvalues of the preconditioned matrix around unity, ensuring the convergence of the Preconditioned Conjugate Gradient (PCG) method with a number of iterations independent of the partial dimensions. Finally, we compare our method with techniques already employed in the literature. A wide range of numerical experiments confirms the effectiveness of the proposed procedure and the adherence to the theoretical analysis.     

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.     

Asymptotically fast solution of toeplitz and related systems of linear equations

We present an inversion algorithm for the solution of a generic N X N Toeplitz system of linear equations with computational complexity O(Nlog²N) and storage requirements O(N). The algorithm relies upon the known structure of Toeplitz matrices and their inverses and achieves speed through a doubling method. All the results are derived and stated in terms of the recent concept of displacement rank, and this is used to extend the scope of the algorithm to include a wider class of matrices than just Toeplitz and also to include block Toeplitz matrices.     

Matrix decomposition algorithms for the <Emphasis Type="Italic">C</Emphasis><Superscript>0</Superscript>-quadratic finite element Galerkin method

Explicit expressions for the eigensystems of one-dimensional finite element Galerkin (FEG) matrices based on C ⁰ piecewise quadratic polynomials are determined. These eigensystems are then used in the formulation of fast direct methods, matrix decomposition algorithms (MDAs), for the solution of the FEG equations arising from the discretization of Poisson’s equation on the unit square subject to several standard boundary conditions. The MDAs employ fast Fourier transforms and require O(N ²log N) operations on an N×N uniform partition. Numerical results are presented to demonstrate the efficacy of these algorithms.     

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.     

Fast wavelet transforms and numerical algorithms I

A class of algorithms is introduced for the rapid numerical application of a class of linear operators to arbitrary vectors. Previously published schemes of this type utilize detailed analytical information about the operators being applied and are specific to extremely narrow classes of matrices. In contrast, the methods presented here are based on the recently developed theory of wavelets and are applicable to all Calderon-Zygmund and pseudo-differential operators. The algorithms of this paper require order O(N) or O(N log N) operations to apply an N × N matrix to a vector (depending on the particular operator and the version of the algorithm being used), and our numerical experiments indicate that many previously intractable problems become manageable with the techniques presented here.     

A fast method to block-diagonalize a Hankel matrix

In this paper, we consider an approximate block diagonalization algorithm of an n×n real Hankel matrix in which the successive transformation matrices are upper triangular Toeplitz matrices, and propose a new fast approach to compute the factorization in O(n ²) operations. This method consists on using the revised Bini method (Lin et al., Theor Comp Sci 315: 511–523, 2004). To motivate our approach, we also propose an approximate factorization variant of the customary fast method based on Schur complementation adapted to the n×n real Hankel matrix. All algorithms have been implemented in Matlab and numerical results are included to illustrate the effectiveness of our approach.     

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.     

A continuation method for solving symmetric Toeplitz systems

A fast algorithm is proposed for solving symmetric Toeplitz systems. This algorithm continuously transforms the identity matrix into the inverse of a given Toeplitz matrix T. The memory requirements for the algorithm are O(n), and its complexity is O(log κ(T)nlogn), where (T) is the condition number of T. Numerical results are presented that confirm the efficiency of the proposed algorithm.     

On neighbouring matrices with quadratic elementary divisors

Summary Algorithms are presented which compute theQR factorization of an order-n Toeplitz matrix inO(n ²) operations. The first algorithm computes onlyR explicitly, and the second computes bothQ andR. The algorithms are derived from a well-known procedure for performing the rank-1 update ofQR factors, using the shift-invariance property of the Toeplitz matrix. The algorithms can be used to solve the Toeplitz least-squares problem, and can be modified to solve Toeplitz systems inO(n) space.     

Representation of real discrete Fourier transform in terms of a new set of functions based upon Möbius inversion

The trigonometric functions sin(2n/N) and cos(2n/N) are transformed into a new set of basis functions using Möbius inversion of certain types of series. The new basis functions are number theoretic series. They are used to represent the real discrete Fourier transform (RDFT) in terms of 2 matrices of factorization. The first matrix, with elements 1, -1 and 0 is obtained by replacing cos(2k/N) and sin(2k/N) by (k/N + 1/4) and (k/N), where (x) is the bipolar rectangular wave function. The second matrix is block-diagonal where each block is a circular correlation and consists of the new basis functions. Some applications of the new representation are discussed.     

A fast solver for linear systems with displacement structure

We describe a fast solver for linear systems with reconstructible Cauchy-like structure, which requires O(rn ²) floating point operations and O(rn) memory locations, where n is the size of the matrix and r its displacement rank. The solver is based on the application of the generalized Schur algorithm to a suitable augmented matrix, under some assumptions on the knots of the Cauchy-like matrix. It includes various pivoting strategies, already discussed in the literature, and a new algorithm, which only requires reconstructibility. We have developed a software package, written in Matlab and C-MEX, which provides a robust implementation of the above method. Our package also includes solvers for Toeplitz(+Hankel)-like and Vandermonde-like linear systems, as these structures can be reduced to Cauchy-like by fast and stable transforms. Numerical experiments demonstrate the effectiveness of the software.     

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.     

Hyperbolic Wavelets and Multiresolution in H^{2}(\mathbb{T})

In signal processing and system identification for H²(\BbbT)H^{2}(\Bbb{T}) and H²(\BbbD)H^{2}(\Bbb{D}) the traditional trigonometric bases and trigonometric Fourier transform are replaced by the more efficient rational orthogonal bases like the discrete Laguerre, Kautz and Malmquist-Takenaka systems and the associated transforms. These bases are constructed from rational Blaschke functions, which form a group with respect to function composition that is isomorphic to the Blaschke group, respectively to the hyperbolic matrix group. Consequently, the background theory uses tools from non-commutative harmonic analysis over groups and the generalization of Fourier transform uses concepts from the theory of the voice transform. The successful application of rational orthogonal bases needs a priori knowledge of the poles of the transfer function that may cause a drawback of the method. In this paper we give a set of poles and using them we will generate a multiresolution in H²(\BbbT)H^{2}(\Bbb{T}) and H²(\BbbD)H^{2}(\Bbb{D}). The construction is an analogy with the discrete affine wavelets, and in fact is the discretization of the continuous voice transform generated by a representation of the Blaschke group over the space H²(\BbbT)H^{2}(\Bbb{T}). The constructed discretization scheme gives opportunity of practical realization of hyperbolic wavelet representation of signals belonging to H²(\BbbT)H^{2}(\Bbb{T}) and H²(\BbbD)H^{2}(\Bbb{D}) if we can measure their values on a given set of points inside the unit circle or on the unit circle. Convergence properties of the hyperbolic wavelet representation will be studied.     

Factoring wavelet transforms into lifting steps

This article is essentially tutorial in nature. We show how any discrete wavelet transform or two band subband filtering with finite filters can be decomposed into a finite sequence of simple filtering steps, which we call lifting steps but that are also known as ladder structures. This decomposition corresponds to a factorization of the polyphase matrix of the wavelet or subband filters into elementary matrices. That such a factorization is possible is well-known to algebraists (and expressed by the formulaSL(n;R[z, z⁻¹])=E(n;R[z, z⁻¹])); it is also used in linear systems theory in the electrical engineering community. We present here a self-contained derivation, building the decomposition from basic principles such as the Euclidean algorithm, with a focus on applying it to wavelet filtering. This factorization provides an alternative for the lattice factorization, with the advantage that it can also be used in the biorthogonal, i.e., non-unitary case. Like the lattice factorization, the decomposition presented here asymptotically reduces the computational complexity of the transform by a factor two. It has other applications, such as the possibility of defining a wavelet-like transform that maps integers to integers. Research Tutorial Acknowledgements and Notes. Page 264.