Univariate Taylor Propagation of Matrices (UTPM)

• UTPM is a framework that lifts matrix variables into a polynomial ring to enable efficient higher-order derivative computation.
• It employs unified forward and reverse automatic differentiation using explicit recurrences for operations like inversion, QR, and eigen-decomposition.
• UTPM offers significant speedup and reduced memory usage compared to elementwise methods, achieving machine precision even in complex systems.

Univariate Taylor Propagation of Matrices (UTPM) is a mathematical and algorithmic framework for the efficient computation of higher-order derivatives of matrix-valued functions. It generalizes univariate Taylor propagation on scalars (UTPS) to the matrix case by lifting variables into a polynomial ring whose coefficients are matrices, enabling forward-mode automatic differentiation (AD) and reverse-mode AD for programs composed of matrix operations. UTPM is particularly effective for linear algebraic routines such as inversion, multiplication, QR decomposition, and eigenvalue decomposition, enabling high-order derivative propagation through composite systems and supporting both analytic matrix functions and complex workflows in scientific computation (0911.4940, Walter et al., 2010, Ramirez et al., 2020).

1. Mathematical Framework

Let $\mathbb{F}$ be a ring (typically $\mathbb{R}$), and let $\mathbb{F}^{M\times N}$ denote the set of $M\times N$ matrices over $\mathbb{F}$. In UTPM, each matrix $X$ is replaced by a truncated Taylor expansion: $$[X]_D = X^{[0]} + X^{[1]}t + \cdots + X^{[D-1]}t^{D-1}$$ where each coefficient $X^{[k]} \in \mathbb{F}^{M\times N}$. Given a matrix-valued function $$F: \mathbb{F}^{M\times N} \to \mathbb{F}^{P\times Q}$$, the Taylor expansion becomes: $$F\left(X(t)\right) = \sum_{d=0}^{D-1} F^{[d]} t^d + O(t^D)$$ where

$\mathbb{R}$0

This lifting enables the algebraic manipulation of Taylor series, supporting higher-order directional derivatives in a single forward sweep (0911.4940, Walter et al., 2010).

The UTPM algebra is closed under the “elementary” matrix functions, with explicit recurrences for each basic operation:

• Matrix Sum: $\mathbb{R}$1
• Matrix Product: $\mathbb{R}$2
• Matrix Inverse: For $\mathbb{R}$3, $\mathbb{R}$4 and for $\mathbb{R}$5,

$\mathbb{R}$6

• Matrix Diagonalization (non-degenerate case): If $\mathbb{R}$7 has simple eigenpairs $\mathbb{R}$8, then the Taylor coefficients $\mathbb{R}$9 admit explicit recursions based on perturbation theory (Ramirez et al., 2020).

2. Algorithmic Realization and Modes

The UTPM workflow consists of a unified forward and reverse pass:

• Forward sweep: The truncated Taylor series for each variable is propagated through the computational graph using the recurrence rules for each supported matrix operation.
• Reverse sweep: Adjoint series are propagated backward using the appropriate reverse-mode rules, which are also defined in coefficient form (e.g., for matrix multiplication or inversion).

This dual formulation allows one sweep to capture all derivatives up to order $\mathbb{F}^{M\times N}$0 for all inputs in the program and immediately obtain gradients, Hessians, and higher-order analogs (0911.4940, Walter et al., 2010).

Specializations include:

• QR decomposition (rectangular): A closed-form system is solved for the Taylor coefficients of $\mathbb{F}^{M\times N}$1 and $\mathbb{F}^{M\times N}$2, propagating the series subject to the defining relations of the decomposition (Walter et al., 2010).
• Eigenvalue decomposition (symmetric, non-degenerate): Taylor series for eigenvalues and eigenvectors are computed recursively, and their derivatives are used in composite objectives or loss functions.

3. Integration with Composite Systems

UTPM is used to compute high-order Taylor expansions in complex algorithms such as the Recursive Scattering Matrix Method (RSMM) for quantum transport (Ramirez et al., 2020). In RSMM, the scattering matrix of a composite system is assembled from the Taylor expansions of subsystem matrices using blockwise sum, product, and inversion: $\mathbb{F}^{M\times N}$3 with $\mathbb{F}^{M\times N}$4, where all series manipulations use UTPM recurrences. For periodic systems, diagonalization and functions of matrices (e.g., the matrix exponential, logarithm, resolvent) admit UTPM recursions, and the same algebraic machinery applies (Ramirez et al., 2020).

4. Complexity, Efficiency, and Numerical Stability

The computational and storage costs of UTPM are substantially reduced compared to elementwise UTPS:

• Forward mode: $\mathbb{F}^{M\times N}$5 for a matrix of size $\mathbb{F}^{M\times N}$6 and polynomial order $\mathbb{F}^{M\times N}$7.
• Memory: $\mathbb{F}^{M\times N}$8 for storing the Taylor coefficients for each matrix in the computation (0911.4940, Walter et al., 2010).

Given a composite algorithm with $\mathbb{F}^{M\times N}$9 matrix-merges (e.g., RSMM), total cost is $M\times N$0 for polynomial order $M\times N$1 (Ramirez et al., 2020). In contrast, UTPS has $M\times N$2 computational overhead due to inefficient scalar-level bookkeeping.

Empirical performance for sample problems (e.g., $M\times N$3, $M\times N$4) shows roughly $M\times N$5 speedup using UTPM over UTPS (e.g., 0.45 s vs. 45 s) and a reduction in memory usage by more than an order of magnitude (0911.4940). UTPM yields derivative coefficients to within $M\times N$6 of analytic or finite-difference values, i.e., machine precision even in ill-conditioned cases (Walter et al., 2010).

5. Analytical Error, Convergence, and Domain of Validity

UTPM Taylor expansions inherit the local error and convergence properties of Taylor series. The Lagrange remainder for the series,

$M\times N$7

is bounded by $M\times N$8. The convergence radius $M\times N$9 is determined by the distance to the nearest complex pole of the analytic function (e.g., a singularity of the matrix being expanded) (Ramirez et al., 2020). In practical algorithms, expansion order $\mathbb{F}$0 is chosen to ensure $\mathbb{F}$1 machine $\mathbb{F}$2 within the domain $\mathbb{F}$3. If $\mathbb{F}$4 exceeds this interval, the expansion center must be shifted and the propagation restarted.

6. Implementation and Applications

UTPM has been implemented in the BSD-licensed Python library ALGOPY (Walter et al., 2010). Key features:

• Matrix-polynomial (UTPM) and scalar-polynomial (UTPS) factor rings for overloading arithmetic
• Explicit push-forward and pullback routines for QR, eigen, inversion, and related matrix functions
• Graph-mode reverse accumulation requiring only storage of function nodes and Jacobians in polynomial form, not all scalar operations

Numerical applications include efficient high-order derivatives for design optimization, mesoscopic transport computations (Landauer conductance), and quantum device simulation. In optimum experimental design, UTPM enables stable machine-precision gradients for complex matrix objectives such as maximization of the largest eigenvalue of composite matrix products (Walter et al., 2010).

7. Generality, Limitations, and Extensions

The UTPM recurrences—sum, product, inverse, diagonalization—apply to any analytic matrix-valued function, including $\mathbb{F}$5, $\mathbb{F}$6, and the matrix resolvent. The theory extends to non-polynomial analytic functions using similar recursions (Ramirez et al., 2020).

UTPM is limited by the need for explicit recurrences for each new matrix operation; for very high polynomial orders ($\mathbb{F}$7), computational cost may become prohibitive. The method assumes fixed dataflow and analytic dependence; control-flow dependent on data and generalized differentiation for non-analytic routines require augmentation by standard AD strategies (0911.4940).

---

Selected References:

• Efficient Higher Order Derivatives of Objective Functions Composed of Matrix Operations (0911.4940)
• Algorithmic Differentiation of Linear Algebra Functions with Application in Optimum Experimental Design (Walter et al., 2010)
• Taylor series of Landauer conductance (Ramirez et al., 2020)