A Split Fast Fourier Transform Algorithm for Block Toeplitz Matrix-Vector Multiplication

La modélisation des phénomènes physiques invariants par translation, dans une base vectorielle spatiale, entraine des systèmes d’équations qui peuvent être représentés par des matrices de Toeplitz par bloc multi-niveaux. C’est notamment le cas des phénomènes liés à l’électromagnétisme ou l’acoustique. La résolution d’équations impliquant ces matrices requiert généralement une méthode de résolution itérative, nécessitant un grand nombre de multiplications de matrice de Toeplitz par bloc avec des vecteurs. Ainsi, pour simuler des problèmes physiques de grande taille, il est indispensable d’améliorer l’efficacité de ces opérations. Il existe un grand nombre de méthodes réalisant ces multiplications, mais chacune présente des limites majeures. Les plus prometteuses de ces méthodes sont celle par insertion de matrices circulantes et celle par décomposition tensorielle en paquet. Malheureusement, ces méthodes présentent des temps de calculs et des besoins en mémoire qui augmentent rapidement avec la taille du système ou le rang du tenseur. Dans le cas de vecteurs sources quelconques, le calcul de la décomposition tensorielle est plus coûteux que la multiplication directe. Aussi, lorsque l’on étend l’insertion de matrices circulantes aux matrices de Toeplitz par bloc en d dimensions, l’insertion entraine l’introduction d’un grand nombre de coefficients redondants. Dans ce cas, seul une proportion de 1/2d du vecteur traité dans le calcul contient des informations utiles. Cette redondance conduit rapidement à une surcharge de mémoire et à une complexité opérationnelle élevée. L’algorithme présenté dans ce mémoire a pour objectif de répondre à ce problème d’inefficacité en reportant au plus tard les insertions et en avançant au plus tôt les projections. Cette stratégie permet de calculer les transformées de Fourier sur un nombre réduit de coefficients, diminuant ainsi la complexité par un facteur de d/ 2 − 2−d+1 et la demande en mémoire maximale de 2/ (d + 1)2−d+1. Pour une multiplication matrice-vecteur en trois dimensions, le ratio entre la complexité théorique de la méthode standard par insertion de matrices circulantes et celle de la version améliorée de cette méthode tend vers 12/7 lorsque la taille de la matrice de Toeplitz par bloc tend vers l’infini pour chaque niveau. De plus, l’algorithme utilise une propriété de la transformée de Fourier rapide afin de séparer le vecteur transformé en deux vecteurs de coefficients d’indices paires et impaires respectivement, et ce pour chaque transformée de Fourier. Cela résulte en une structure arborescente des tâches de l’algorithme.

Abstract

Modeling translationally invariant physics, such as electromagnetism, in a spatial basis typ-ically results in systems of equations involving multi-level block-Toeplitz matrices. These matrices often necessitate iterative solution methods, which in turn require numerous block-Toeplitz matrix-vector multiplications. Therefore, enhancing the efficiency of these multi-plications is essential for simulating a wide range of large-scale physical systems. While various methods for performing these multiplications currently exist, each comes with major drawbacks. Specifically, the most promising methods—circulant embedding and tensor train decomposition—suffer from poor scalability with respect to system size or tensor rank. For general source vectors, computing the tensor decomposition is more computationally inten-sive than direct multiplication. When extended to d-dimensional block-Toeplitz matrices, embedding leads to the introduction of a large number of redundant coefficients, so that only 1/2d of the vector treated in computation contains useful information. This scaling quickly leads to memory overload and high operational complexity. The algorithm introduced in this thesis aims to address this latter inefficiency through lazy embeddings and eager projections, deferring embeddings as late as possible and performing projections as early as possible. This strategy allows for the computation of Fourier trans-forms on a reduced number of coefficients, lowering complexity by a factor of d/ 2 − 2−d+1and peak memory usage by a factor of 2/ (d + 1)2−d+1. For three-dimensional matrix-vector multiplications, the theoretical complexity ratio of the standard circulant embedding method over the enhanced version approaches 12/7 as the size of the block-Toeplitz matrix tends to infinity for each level. Additionally, the algorithm exploits a property of the fast Fourier transform (FFT) to divide the transformed vector into even and odd coefficients, fol-lowing a tree branch structure. This structure offers greater flexibility in managing memory and computational speed through various parallelization strategies along each branch using several GPUs. Running these branches independently on separate GPUs can boost compu-tational speed at the expense of increased memory usage. On a single GPU, the memory usage ratio between the standard and enhanced methods tends to 4/3 for three-dimensional matrices and vectors as the sizes of each level increase. In cases where the Toeplitz data is entirely symmetric or anti-symmetric—meaning each block at each level is a symmetric or anti-symmetric Toeplitz matrix–—a clever embedding can reduce memory consumption. The theoretical memory ratio between the two methods then becomes (2d +1)/(d+2).