Fractional Differentiation Matrices

• Fractional differentiation matrices are finite-dimensional operators that discretize non-integer order derivatives using binomial coefficients and spectral collocation methods.
• They employ various matrix constructions, including Toeplitz, banded, and Kronecker products, to accurately approximate both fractional ODEs and PDEs.
• Advanced applications span neural network optimization and matrix-order differentiation, unifying continuous fractional operators within a robust computational framework.

Fractional differentiation matrices are finite-dimensional operators designed to discretize and approximate fractional derivatives—operators of non-integer order—on grids, polynomial bases, or function spaces. These matrices underpin a wide range of numerical methods for fractional ordinary and partial differential equations, including approaches based on finite differences, Toeplitz structure, spectral collocation, and operational matrices. In recent developments, they also serve as key algorithmic tools in neural network optimization with fractional-order derivatives and as discrete analogues for more general operator powers, including matrix-order differentiation.

1. Classical Discretizations: Binomial, Grünwald-Letnikov, and Toeplitz Structures

Foundational approaches to fractional differentiation matrices originate from discrete analogues of the Grünwald-Letnikov formula and the properties of Toeplitz matrices.

Discrete Binomial Approximations:

Let α > 0 be the fractional order and t ∈ {0,1,…,N}. Discrete Grünwald-type definitions result in lower- and upper-triangular (N+1) × (N+1) matrices with binomial weights:

• Left fractional difference (delta or nabla):

$$D_{i,j}^{(\alpha)} = (-1)^{i-j}\binom{\alpha}{i-j}\mathbf{1}_{j\le i}$$

• Right fractional difference:

$$D_{i,j}^{(\alpha)} = (-1)^{j-i}\binom{\alpha}{j-i}\mathbf{1}_{j\ge i}$$

where $$\binom{\alpha}{k} = \alpha(\alpha-1)\cdots(\alpha-k+1)/k!$$ for $k \ge 0$.

For a grid function $\vec{f} = (f_0, ..., f_N)^T$, the product $D^{(\alpha)}\vec{f}$ computes the fractional difference at each node. These matrices, and their definitions, extend to both delta and nabla operators, with practical boundary conventions (e.g., $f_j=0$ outside the domain) and direct implementation pathways (Abdeljawad et al., 2012).

Toeplitz Matrix Discretizations:

Fractional differentiation and integration operators are also realized through Toeplitz matrix constructions. Given a symbol $\varphi_{α,R}(θ) = (1-Re^{iθ})^α(1+Re^{-iθ})^α$, the Toeplitz matrix $T_N(\varphi_{α,R})$ with Fourier-coefficient entries approximates (on the uniform grid) the action of the Marchaud/Riemann–Liouville derivative. The explicit asymptotic decay of the symbol’s coefficient $$\varphi_{α}(n) \sim 2^α/\Gamma(-\alpha)(-1)^n|n|^{-\alpha-1}$$ ensures the convergence of $N^\alpha T_N(\varphi_\alpha)\vec{f}$ to the continuous fractional derivative. Inversion of this Toeplitz matrix produces a discrete fractional integral operator corresponding to the Riemann–Liouville kernel, and all classical analytic relations (such as those for Green’s kernels and integer-order reduction) are recovered in the matrix framework (Rambour et al., 2018).

2. Fractional Differentiation Matrices via Orthogonal Polynomials and Spectral Collocation

High-order discretizations leverage the spectral properties of orthogonal polynomials, especially Jacobi-type polynomials on Jacobi–Gauss–Lobatto (JGL) grids.

Given a polynomial interpolation basis $F_j(x)$ (Lagrange at JGL points $x_k$), any function $u(x)$ is expanded as $$I_N u(x) = \sum_{j=0}^N u(x_j) F_j(x) = \sum_{k=0}^N p_k P_k^{a,b}(x)$$. The construction of the fractional Caputo or Riemann-Liouville differentiation matrix involves several recurrence and expansion steps:

• Compute recurrence relations for fractional integrals of Jacobi polynomials,
• Form integral/differentiation matrices (e.g., $_C\widehat D_L^{(a,b,\alpha)}$ for left Caputo),
• Apply conversion matrices from nodal values to the polynomial basis (matrix $M$),
• Combine to yield $$_C D_L^{(a,b,\alpha)} = _C\widehat D_L^{(a,b,\alpha)} M$$.

These matrices, when applied, realize high-accuracy spectral collocation methods. Spectral-radius analysis shows that the interior differentiation matrices $M_i^{(\alpha,a,b)}$ scale as $N^{2\alpha}$ in spectral radius with uniform $O(1)$ proportionality constants, indicating stable and spectrally accurate discretizations. Empirical tests on Baglay–Torvic and fractional diffusion equations demonstrate spectral (exponential) convergence and surpass the accuracy of alternative operational-matrix and finite-difference schemes for moderate polynomial degree $N$ (Zeng et al., 2014).

3. Matrix Formulations for PDEs: Banded-Strip and Operational Matrices

For time- and space-fractional PDEs, matrix discretizations generalize to Kronecker-product structures, banded matrices, and precomputable operational matrices.

Grünwald-Type Banded Strip Matrices:

On a uniform time grid $t_j = j\tau$, the backward-difference approximation for a Caputo or Riemann–Liouville fractional derivative yields a banded strip matrix:

$$B_n^{(\alpha)} = \tau^{-\alpha} \begin{pmatrix} \omega_0 & \omega_1 & \cdots & \omega_n\ 0 & \omega_0 & \cdots & \omega_{n-1}\ \vdots & \ddots & \ddots & \vdots\ 0 & \cdots & 0 & \omega_0\ \end{pmatrix}$$

with $\omega_k^{(\alpha)} = (-1)^k \binom{\alpha}{k}$. The discrete operator for a space-fractional Riesz derivative is built as a symmetric banded matrix by taking specific linear combinations or direct weights (e.g., Ortigueira’s formula) (0811.1355).

Kronecker Product Construction:

To handle grids in multi-dimensions (e.g., for $u(x_i,t_j)$), the full operator is assembled as $$A = B_n^{(\alpha)} \otimes I_{m+1} - \chi I_{n+1} \otimes R_m^{(\beta)}$$, enabling direct encoding of time and space derivatives, with boundary and delay conditions accommodated via eliminators and shifters.

Caputo Operational Matrices for PINNs:

The Caputo derivative can be discretized on arbitrary grids using telescoping auxiliary weights, producing a lower-triangular operational matrix $A$ such that for vectorized function values $\vec{f}$, the approximation $\vec{f}^{(\alpha)} \approx A\vec{f}$ holds. This supports non-uniform domains and is particularly suitable for embedding in physics-informed neural networks (PINNs), enabling matrix-vector multiplications in place of convolution quadrature steps, and admits integration with polynomial neural blocks (e.g., Legendre Block) (Taheri et al., 2024).

4. Extensions: Matrix-Order Fractional Differentiation and Jordan Functional Calculus

Fractional calculus admits further extension where the order of differentiation itself is an $n \times n$ complex matrix $A$, with differentiations defined via matrix-valued gamma and power functions.

Consider the Riemann–Liouville matrix-order derivative for a vector function $F(t)$:

$$D_{0+}^A F(t) = \frac{1}{\Gamma(-A)} \int_0^t (t-\tau)^{-A-I} F(\tau) d\tau,$$

with Gamma and power functions applied via the Jordan decomposition of $A$. Discrete analogues leverage block Toeplitz matrices assembled from the matrix-valued binomial coefficients:

$$\binom{A}{k} = \frac{A(A-I)\cdots (A-(k-1)I)}{k!}\,,$$

enabling the construction of matrix-difference operators $D_h^A$ approximating $D_{0+}^A$. Applications involve systems of viscoelasticity, coupled anomalous diffusions, and vector-valued control and econometric models, with all analytical expressions (integral representations, inverse relations, semigroup laws) extending naturally to the matrix-order case (Porciuncula, 2020).

5. Algebraic and Implementation Properties

Triangular Structures and Spectral Stability:

• Left-fractional matrices (lower-triangular) and right-fractional matrices (upper-triangular) are well-posed for forward and backward time-stepping.
• For Toeplitz and banded matrices, spectral analysis and factorization via Hankel–Toeplitz techniques guarantee invertibility and correct analytic behavior as grid size increases (Rambour et al., 2018).

Convergence and Accuracy:

• Grünwald-type time discretizations yield $O(\tau)$ accuracy for $0<\alpha<1$; Riesz spatial approximations provide $O(h^2)$ accuracy.
• Spectral collocation with fractional differentiation matrices achieves exponential convergence in $N$ for smooth problems and outperforms low-order schemes (Zeng et al., 2014).

Numerical Recipes:

• For uniform or non-uniform grids, weights and matrix entries are constructed via binomial coefficients, spectral polynomial recurrences, or integral formulas.
• Implementation details (pseudocode) are available for a variety of scenarios—see (0811.1355, Kolokoltsov et al., 11 Dec 2025, Taheri et al., 2024).

6. Applications and Computational Impact

Fractional differentiation matrices are employed in:

Domain	Matrix Approach	Reference
Fractional ODE/PDE Solvers	Banded, Toeplitz, Spectral Collocation, Kronecker	(0811.1355, Zeng et al., 2014, Rambour et al., 2018)
PINNs for FDEs	Caputo operational matrices in neural blocks	(Taheri et al., 2024)
Matrix-order FDEs	Block Toeplitz and Jordan-calculus based matrices	(Porciuncula, 2020)
Operator Powers	Diagonalization, Balakrishnan’s formulae	(Kolokoltsov et al., 11 Dec 2025)
Autograd in DNN training	Fractional-order Jacobians (block matrices)	(zhou et al., 9 Jun 2025)

In PINNs and deep-learning contexts, fractional-order differentiation matrices provide efficient means for integrating non-integer order constraints directly into the training loop, lead to intrinsic regularization and memory properties, and are shown empirically to accelerate or improve test set metrics for certain time-series tasks (zhou et al., 9 Jun 2025, Taheri et al., 2024).

7. Synthesis and Theoretical Unification

Fractional differentiation matrices unify and discretize all classical continuous fractional operators (Riemann–Liouville, Caputo, Weyl, Grünwald–Letnikov, Marchaud), as well as their integer-order limits, within a single matrix-based numerical and analytic framework. This extends to matrix-order and operator-function settings, spectral methods, neural operator learning architectures, and iterative solvers leveraging Kronecker and Toeplitz structures. Their properties, stability, and convergence have been analyzed across finite-difference, spectral, and functional-analytic perspectives, demonstrating the adaptability and efficacy of the matrix approach for fractional calculus in numerical and computational analysis (Rambour et al., 2018, 0811.1355, Zeng et al., 2014, Porciuncula, 2020, Taheri et al., 2024, zhou et al., 9 Jun 2025, Kolokoltsov et al., 11 Dec 2025, Abdeljawad et al., 2012).