Computing the action of a matrix exponential on an interval via the $\star$-product approach

Probleemstelling

Het berekenen van de actie van de matrixexponentiële op een vector, $e^{At}v$, is een fundamentele taak in vele wetenschappelijke domeinen, waaronder kwantumdynamica, besturingstheorie en de tijdsintegratie van grote lineaire differentiaalvergelijkingen (ODE's) die voortvloeien uit ruimtelijke discretisaties van partiële differentiaalvergelijkingen (PDE's).

De bestaande methoden (zoals expv of expmv_tspan) zijn vaak geoptimaliseerd voor het berekenen van de oplossing op specifieke tijdstippen of op een gelijkmatig verdeeld rooster. Een uitdaging is echter het efficiënt en nauwkeurig evalueren van de oplossing voor elk tijdstip $t$ binnen een vooraf bepaald, begrensd interval. Traditionele methoden vereisen vaak herhaalde berekeningen of interpolatie als men de oplossing op willekeurige punten binnen het interval nodig heeft.

Methodologie

De auteurs introduceren een nieuwe numerieke methode, de $\star$-methode, die de oplossing uitbreidt in een reeks van orthogonale polynomen (specifiek Legendre-polynomen) over het interval. De kern van de aanpak is als volgt:

1. $\star$-Algebra en Distributies:
   De methode maakt gebruik van de recente theorie van de $\star$-algebra, een algebra van bivariate distributies. De matrixexponentiële wordt hierin voorgesteld via een veralgemening van de Volterra-convolutie (de $\star$-product). De differentiaalvergelijking $\dot{u}(t) = Au(t)$ wordt herschreven in termen van een bivariate operator $A(t,s)$ die de Heaviside-stapfunctie $\Theta(t-s)$ bevat.

2. Spectrale Benadering:
   De oplossing wordt benaderd door een afgekapt dubbel Legendre-polynoom expansie. De coëfficiënten van deze expansie worden verzameld in een matrix $U_M$. Door het $\star$-product te discretiseren met behulp van Legendre-polynomen, wordt het probleem van het benaderen van de matrixexponentiële omgezet in een lineair algebraïsch probleem.

3. Stein-vergelijking:
   Het resulterende systeem leidt tot een lineaire vergelijking die equivalent is aan een Stein-type matrixvergelijking (ook wel een discrete Lyapunov-vergelijking genoemd):
   $$X - \frac{1}{2}T_M X \tilde{A}^T = \phi_M(s)v^T$$
   Hierbij is $T_M$ de coëfficiëntenmatrix van de Heaviside-functie in de Legendre-basis, en $\phi_M(t)$ de vector van de Legendre-polynomen.

   • Voor kleine tot middelgrote matrices wordt deze vergelijking direct opgelost (bijv. via Schur-decompositie en de dlyap-routine in MATLAB).
   • Voor grote, schaarse matrices wordt de methode gecombineerd met de Arnoldi-algoritme (Krylov-ruimteprojectie). Hierbij wordt de vergelijking geprojecteerd op een lagere dimensie Krylov-ruimte, wat de rekentijd drastisch verlaagt zonder in te leveren op de nauwkeurigheid.

4. Foutanalyse:
   De auteurs leiden een bovengrens af voor de afkappelfout (truncation error). Ze tonen aan dat de fout super-exponentieel afneemt met het aantal termen $M$ in de expansie, mits bepaalde voorwaarden voor de spectrale radius van de matrix $A$ en de parameter $c$ in de Taylor-reeks worden voldaan.

Belangrijkste Bijdragen

• Nieuwe Representatie: Een nieuwe formulering van de matrixexponentiële in de $\star$-algebra die leidt tot een Stein-vergelijking.
• Intervaloplossing: De methode levert een continue representatie van de oplossing op via Legendre-coëfficiënten. Dit betekent dat de oplossing voor elk $t$ in het interval kan worden berekend zonder de onderliggende lineaire vergelijking opnieuw op te hoeven lossen.
• Schalbaarheid: De integratie met Krylov-ruimte methoden (Arnoldi) maakt de methode toepasbaar op grote, schaarse matrices, wat essentieel is voor praktische toepassingen zoals FEM-discretisaties.
• Foutbound: Een theoretische bovengrens voor de afkappelfout die de convergentie van de methode garandeert.

Resultaten

De auteurs hebben de methode getest op een breed scala aan matrices, waaronder:

• 2D Poisson-matrices (FEM-discretisaties).
• Complexe tridiagonale matrices.
• Dichte matrices met eigenwaarden op Chebyshev-nodes.
• Schrödinger-vergelijkingen (FEM semidiscretisatie).

Vergelijking met state-of-the-art:
De resultaten zijn vergeleken met expv (Krylov-projectie voor één tijdstip) en expmv_tspan (Schaling-en-kwadrering met Taylor-reeks voor een rooster).

• Nauwkeurigheid: De $\star$-methode (zowel met als zonder Arnoldi) bereikt een nauwkeurigheid die vergelijkbaar is met de beste bestaande methoden (relatieve fouten in de orde van $10^{-10}$ tot $10^{-15}$).
• Efficiëntie:
  • De pure $\star$-methode (zonder Krylov) is traag voor zeer grote matrices.
  • De $\star$-methode + Arnoldi is zeer competitief. In veel gevallen is het sneller dan of even snel als expmv_tspan, terwijl het het voordeel biedt dat de oplossing op willekeurige tijdstippen beschikbaar is zonder extra kosten.
  • De methode is bijzonder efficiënt wanneer de oplossing op veel of onbekende tijdstippen nodig is, omdat de "rekenwerk" (het oplossen van de Stein-vergelijking) slechts één keer nodig is.

Betekenis en Toekomstperspectief

Deze paper biedt een krachtig alternatief voor traditionele methoden voor het berekenen van matrixexponentiële acties over een interval. De belangrijkste meerwaarde is de flexibiliteit: de oplossing wordt niet gebonden aan een vooraf gedefinieerd rooster, maar is een continue functie (benaderd door polynomen) die direct uit de berekende coëfficiënten kan worden afgeleid.

Dit maakt de methode ideaal voor toepassingen waar:

1. De evaluatietijdstippen niet a priori bekend zijn.
2. De oplossing op een zeer fijn tijdsrooster nodig is (bijv. voor visualisatie of controle).
3. Men werkt met grote, schaarse systemen uit PDE-discretisaties.

De auteurs wijzen erop dat de huidige foutbound theoretisch correct is maar in de praktijk conservatief kan zijn, en dat toekomstig onderzoek gericht zal zijn op het ontwikkelen van scherpere bounds en het ontwikkelen van nieuwe preconditioneringstechnieken voor de onderliggende lineaire systemen.