Chebyshev Matrix Approximation

• Chebyshev approximation is a method to represent matrix functions using polynomial expansions with near-minimax error bounds and efficient three-term recurrences.
• It provides robust low-rank approximations via alternating minimization and generalized Remez algorithms, achieving entrywise precision even with slow singular value decay.
• The approach extends to spectral graph convolution, rational approximations for matrix exponentials, and operator learning, offering scalable and accurate computational frameworks.

Chebyshev approximation in the context of matrices refers to representing functions of matrices—especially polynomials and rational functions—using expansions in Chebyshev polynomials, and to the associated computational frameworks for low-rank, entrywise uniform, or spectral function approximation. This approach exploits the near-minimax error properties and optimal convergence rates of Chebyshev polynomials, provides scalable algorithms for large, sparse matrices, and yields theoretically sharp error bounds even in situations where singular value decay is slow. Chebyshev approximation spans polynomial expansion for computing matrix functions, rational partial fraction schemes for matrix exponentials, entrywise uniform norm (Chebyshev norm) low-rank approximation, tropical optimization for log-Chebyshev rating aggregation, and spectral operator discretizations.

1. Chebyshev Polynomial Expansion for Matrix Functions

Given a matrix $A$ (typically Hermitian or symmetric) with spectrum in $[a,b]$, the Chebyshev expansion approximates a target scalar function $f(x)$ by a truncated sum $p_K(x) = \sum_{k=0}^K c_k T_k(x)$, where $T_k$ are the Chebyshev polynomials of the first kind, defined by $$T_0(x) = 1, T_1(x) = x, T_{k+1}(x) = 2xT_k(x) - T_{k-1}(x)$$ for $k\geq1$ (Pereira et al., 2018). Mapping $A$ to $[-1,1]$ via affine scaling, $f(A)$ is approximated as $p_K(\tilde{A})$ and efficiently computed via the three-term recurrence for matrix-vector products. This method yields approximations for functions such as matrix inverse, matrix exponential, fractional powers, and more—without explicit diagonalization.

Error bounds for analytic $f$ are spectral: for $f$ analytic in a Bernstein ellipse $E_\rho$, the expansion error decays as $\mathcal{O}(\rho^{-K})$, yielding asymptotic machine precision with moderate $K$ (Pereira et al., 2018, Han et al., 2016). Applications include scalable simulation of Gaussian Markov random fields, spectral sums such as $\text{tr}\, f(A)$ (log-determinant, Estrada index), and solution of function-of-matrix equations.

2. Chebyshev Norm Low-Rank Approximation

For $A\in\mathbb{R}^{m\times n}$ and target rank $k$, Chebyshev norm approximation seeks $X$ of rank $\leq k$ with minimal entrywise error:

$$\min_{rank(X)\leq k}\|A - X\|_\infty = \min_{U\in\mathbb{R}^{m\times k}, V\in\mathbb{R}^{n\times k}} \max_{i,j}|A_{ij} - (UV^T)_{ij}|$$

A key methodology is alternating minimization: fix $V$, solve the $U$-subproblem (or vice versa), where each step involves $m$ (or $n$) independent vector Chebyshev approximation problems, efficiently solved via generalized Remez (alternance) algorithms (Morozov et al., 2022, Morozov et al., 2024). In practice, these algorithms converge monotonically, decouple in row/column blocks, and are robust to initialization for large $n$, with complexity $O((m+n)k^4 + 2mnk)$ for each iteration.

Crucially, Chebyshev-norm approximation does not require rapid singular value decay—entrywise precision $\epsilon$ is achieved at rank $k=\mathcal{O}(\log n/\epsilon^2)$ (Morozov et al., 2022). This sharply contrasts with Frobenius or spectral norm approximations, which require fast decay for small $k$, making Chebyshev approaches superior for e.g. recommender systems or kernel matrices with flat spectra.

Recent advances introduce the concept of "2-way alternance"—a necessary optimality criterion for low-rank Chebyshev solutions, and develop accelerated block solvers with $O(r^3+nrI)$ cost per update, where $I$ is the exchange count (Morozov et al., 2024). For symmetric positive definite $A$, improved entrywise error bounds are available, interpolating between $(r+1)\sigma_{r+1}(A)$ and $H(\lambda_1,\ldots,\lambda_{r+1})/\lambda_{r+1}$ depending on spectral decay rate (Hoog et al., 2021).

3. Chebyshev Interpolation, Filtering, and Spectral Graph Convolutions

In graph neural networks, spectral graph convolution can be efficiently approximated with Chebyshev polynomials of the rescaled graph Laplacian $L$ (He et al., 2022). The standard ChebNet uses polynomial filtering $$h(\lambda) \approx \sum_{k=0}^K w_k T_k(\hat{\lambda})$$ and applies convolutions via repeated multiplications of $T_k(\tilde{L})x$. While theoretically optimal for analytic filters, unconstrained learning of $w_k$ may violate coefficient decay and induce over-fitting (Runge phenomenon).

ChebNetII generalizes this by using Chebyshev interpolation: filter values are learned at $K+1$ Chebyshev nodes $x_j = \cos\left(\frac{2j+1}\pi{2K+2}\right)$, and weights reconstructed via

$w_k = \frac{2}{K+1}\sum_{j=0}^K \gamma_j T_k(x_j)$

This enforces analytic decay, reduces polynomial oscillation, and yields near-minimax uniform-norm approximation with complexity $O(K|E|)$ per pass. Empirical results show ChebNetII matches or exceeds performance of GCN, GPR-GNN, and BernNet in node classification, scaling up to billion-node graphs (He et al., 2022).

4. Chebyshev Rational and Fractional Approximation for Matrix Exponentials and Fractional Operators

Chebyshev rational approximation provides near-minimax uniform accuracy for $e^x$ on $(−\infty,0]$; the partial fraction decomposition enables efficient computation of $e^{A t}$ for matrices with spectrum near the negative real axis (Pusa, 2012). With corrected coefficients, the uniform error satisfies $E_n = \mathcal{O}(H^{-n})$ (Halphen constant $H \approx 9.289$). Evaluation for $A$ involves $n$ shifted solves $(A t - \theta_j I)^{-1}x_0$ for the precomputed poles $\theta_j$ and residues $\alpha_j$:

$$e^{A t}x_0 \approx C_n x_0 + 2\,\mathrm{Re} \sum_{j=1}^{n/2} \alpha_j (A t - \theta_j I)^{-1}x_0$$

This algorithm is widely used in nuclear burnup and parabolic PDE solvers.

For fractional integral operators,

$$\mathcal{I}^{\alpha} u(x) = \frac{1}{\Gamma(\alpha)} \int_{-1}^x \frac{u(t)}{(x-t)^{1-\alpha}} dt$$

mapped Chebyshev (fractional Jacobi) polynomials are used to build well-conditioned matrices approximating $\mathcal{I}^\alpha$ (Liu et al., 24 Jun 2025). A three-term recurrence for fractional integrals combined with a banded-plus-low-rank structure ensures spectral convergence and $O(N^2)$ complexity, applicable to boundary-value, fractional differential, and eigenvalue problems.

5. Log-Chebyshev Matrix Approximation and Tropical Optimization

In pairwise-comparison/rating aggregation, log-Chebyshev matrix approximation seeks a consistent rank-1 matrix $X_{ij}=x_i/x_j$ close in the worst-case log sense to given reciprocal matrices $A^{(k)}$ (Krivulin, 25 Jan 2026, Krivulin, 2022). The objective

$$\min_{x>0} \max_{i,j,k} \left| \log A^{(k)}_{ij} - \log \frac{x_i}{x_j} \right|$$

is reformulated using max-plus tropical algebra, where optimization reduces to minimizing $x^{-}A x$ over $x>0$. Analytical solutions are constructed using the tropical spectral radius and Kleene-star $\left(\lambda^{-1}A\right)^*$, yielding full solution cones and best/worst differentiating vectors. This approach directly controls worst-case deviation and supports multicriteria analysis with box constraints, contrasting with AHP and weighted geometric mean which optimize sum-of-squares criteria.

6. Spectrum-Adapted Chebyshev Approximation and Discrete Operator Realization

Classical Chebyshev truncation is spectrum-blind: accuracy is uniform across the interval, possibly wasting polynomial order away from densely clustered eigenvalues. Spectrum-adapted Chebyshev schemes first estimate the spectral density of $A$ (stochastic trace/Heaviside approximation), then warp nodes or fit weights where eigenvalues concentrate (Fan et al., 2018). Two constructions are effective:

• Inverse-CDF warping of Chebyshev nodes to regions of dense eigenvalues, followed by barycentric interpolation.
• Weighted least squares using orthogonal polynomials formed with weights set to the local spectral density.

This yields substantially better spectral error, especially where $A$ has interior clusters of eigenvalues. In graph signal processing and large-scale matrix computation, these methods improve over Lanczos and classic Chebyshev approaches at moderate $K$.

7. Koopman Matrix Construction by Chebyshev Expansion

Chebyshev polynomial projection is utilized for operator learning—e.g., discovery of governing PDEs via Koopman operator approximation (Sisaykeo et al., 27 Oct 2025). Sampled data or differential operators are projected onto a Chebyshev basis to build finite-dimensional Koopman matrices:

1. Sample at $N+1$ Gauss–Lobatto nodes: $x_j = \cos(\pi j/N)$.
2. Compute basis transforms (DCT-II), differentiation matrices, and operator representations by Galerkin or collocation.
3. Estimate the Koopman matrix from data: $\hat{K} = A_1 A_0^\dagger$.
4. Spectral convergence is achieved exponentially in $N$ for analytic problems, allowing direct spectrum linking between discrete matrices and continuous operators.

This approach links finite matrix approximation to PDE identification, operator learning, and modal analysis.

---

Summary Table: Chebyshev Matrix Approximation Modes

Mode	Key Algorithmic Elements	Main Paper
Matrix Polynomial Function Approximation	Three-term recurrence, scaling	(Pereira et al., 2018)
Low-Rank $\|\cdot\|_\infty$ Approximation	Alternating minimization, Remez	(Morozov et al., 2022, Morozov et al., 2024)
Spectral Graph Filter Approximation	Chebyshev interpolation, node-wise learning	(He et al., 2022)
Rational Matrix Exponential	PFD, shifted linear solves	(Pusa, 2012)
Fractional Integral Operator	Fractional Chebyshev, band+sparse	(Liu et al., 24 Jun 2025)
Log-Chebyshev/Tropical Ranking	Spectral radius, Kleene-star	(Krivulin, 25 Jan 2026, Krivulin, 2022)
Spectrum-Adapted Polynomial	Spectral density fitting, warping	(Fan et al., 2018)
Koopman Operator Construction	DCT, collocation, Galerkin	(Sisaykeo et al., 27 Oct 2025)

Impact and Current Directions

Chebyshev matrix approximation frameworks underpin major algorithmic advances in matrix computations: enabling efficient, high-accuracy simulation, scalable function evaluation in high dimensions, robust entrywise error control unconstrained by spectral decay, and interpretable operator learning. Current research directions include enhanced tropical optimization for multicriteria aggregation, spectral adaptation strategies, fast banded operators for fractional calculus, and deeper integration with data-driven modal analysis in dynamical systems. The area continues to evolve with new fast solvers, optimality characterizations (alternance), and finer error analysis for both analytic and structured matrices.