Fast and numerically stable algorithms for discrete Hartley transforms and applications to preconditioning

Abstract

Abstract. The discrete Hartley transforms (DHT) of types I -IV and the related matrix algebras are discussed. We prove that any of these DHTs of length N = 2 t can be factorized by means of a divide-and-conquer strategy into a product of sparse, orthogonal matrices where in this context sparse means at most two nonzero entries per row and column. The sparsity joint with orthogonality of the matrix factors is the key for proving that these new algorithms have low arithmetic costs equal to 5 2 N log 2 (N )+ O(N ) arithmetic operations and an excellent normwise numerical stability. Further, we consider the optimal Frobenius approximation of a given symmetric Toeplitz matrix generated by an integrable symbol in a Hartley matrix algebra. We give explicit formulas for computing these optimal approximations and discuss the related preconditioned conjugate gradient (PCG) iterations. By using the matrix approximation theory, we prove the strong clustering at unity of the preconditioned matrix sequences under the sole assumption of continuity and positivity of the generating function. The multilevel case is also briefly treated. Some numerical experiments concerning DHT preconditioning are included.

On Block Matrices Associated with Discrete Trigonometric Transforms and Their Use in the Theory of Wave Propagation

“…where ⊕ denotes the orthogonal sum of linear subspaces of C n×n with respect to the Frobenius inner product [16]. In particular, from (2.8), Diag(H X n ) are subspaces of T n + HA n , and hence every matrix of Diag(H X n ) may be expressed as a symmetric Toeplitz plus Hankel matrix.…”

“…where M is a Fourier, Hartley, Cosine, or Sine square matrix of order n. The respective DTT is named the Discrete Fourier, Hartley, Cosine, or Sine Transform (DFT, DHT, DCT, or DST) of length n. We collect the 4 types I, II, III, and IV of each one of the n × n Fourier, Hartley, Cosine, and Sine matrices, commonly used in the literature (see, e.g., [6,16,18]) by (j, k ∈ {0, . .…”

“…To formalize a unified situation, we introduce, in a similar way to [16], for a fixed unitary n × n complex matrix U, the term matrix algebra of U-diagonalizable n × n matrices:…”

“…5.1a and 5.1b depict the CPU time for the inversions with respect to mn of the present method and the LU decomposition for m × m block matrices with F Furthermore, we consider matrices with symmetric Toeplitz-plus-Hankel blocks, composed of random real entries, chosen from a normal distribution with mean zero, variance one, and standard deviation one. The Toeplitz and the Hankel counterparts of these blocks, belonging to H X n matrix algebras, are constructed by following the general guidelines of [16], modified into the present context in Table 5.1. III ndiagonalizable Toeplitz-plus-Hankel blocks of order n for n=300, 400, 500.…”

“…These unitary matrices offer simultaneous diagonalizations of specific matrix algebras, including circulants, skew-circulants, Toeplitz-plus-Hankel, tridiagonal [16]. These matrix algebras are unified by considering the algebra Diag(U ) of all U -diagonalizable matrices, for U a fixed unitary matrix.…”