An extra-component method for evaluating fast matrix-vector multiplication with special functions

Abstract

In calculating integral or discrete transforms, use has been made of fast algorithms for multiplying vectors by matrices whose elements are specified as values of special (Chebyshev, Legendre, Jacobi, Laguerre, etc.) functions. To achieve higher efficiency, a convenient general approach for calculating matrix-vector products for some class of problems is proposed. A series of fast simple-structure algorithms developed under this approach can be efficiently implemented with software based on modern microprocessors. Computational experiments show that the method has a precomputation complexity of $O(N^2 \log N)$ and an execution complexity of $O(N \log N)$ in general. This agrees with the analytical estimates obtained for trigonometric functions. In addition, for trigonometric functions, the number of arithmetic operations at the preliminary stage can be decreased to $O(N \log N)$. The results of computational experiments with the algorithms show that these procedures can decrease the calculation time by several orders of magnitude compared with a conventional direct method of matrix-vector multiplication.

Appendix: An efficient precomputation procedure for trigonometric transforms

In the case of calculation of nonuniform trigonometric transformations, the precomputation costs of Algorithms 1 and 2 can be reduced from $O(N^2\log N)$ to $O(N \log N +Nq)$ arithmetic operations, where q is determined by the required computation accuracy. As an example, let us consider an efficient procedure for calculating a row of the compressed matrix (7).

First, according to the inverse convolution theorem, for any $\mathbf {X},\mathbf {Y}\in \mathbb {C}^{N}$ we have

$$\begin{aligned} \mathcal {F}\left( \mathbf {X} \odot \mathbf {Y} \right) = \mathcal {F}\mathbf {X} *\mathcal {F}\mathbf {Y}=\tilde{\mathbf {X}} *\tilde{\mathbf {Y}}, \end{aligned}$$

where $\odot$ denotes the component-wise product and

$$\begin{aligned} \left( \tilde{\mathbf {X}} *\tilde{\mathbf {Y}}\right) _n\triangleq \frac{1}{\sqrt{N}}\sum \limits _{m=0}^{N-1}\tilde{x}_m\tilde{y}_{(n-m)\text {mod}\ N}, \quad n=0,\ldots ,N-1. \end{aligned}$$

(A.1)

Figure 2 shows that the Kaiser window can be approximated by a small number of coefficients of the Fourier series denoted here as $\tilde{x}_m$. Since in the discrete case for the Kaiser window no efficient formula for calculating individual Fourier components is known, the calculation of several components of the spectrum $\tilde{x}_m$ using FFT will require $O(N\log N)$ operations. This operation is performed once for all columns or rows of the transformation matrix.

Second, Fig. 3b shows that there is no need to calculate all elements of a compressed matrix row by formula (A.1). It is sufficient to calculate about 24 convolutions of the form (A.1) for each row to ensure an accuracy of transformation of the order $10^{-15}$.

Third, depending on the trigonometric transformation type given by the matrix A, the individual DFT components for the rows of the matrix A can be calculated by one of the following formulas:

$$\begin{aligned}&\tilde{y}^{\mathrm {cos}}_k(\theta )=\frac{1}{\sqrt{N}}\sum\nolimits_{j=0}^{N-1}\cos \left( j\theta \right) \omega _N^{jk}\nonumber \\= & {} \frac{1}{\sqrt{N}}\frac{\omega _N^k\left( (\left( \cos \left( N \theta \right) -1)\cos \left( \theta \right) +\sin \left( \theta \right) \sin \left( N \theta \right) \right) -(\cos \left( N \theta \right) -1)\omega _N^{k}\right) }{\omega _N^{2k}-2 \omega _N^{k} \cos \left( \theta \right) +1}, \end{aligned}$$

(A.2)

$$\begin{aligned} \tilde{y}^{\mathrm {sin}}_k(\theta )= & {} \frac{1}{\sqrt{N}}\sum \nolimits _{j=0}^{N-1}\sin \left( j\theta \right) \omega _N^{jk}\nonumber \\= & {} \frac{1}{\sqrt{N}}\frac{\omega _N^{k}\left( \left( \sin \left( \theta \right) (\cos \left( N \theta \right) -1)-\cos \left( \theta \right) \sin \left( N\theta \right) \right) +\sin \left( N \theta \right) \omega _N^{k}\right) }{-\omega _N^{2k}+2 \omega _N^{k} \cos \left( \theta \right) -1},\end{aligned}$$

(A.3)

$$\begin{aligned} \tilde{y}^{\mathrm {exp}}_k(\theta )= & {} \frac{1}{\sqrt{N}}\sum \nolimits _{j=0}^{N-1}{\mathrm e}^{\mathrm {i} j\theta }\omega _N^{jk}=\frac{1}{\sqrt{N}}\frac{{\mathrm e}^{\mathrm {i} N \theta }-1}{{\mathrm e}^{\mathrm {i}\theta }\omega _N^{k}-1}, \end{aligned}$$

(A.4)

where $\omega _N=\exp {\left( -{2\pi \mathrm {i}}/{N}\right) }$. The variables $\theta$ and k in formulas (A.2)-(A.4) are separated, which reduces the precomputation time. As a result, the total cost of calculating all rows of the matrix (7) will be proportional to $O(N \log N+Nq)$, where the first term is the number of operations required for a single FFT calculation for the Kaiser window, and the second term is the number of operations of N convolutions of the form (A.1) to calculate the compressed matrix elements.