Fast and Numerically Stable Approximate Solution of Trummer’s Problem

Trummer’s problem is the problem of multiplication of an n × n Cauchy matrix C by a vector. It serves as the basis for the solution of several problems in scientific computing and engineering [1]. The straightforward algorithm solves Trummer’s problem in O(n2) flops. The fast algorithm solves the problem in O(nlog2n) flops [2] but has poor numerical stability. The algorithm we discuss here in this paper is the celebrated multipoint algorithm [3] which has been studied by Pan et al. The algorithm approximates the solution in O(nlogn) flops in terms of n but its cost estimate depends on the bound of the approximation error and also depends on the correlation between the entries of the pair of n-dimensional vectors defining the input matrix C.


Introduction
Computations with dense structured matrices have many applications in sciences, communications and engineering.The structure enables dramatic acceleration of the computations and major decrease in memory space but sometimes leads to numerical stability problems.The best well-known classes of structured matrices are Toeplitz, Hankel, Cauchy and Vandermonde matrices.
The computations with such matrices are widely applied in the areas of algebraic coding, control, signal processing, solution of partial differential equations and algebraic computing.For example, Toeplitz matrices arise in some major signal processing computations and the problem of multiplying Vandermonde matrix by a vector is equivalent to polynomial evaluation, whereas solving a Vandermonde system is equivalent to polynomial interpolation.Moreover, Cauchy matrices appear in the study of integral equations and conformal mappings.The complexity of computations with n × n dense structured matrices dramatically decreases in compari-son with the general n × n matrices, that is, from the order of n 2 words of storage space and n α arithmetic op- erations (ops) with 2.37 3 α < ≤ in the best algorithms, to ( ) O n words of storage space and ( ) log O n n ops (see Table 1 below for more details).

Some Basic Definitions
Throughout this paper, we use the following notations;  denotes the set of nonnegative integers, +  is the set of positive real numbers, Z + denotes the set of positive integers, and R denotes the set of real numbers.for every pair of its entries , i j t and 1, 1 for every pair of its entries , i j h and 1, 1 i j h − + .
Definition 2.2.For a given vector , the matrix ( ) Definition 2.3.Given two vectors s and t such that i j s t ≠ for all i and j, the n × n matrix ( ) s t is a Cauchy (generalized Hilbert) where ( ) For more details regarding the four classes of structured matrices, see Table 2 below.
It is quite easy to verify that TJ and JT are Hankel matrices if T is a Toeplitz matrix, and HJ and JH are Toeplitz matrices if H is a Hankel matrix where J is the following reflection matrix, 0 1

The Displacement Operators of Dense Structured Matrices
The concept of displacement operators and displacement rank which was introduced by T. Kailath, S. Y. Kung, and M. Morf in 1979 and studied by Pan, Bini, and other authors is one of the powerful tools for studying and dealing with matrices that have structure.The displacement rank approach, when it was initially introduced, was intended for more restricted use [4], namely, to measure how "close" to Toeplitz a given matrix is.Then the idea turned out to be even more powerful, thus it was developed, generalized and extended to other structured matrices.In this section, we consider the most general and modern interpretation of the displacement of a matrix.The main idea is, for a given structured matrix A, we need to find an operator L that transforms the matrix into a low rank matrix ( ) L A such that one can easily recover A from its image ( ) L A and operate with low rank matrices instead.Such operators that shift and scale the entries of the structured matrices turn out to be appropriate tools for introducing and defining the matrices of Toeplitz-like, Hankel-like, Vandermonde-like, and Cauchy-like types [5].Definition 3.1.For any fixed field  such as the complex field  and a fixed pair { } , M N of operator matrices, we define the linear displacement operators :

A MA AN M MM AN M A M AN M A A MA AN MAN N AN MAN A N A MAN A N
The operator matrices that we will be using are the matrices Z f , T f Z , and ( ) We may use the operator matrices Z 1 and Z 0 in the case of Toeplitz matrices, Z 1

The Correlation of Structured Matrices to Polynomials
The product ( ) represents the vector of the values of the polynomial ( ) is the vector of the n th roots of unity, 2π 2π 2π e cos sin , where 1, and and multipoint evaluation turns into discrete Fourier transform which takes only ( ) O n n ops and allows numerically stable implementation according to [6].If we express ( ) V x in Equation (1) via Cauchy matrices we will get ( ) ( ) ( ) ( ) Note that the numerical stability is very important in approximation algorithm for multipoint polynomial evaluation.It relies on expressing ( ) V x in terms of Cauchy matrices as in Equation (2).Clearly in Equation ( 2), the product ( ) V x by a vector has been reduced to ones with Cauchy, which brings us to Trummer's problem, that is, the problem of multiplication of Cauchy matrix by a vector.Its solution by multipoint algorithm ( [6], pp.261-262) leads to multipoint polynomial evaluation based on Equation (2) which is fast in terms of ops and numerically stable as it was proved by Pan.
We may vary the vector x by linearly mapping it to the vector a b = + y x e where we can take

New Transformation of Cauchy Matrices
As we mentioned earlier, Trummer's problem is the problem of multiplication of an n × n Cauchy matrix C by a vector which is the basis for the solution of many important problems of scientific computing and engineering.The straightforward algorithm solves Trummer's problem in ( )

2
O n flops.The fast algorithm solves the problem in ( ) n flops but has poor numerical stability.
The algorithm we presenting in this paper approximates the solution in ( ) log O n n flops in terms of n.However, its cost estimate depends on the bound of the approximation error and on the correlation between the entries of the pair of n-dimensional vectors defining the input matrix C.This algorithm is numerically stable as we will see throughout this section and the next section.
The main goal in this paper is to enrich the power of the multipoint algorithm by introducing and proving some new expressions for Cauchy matrix via other Cauchy matrices [7], which we may vary by changing one of their basis vectors.Under a certain choice of such a vector, the solution of Trummer's problem will be simplified; thus, the power of the multipoint algorithm can be developed as we will see in the next section.
Therefore, we will achieve our goal by using a simple transformation of the useful basic formula of [8], and the resulting expressions for C will give us further algorithmic opportunities.Definition 5.1.For a pair of n-dimensional vectors ( ) ( ) (See [8]) The main idea of the transformation of the basic vectors defining the problem is taken from [9], where this idea was used for multipoint polynomial evaluation and interpolation.( ) , where ( ) is the vector filled with the values one and δ is a scalar parameter.Proof: due to Lemma 5.1.

3) Proof of Equation (7):
From Equation ( 4 we have also Start with Equation (4) which is ( )  c d  d c  d c  d  c  c   d c  d  c  c   d c  d  c This implies that ( ) ( ) ( ) ( ) ( ) ( ) which is Equation (7).

4) Proof of Equation (8):
From Equation ( 4 and take the transpose of both side of the last equation to get: Substitute Equation (10) and the matrix equation ( )

Approximate Stable Solutions of Trummer's Problems
The algorithm we are studying and presenting in this section depends on the multipoint algorithm which approximates the solution of Trummer's problem in ( ) O n ops in terms of n, and works efficiently for a large classes of input vectors but sometimes has problems with some of the input vectors, especially if the ratio of the input vectors is close to 1.
Recall, the power series expansion: The basis for the algorithm is the following expressions: Once again for large M the expression The product of Cauchy matrix by a vector is ( ) . This is just Trummer's problem.If we simplify this expression, we get the following approximations: For any n × n Cauchy matrix ( ) , C s t , the approximation requires ( ) O nM ops for all i and it is numerically stable.
If either of the ratios i j s t or j i t s is small, then the approximate error will be small for large M.However, there will be a problem whenever one of the ratios is close to 1, in this case, the error will be large.

Discussions and Conclusions
Recall the following two formulas from Section 5: Equations ( 13) and ( 14) are Vandermonde-free and Hankel-free, but they enable us to transform the basis vectors s and t for ( ) , C s t into the two pairs of basis vectors s, q and q, t for any choice of the vector ( ) . Then Trummer's problem is reduced to the evaluation of the diagonal matrices ( ) for ( ) , s t , ( ) , q t , ( ) , s q , ( ) , q s , ( ) , t q and/or ( ) , t s and also reduced to recursive multiplication of the above matrices and ( ) , C q t and ( ) . We compute the coefficients by simply pairwise multiply the linear factors and the matrices ( ) , D u q and ( ) D q can be immediately evaluated in ( ) log O n n flops.In addition, any polynomial ( ) p x of degree n can be evaluated at the scaled n th roots of 1 in ( ) log O n n ops by means of Fast Fourier Transform (FFT).Trummer's problem is the multiplication of ( ) , C q t by a vector or ( ) , C s q by a vector.Its solution can be simplified under appropriate choice of the vector q.One way to do it is to restrict q to the above choice in Equation (15).Even with this particular choice, yet the scalar a allows faster convergence of the power series of the Multipole Algorithm presented in Section 7.This can be extended to the Equations ( 5) and (7).On the other hand, one can linearly map the vector q into a b = + y q e, where ( ) , and 0 a ≠ and b are any scalars.In addition, the computations of the diagonal matrices will be simplified if our choice of the vector q is the scaled nth root of unity. ,δ is a scalar parameter.Hence, ( ) ( )

L
A of the operator L is called the displacement of the matrix A. The operators of Sylvester and Stein types can be transformed easily into one another if at least one of the two associated operator matrices is non-singular.The following theorem explains this fact.

and T 1 Z
in the case of Han-kel matrices, Z 0 and ( ) D v in the case of Vandermonde matrices, and ( ) D x and( ) D y in the case of Cauchy matrices.However, there are other choices of operator matrices that can transform these matrices to low rank.
≠ and b are any scalars.

≠
denote the asso- ciated n × n Cauchy, Vandermonde, and triangular Hankel matrices, respectively.For a vector ( ) a Cauchy degenerate matrix ( ) C a has the diagonal entries zeros and the ( ) of n × n diagonal matrices, defined by the vectors a and b.Theorem 5.1.

Definition 5 . 2 .
Trummer's problem is the problem of computing the vector ( ) , C a b v for three given vectors ( ) for all pairs i, j.Trummer's degenerate problem is the problem of computing the vector ( ) C a v for two given vectors solution of Trummer's degenerate problem can be reduced to Trummer's problem due to the next simple result.
This is done by taking the inverse of Equation (3) and replacing the vectors c, d by b, d.Then substitute the equation

(
Equation (3) for ( ) , C c d into the following matrix identity ( ) ( ) ( ) ( ) the vector d by b and obtain (that is, use the identity , recursively, the computed products.The computation is numerically stable and uses (ops), but it is not numerically stable; therefore the fast and numerically stable approximation techniques of[10] can be used instead.If we choose any vector n th roots of unity for a scalar a and

Remark 8 . 1 .
Trummer's problem frequently arises for Cauchy degenerate matrices that are defined as follows:

Table 1 .
Parameter and flops count.
Matrices A size n × n Number of parameters for A

Table 2 .
General definition of the four classes of structured matrices.