Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

Mixed-Precision Paterson--Stockmeyer Method for Evaluating Polynomials of Matrices

Xiaobo Liu

TL;DR

This work introduces a mixed-precision Paterson–Stockmeyer framework for evaluating matrix polynomials with decaying scalar coefficients, enabling lower-precision computations for small-magnitude terms while maintaining high-precision accuracy. It provides a formal rounding-error analysis and practical algorithms (for Taylor, Padé, and cosine-type polynomials) that integrate with scaling-and-squaring strategies for the matrix exponential. The authors demonstrate that, under reasonable decay and norm conditions, the mixed-precision PS method achieves comparable accuracy to fixed-precision PS but with substantial speed-ups, especially for large, high-precision problems. They validate the approach through extensive numerical experiments and profiling, highlighting how precision scheduling and reduced inner multiplications yield notable performance gains in arbitrary-precision contexts. The framework is broadly applicable to other matrix polynomials and paves the way for further refinements, such as extending mixed-precision techniques to the computation of the first $s$ powers of $X$ and exploring more aggressive precision hierarchies.

Abstract

The Paterson--Stockmeyer method is an evaluation scheme for matrix polynomials with scalar coefficients that arise in many state-of-the-art algorithms based on polynomial or rational approximation, for example, those for computing transcendental matrix functions. We derive a mixed-precision version of the Paterson--Stockmeyer method that is particularly useful for evaluating matrix polynomials with scalar coefficients of decaying magnitude. The new method is mainly of interest in the arbitrary precision arithmetic, and it is attractive for high-precision computations. The key idea is to perform computations on data of small magnitude in low precision, and rounding error analysis is provided for the use of lower-than-the-working precisions. We focus on the evaluation of the Taylor approximants of the matrix exponential and show the applicability of our method to the existing scaling and squaring algorithms. We also demonstrate through experiments the general applicability of our method to other problems, such as computing the polynomials from the Padé approximant of the matrix exponential and the Taylor approximant of the matrix cosine. Numerical experiments show our mixed-precision Paterson--Stockmeyer algorithms can be more efficient than its fixed-precision counterpart while delivering the same level of accuracy.

Mixed-Precision Paterson--Stockmeyer Method for Evaluating Polynomials of Matrices

TL;DR

This work introduces a mixed-precision Paterson–Stockmeyer framework for evaluating matrix polynomials with decaying scalar coefficients, enabling lower-precision computations for small-magnitude terms while maintaining high-precision accuracy. It provides a formal rounding-error analysis and practical algorithms (for Taylor, Padé, and cosine-type polynomials) that integrate with scaling-and-squaring strategies for the matrix exponential. The authors demonstrate that, under reasonable decay and norm conditions, the mixed-precision PS method achieves comparable accuracy to fixed-precision PS but with substantial speed-ups, especially for large, high-precision problems. They validate the approach through extensive numerical experiments and profiling, highlighting how precision scheduling and reduced inner multiplications yield notable performance gains in arbitrary-precision contexts. The framework is broadly applicable to other matrix polynomials and paves the way for further refinements, such as extending mixed-precision techniques to the computation of the first powers of and exploring more aggressive precision hierarchies.

Abstract

The Paterson--Stockmeyer method is an evaluation scheme for matrix polynomials with scalar coefficients that arise in many state-of-the-art algorithms based on polynomial or rational approximation, for example, those for computing transcendental matrix functions. We derive a mixed-precision version of the Paterson--Stockmeyer method that is particularly useful for evaluating matrix polynomials with scalar coefficients of decaying magnitude. The new method is mainly of interest in the arbitrary precision arithmetic, and it is attractive for high-precision computations. The key idea is to perform computations on data of small magnitude in low precision, and rounding error analysis is provided for the use of lower-than-the-working precisions. We focus on the evaluation of the Taylor approximants of the matrix exponential and show the applicability of our method to the existing scaling and squaring algorithms. We also demonstrate through experiments the general applicability of our method to other problems, such as computing the polynomials from the Padé approximant of the matrix exponential and the Taylor approximant of the matrix cosine. Numerical experiments show our mixed-precision Paterson--Stockmeyer algorithms can be more efficient than its fixed-precision counterpart while delivering the same level of accuracy.
Paper Structure (17 sections, 3 theorems, 65 equations, 3 figures, 2 tables)

This paper contains 17 sections, 3 theorems, 65 equations, 3 figures, 2 tables.

Table of Contents

  1. Introduction
  2. Rounding error analysis and evaluation scheme
  3. Powers and polynomials of matrices
  4. Taylor approximants to the matrix exponential
  5. Analysis with the norm constraint
  6. Applicability in the scaling and squaring algorithms for the matrix exponential
  7. The mixed-precision Paterson--Stockmeyer algorithm
  8. Numerical experiments
  9. Distributions of the used lower-than-working-precisions
  10. Accuracy and theoretical speed-up
  11. Code profiling and runtime comparison
  12. A mixed-precision Paterson--Stockmeyer algorithm for general polynomials of matrices
  13. Padé approximants to the matrix exponential
  14. Low-precision variant
  15. Taylor approximants to the matrix cosine
  16. ...and 2 more sections

Key Result

Theorem 1

\newlabelthm:ps-mp-err-norm0 If $\lVert\widehat{B}_i-B_i\rVert\le u_i\lVert B_i\rVert$, $i=\nu-1\colon r$ and $\lVert\widehat{Y}-Y\rVert\le u_{\nu}\lVert Y\rVert$ where $Y\equiv X^s$, then for the matrix $\widehat{q}$, computed in finite precision, for $q(X)$ in eq:ps-Horner-part, the evaluation s satisfies where

Figures (3)

  • Figure 1: Left: The relative error $\epsilon_v=\lVert\widehat{p}_m-{p}_m(X)\rVert_1/\lVert{p}_m(X)\rVert_1$ produced by Algorithm \ref{['alg.mpps']} compared with the relative error $\epsilon_f$ produced by the fixed-precision counterpart from fahi19 with $s=\lceil \sqrt{m} \rceil$ on various matrices with $2\le n\le 100$. Right: The associated approximate ratio of cost $C_r$ in \ref{['eq:Cr']}.
  • Figure 1: Left: The relative error $\epsilon_v=\lVert\widehat{p}_m-{p}_m(X)\rVert_1/\lVert{p}_m(X)\rVert_1$ produced by Algorithm \ref{['alg.mpps.gen']} in precision $u = 10^{-64}$ compared with the relative error $\epsilon_f$ produced by the fixed-precision counterpart from ahl22 or fahi19 with $s=\lceil \sqrt{m} \rceil$ on various matrices with $2\le n\le 100$. Right: The associated approximate ratio of cost $C_r$ in \ref{['eq:Cr']}.
  • Figure 2: The ratio between the runtime of Algorithm \ref{['alg.mpps']} and the runtime of its fixed-precision counterpart from fahi19 on matrices of different sizes, denoted by $\texttt{time[size]}$. The matrices are sorted in ascending order of their sparsity measured by the proportion of nonzero elements, characterized by the curves $\texttt{spar[size]}$.

Theorems & Definitions (4)

  • Theorem 1
  • Lemma 2
  • Proof 1
  • Theorem 3