An estimate of approximation of a matrix-valued function by an interpolation polynomial

TL;DR

This paper provides an upper bound estimate for the approximation error when representing a matrix-valued function by an interpolation polynomial, involving derivatives and the spectrum of the matrix.

Contribution

It introduces a new matrix norm bound for the approximation error of matrix functions by interpolation polynomials, extending classical scalar polynomial approximation results.

Findings

Derived an explicit error bound involving matrix derivatives and spectrum.

Generalized scalar polynomial approximation bounds to matrix-valued functions.

Provides a theoretical tool for analyzing polynomial approximations of matrix functions.

Abstract

Let $A$ be a square complex matrix, $z_1$, ..., $z_{n}\in\mathbb C$ be (possibly repetitive) points of interpolation, $f$ be analytic in a neighborhood of the convex hull of the union of the spectrum of $A$ and the points $z_1$, ..., $z_{n}$, and $p$ be the interpolation polynomial of $f$, constructed by the points $z_1$, ..., $z_{n}$. It is proved that under these assumptions $$\Vert f(A)-p(A)\Vert\le\frac1{n!} \max_{t\in[0,1];\,\mu\in\text{co}\{z_1,z_{2},\dots,z_{n}\}}\bigl\Vert\Omega(A)f^{{(n)}} \bigl((1-t)\mu\mathbf1+tA\bigr)\bigr\Vert,$$ where $\Omega(z)=\prod_{k=1}^n(z-z_k)$.

Equations

$$\lVert f(A)-p(A)\rVert\leq\frac{1}{n!}\max_{\begin{subarray}{c}t\in[0,1]\\ \mu\in\operatorname{co}\{z_{1},z_{2},\dots,z_{n}\}\end{subarray}}\bigl{\lVert}\Omega(A)f^{{(n)}}\bigl{(}(1-t)\mu\mathbf{1}+tA\bigr{)}\bigr{\rVert},$$

$$|f(z)-p(z)|\le \frac{1}{n!}|\Omega (z)|\underset{\lambda \in \mathrm{co}\{z_1,z_2,\dots ,z_n,z\}}{\max }|f^{(n)}(\lambda )|,$$

$$\Omega (z)=\prod _{k=1}^n(z-z_k),$$

$$\lVert f(A)-p(A)\rVert\leq\frac{1}{n!}\max_{\begin{subarray}{c}t\in[0,1]\\ \mu\in\operatorname{co}\{z_{1},z_{2},\dots,z_{n}\}\end{subarray}}\bigl{\lVert}\Omega(A)f^{{(n)}}\bigl{(}(1-t)\mu\mathbf{1}+tA\bigr{)}\bigr{\rVert},$$

$$\Omega (z)=\prod _{k=1}^n(z-z_k).$$

$$f(z)-p(z)=\Omega (z)f[z_1,z_2,\dots ,z_n,z],$$

$$f[z_{1},z_{2},\dots,z_{n},z]=\int_{0}^{1}\int_{0}^{t_{1}}\dots\int_{0}^{t_{n-1}}f^{(n)}\bigl{(}z_{1}+(z_{2}-z_{1})t_{1}+\dots{\\ }+(z_{n}-z_{n-1})t_{n-1}+(z-z_{n})t_{n}\bigr{)}\,dt_{n}dt_{n-1}\dots dt_{1}.$$

$$f[z_{1},z_{2},\dots,z_{n},z]=\int_{0}^{1}\int_{0}^{t_{1}}\dots\int_{0}^{t_{n-1}}f^{{(n)}}\bigl{(}(1-t_{1})z_{1}+(t_{1}-t_{2})z_{2}+\dots{\\ }+(t_{n-1}-t_{n})z_{n}+t_{n}z\bigr{)}\,dt_{n}dt_{n-1}\dots dt_{1}.$$

$$(1-t_1)z_1+(t_1-t_2)z_2+\cdots +(t_{n-1}-t_n)z_n+t_{n}z$$

$$f(A)-p(A)=\Omega(A)\int_{0}^{1}\int_{0}^{t_{1}}\dots\int_{0}^{t_{n-1}}f^{{(n)}}\bigl{(}(1-t_{1})z_{1}\mathbf{1}+(t_{1}-t_{2})z_{2}\mathbf{1}+\dots{\\ }+(t_{n-1}-t_{n})z_{n}\mathbf{1}+t_{n}A\bigr{)}\,dt_{n}dt_{n-1}\dots dt_{1}.$$

$$\lVert f(A)-p(A)\rVert=\xi\bigl{(}f(A)-p(A)\bigr{)}.$$

$$\begin{split}\lVert f(A)-p(A)\rVert&=\xi\biggl{(}\int_{0}^{1}\int_{0}^{t_{1}}\dots\int_{0}^{t_{n-1}}\Omega(A)f^{{(n)}}\bigl{(}(1-t_{1})z_{1}\mathbf{1}+\dots\\ &\phantom{fghfg}+(t_{n-1}-t_{n})z_{n}\mathbf{1}+t_{n}A\bigr{)}\,dt_{n}dt_{n-1}\dots dt_{1}\biggr{)}\\ &\leq\biggl{\lVert}\int_{0}^{1}\int_{0}^{t_{1}}\dots\int_{0}^{t_{n-1}}\Omega(A)f^{{(n)}}\bigl{(}(1-t_{1})z_{1}\mathbf{1}+\dots\\ &\phantom{fghfg}+(t_{n-1}-t_{n})z_{n}\mathbf{1}+t_{n}A\bigr{)}\,dt_{n}dt_{n-1}\dots dt_{1}\biggr{\rVert}\\ &\leq\int_{0}^{1}\int_{0}^{t_{1}}\dots\int_{0}^{t_{n-1}}\bigl{\lVert}\Omega(A)f^{{(n)}}\bigl{(}(1-t_{1})z_{1}\mathbf{1}+\dots\\ &\phantom{fghfg}+(t_{n-1}-t_{n})z_{n}\mathbf{1}+t_{n}A\bigr{)}\bigr{\rVert}\,dt_{n}dt_{n-1}\dots dt_{1}\\ &\leq\int_{0}^{1}\int_{0}^{t_{1}}\dots\int_{0}^{t_{n-1}}\max_{0\leq t_{n}\leq\dots\leq t_{1}\leq 1}\bigl{\lVert}\Omega(A)f^{{(n)}}\bigl{(}(1-t_{1})z_{1}\mathbf{1}+\dots\\ &\phantom{fghfg}+(t_{n-1}-t_{n})z_{n}\mathbf{1}+t_{n}A\bigr{)}\bigr{\rVert}\,dt_{n}dt_{n-1}\dots dt_{1}.\end{split}$$

$$\frac{1}{1-t_{n}}\bigl{(}(1-t_{1})z_{1}+(t_{1}-t_{2})z_{2}+\dots+(t_{n-1}-t_{n})z_{n}\bigr{)}$$

$${\int }_0^1{\int }_0^{t_1}\cdots {\int }_0^{t_{n-1}}d{t}_n\dots d{t}_1=\frac{1}{n!}.$$

$$\lVert f(A)-p(A)\rVert\leq\frac{1}{n!}\max_{\begin{subarray}{c}t\in[0,1]\\ \mu\in\operatorname{co}\{z_{1},z_{2},\dots,z_{n}\}\end{subarray}}\bigl{\lVert}\Omega(A)f^{{(n)}}\bigl{(}(1-t)\mu\mathbf{1}+tA\bigr{)}\bigr{\rVert}.\qed$$

$$\max_{\mu\in\operatorname{co}\{z_{1},z_{2},\dots,z_{n}\}}\bigl{\lVert}\Omega(A)f^{{(n)}}\bigl{(}(1-t)\mu\mathbf{1}+tA\bigr{)}\bigr{\rVert}\\ =\max_{\mu\in\partial\operatorname{co}\{z_{1},z_{2},\dots,z_{n}\}}\bigl{\lVert}\Omega(A)f^{{(n)}}\bigl{(}(1-t)\mu\mathbf{1}+tA\bigr{)}\bigr{\rVert}.$$

$$\bigl{\lVert}\Omega(A)f^{{(n)}}\bigl{(}(1-t)\mu\mathbf{1}+tA\bigr{)}\bigr{\rVert}=\max_{\lVert\xi\rVert\leq 1}\xi\bigl{[}\Omega(A)f^{{(n)}}\bigl{(}(1-t)\mu\mathbf{1}+tA\bigr{)}\bigr{]},$$

$$\xi\mapsto\xi\bigl{[}\Omega(A)f^{{(n)}}\bigl{(}(1-t)\mu\mathbf{1}+tA\bigr{)}\bigr{]}$$

$$\max_{\mu\in\operatorname{co}\{z_{1},z_{2},\dots,z_{n}\}}\bigl{|}\xi\bigl{[}\Omega(A)f^{{(n)}}\bigl{(}(1-t)\mu\mathbf{1}+tA\bigr{)}\bigr{]}\bigr{|}\\ =\max_{\mu\in\partial\operatorname{co}\{z_{1},z_{2},\dots,z_{n}\}}\xi\bigl{[}\Omega(A)f^{{(n)}}\bigl{(}(1-t)\mu\mathbf{1}+tA\bigr{)}\bigr{]}.$$

$$f(\lambda )=\sum _{k=0}^{\infty }{c}_k(\lambda -z_1{)}^k,$$

$$\biggl{\|}f(A)-\sum_{k=0}^{n-1}c_{k}(A-z_{1}\mathbf{1})^{k}\biggr{\|}\leq\frac{1}{n!}\max_{t\in[0,1]}\bigl{\lVert}(A-z_{1}\mathbf{1})^{n}f^{(n)}\bigl{(}(1-t)z_{1}\mathbf{1}+tA\bigr{)}\bigr{\rVert}.$$

$$\alpha$$

$$\beta$$

$$\gamma$$

$$\Vert e^A-p(A)\Vert \le \frac{1}{n!}\underset{t\in [0,1]}{\max }{e}^{(1-t)\beta }\Vert \Omega (A)e^{tA}\Vert .$$

$$f^{{(n)}}\bigl{(}(1-t)\mu\mathbf{1}+tA\bigr{)}=e^{(1-t)\mu}e^{tA}.$$

$$\underset{\mu \in \mathrm{co}\{z_1,z_2,\dots ,z_n\}}{\max }|e^{(1-t)\mu }|=e^{(1-t)\beta }.\text{qed}$$

$$\Vert e^A-p(A){\Vert }_{2\to 2}\le e^{\gamma }\frac{\Vert \Omega (A){\Vert }_{2\to 2}}{n!}\sum _{j=0}^{n-1}\frac{(2\Vert A{\Vert }_{2\to 2}{)}^j}{j!}.$$

$$\Vert e^A-p(A){\Vert }_{2\to 2}\le \frac{\Vert \Omega (A){\Vert }_{2\to 2}}{n!}\underset{t\in [0,1]}{\max }{e}^{(1-t)\beta }\Vert e^{tA}{\Vert }_{2\to 2}.$$

$$\Vert e^{At}{\Vert }_{2\to 2}\le e^{\alpha t}\sum _{j=0}^{n-1}\frac{(2t\Vert A{\Vert }_{2\to 2}{)}^j}{j!},t\ge 0.\text{qed}$$

$$\Vert e^A-p(A){\Vert }_{2\to 2}\le e^{\gamma }\frac{\Vert \Omega (A){\Vert }_{2\to 2}}{n!}\sum _{j=0}^{n-1}\frac{\Vert N{\Vert }_{2\to 2}^j}{j!},$$

Full text

[TABLE]

where $\Omega(z)=\prod_{k=1}^{n}(z-z_{k})$ and the symbol $\operatorname{co}$ means the convex hull.

Introduction

An approximate calculation of analytic functions of matrices [4, 12] arises in many applications. One of the often used methods for the approximate calculation of a function $f$ of a large matrix $A$ is the replacement of $f$ by its polynomial approximation $p$. For the approximation by the Taylor polynomial, it is known a good estimate of accuracy [16], see Corollary 2. In this paper, we propose an estimate of the norm $\lVert f(A)-p(A)\rVert$, which is a generalization of the estimate from [16] for the case when $p$ is an interpolation polynomial of $f$. This estimate may help to choose an interpolation polynomial for an approximate calculation of a matrix function in an optimal way. Our estimate can be considered as a matrix analogue of the estimate [3, Theorem 3.1.1]

[TABLE]

where

[TABLE]

for the difference between an analytic function $f$ and its interpolation polynomial $p$ with respect to the points of interpolation $z_{1}$, …, $z_{n}$ provided that $f$ is analytic in a neighborhood of the convex hull of the points $z_{1}$, …, $z_{n}$ and $z$.

It is known many estimates of $\lVert f(A)\rVert$, see, e.g., [2, 7, 8, 9, 11, 15, 17, 18, 19]. All of them can be equivalently written as estimates of $\lVert f(A)-p(A)\rVert$, see, e.g., [10, Theorem 11.2.2]. The difference between these estimates and the proposed one (Theorem 1) is that the latter is adapted for the approximation by an interpolation polynomial.

In Section 1, we prove our estimate (Theorem 1) and describe some of its variants for the case of the matrix exponent. In Section 2, we give a numerical application.

1. The estimate

Theorem 1**.**

Let $A$ be a square complex matrix, $z_{1}$, …, $z_{n}\in\mathbb{C}$ be arbitrary (possibly repetitive) points of interpolation, $f$ be an analytic function defined in a neighborhood of the convex hull of the union of the spectrum of $A$ and the points $z_{1}$, …, $z_{n}$, and $p$ be the interpolation polynomial of $f$ constructed by the points $z_{1}$, …, $z_{n}$ (taking into account their multiplicities). Then (for any norm on the space of matrices)

[TABLE]

where $\mathbf{1}$ is the identity matrix, the symbol $\operatorname{co}$ means the convex hull, and

[TABLE]

Proof.

It is well-known (see, e.g., [3, Theorem 3.4.1] or [6, formula (52)]) that

[TABLE]

where $f[z_{1},z_{2},\dots,z_{n},z]$ is the divided difference [3, 6, 14]. On the other hand, by [6, formula (47)], we have

[TABLE]

Or

[TABLE]

Clearly, the complex numbers

[TABLE]

form the convex hull of the set $\{z_{1},z_{2},\dots,z_{n},z\}$ when $t_{1},\dots,t_{n}$ run through the set specified by the inequalities $0\leq t_{n}\leq\dots\leq t_{1}\leq 1$. Thus, considering integral (1), we use the fact that $f^{(n)}$ is defined on the convex hull of the points $z_{1}$, …, $z_{n}$, and $z$.

Substituting $A$ for $z$ into the previous formulas (thus, we assume that any point of the spectrum of $A$ can be taken as $z$, which can be done, since $f$ is analytic in a neighborhood of the convex hull of the union of the spectrum of $A$ and the points $z_{1}$, …, $z_{n}$), we obtain

[TABLE]

Let $\xi$ be a linear functional on the space of matrices (equipped by an arbitrary norm) such that $\lVert\xi\rVert=1$ and

[TABLE]

Such a functional exists by the Hahn-Banach theorem [13, Theorem 2.7.4]. Then we have the estimate

[TABLE]

We observe that the complex numbers

[TABLE]

form the convex hull of $\{z_{1},z_{2},\dots,z_{n}\}$ when $t_{1},\dots,t_{n}$ run through the set specified by the inequalities $0\leq t_{n}\leq\dots\leq t_{1}\leq 1$. Besides,

[TABLE]

Therefore from estimate (2) it follows that

[TABLE]

Remark 1*.*

For numerical calculations, it may be useful to note that the maximum can be taken over the boundary $\partial\operatorname{co}$ of a convex hull instead of the whole convex hull:

[TABLE]

Indeed, by the Hahn-Banach theorem,

[TABLE]

where the functional $\xi$ runs over the unit ball of the dual space of the space of all matrices. The function

[TABLE]

is analytic. Therefore, by the maximum modulus principle,

[TABLE]

Taking maximum over all functionals $\xi$ of unit norm, we arrive at the equality being proved.

Our Theorem 1 was inspired by the following result.

Corollary 2** ([16, Corollary 2], [12, Theorem 4.8]).**

Let the Taylor series

[TABLE]

where $c_{k}\in\mathbb{C}$, converges on an open circle of radius $r$ with the center at $z_{1}$, and the spectrum of a square matrix $A$ is contained in this circle. Then

[TABLE]

In corollaries below, we simplify the estimate from Theorem 1 for the case of the most important function $f(z)=e^{z}$.

In notation of Theorem 1, we set

[TABLE]

Corollary 3**.**

Let the assumptions of Theorem 1 be fulfilled and $f(z)=e^{z}$. Then

[TABLE]

Proof.

Clearly, $f^{(n)}(z)=e^{z}$. Therefore

[TABLE]

It remains to observe that

[TABLE]

The following three corollaries are more effective (but rougher) versions of the previous one.

We denote by $\lVert\cdot\rVert_{2\to 2}$ the matrix norm induced by the Euclidian norm on $\mathbb{C}^{n}$.

Corollary 4**.**

Let the assumptions of Theorem 1 be fulfilled and $f(z)=e^{z}$. Then

[TABLE]

where the matrix $A$ has the size $n\times n$.

Proof.

From Corollary 3 it follows that

[TABLE]

Next, we make use of the estimate [1, p. 131, Lemma 10.2.1], [5, p. 68, formula (13)]

[TABLE]

Corollary 5**.**

Let the assumptions of Theorem 1 be fulfilled and $f(z)=e^{z}$. Let the matrix $A$ be represented in the triangular Schur form [10] $A=Q^{-1}BQ$, where $B$ is triangular and $Q$ is unitary. Further, let $B=D+N$, where $D$ is diagonal and $N$ is strictly triangular. Then

[TABLE]

where the matrix $A$ has the size $n\times n$.

Proof.

The proof is similar to that of Corollary 4 and based on the estimate [18]

[TABLE]

Corollary 6**.**

Let the assumptions of Theorem 1 be fulfilled and $f(z)=e^{z}$. Let the matrix $A$ be normal. Then

[TABLE]

Proof.

For the normal matrix $A$, we have $N=0$. Therefore the proof follows from Corollary 5. ∎

Let us discuss whether the estimate is close to real accuracy.

Example 1*.*

Let the points of interpolation $z_{1}$, …, $z_{n}$ be taken coinciding with the points of the spectrum of $A$ (counted according to their algebraic multiplicities). Then $\Omega$ is the characteristic polynomial of a matrix $A$. By the Cayley–Hamilton theorem, $\Omega(A)=0$. Thus, in this case, Theorem 1 implies the well-known identity $p(A)=f(A)$. Similarly, if $\Omega$ and its derivatives are small on the spectrum of $A$, the factor $\Omega(A)$ is also small.

Example 2*.*

Let $A$ be a Hermitian matrix with the spectrum lying in $[-1,1]$. Let the points of interpolation be the zeroes of the Chebyshev polynomial of the first kind [3, § 3.3] of degree $n$ on $[-1,1]$. In this case, $\Omega$ is this Chebyshev polynomial; if its leading coefficient is 1, then the maximal absolute value of $\Omega$ on $[-1,1]$ is $\frac{1}{2^{n-1}}$. Therefore, Corollary 6 implies that

[TABLE]

If the spectrum of $A$ is not known exactly, the sharp estimate (for this polynomial) is

[TABLE]

We compare these two estimates for $n=10$: we have $\frac{1}{n!}\frac{e}{2^{n-1}}=1.46\cdot 10^{-9}$ and $\max_{\lambda\in[0,1]}|e^{\lambda}-p(\lambda)|=0.60\cdot 10^{-9}$. The comparison shows that the estimate from Corollary 6 is rather close to sharp one.

2. Numerical experiment

Theorem 1 can help to estimate whether the accuracy of the approximation of the matrix function $f(A)$ by a matrix polynomial $p(A)$ is good enough for given points of interpolation. We give an example of such a verification based on Corollary 3.

We put $N=1024$. We take complex numbers $d_{i}$, $i=1,\dots,N$, uniformly distributed in $[-1,0]\times[-i\pi,i\pi]$. We consider the diagonal matrix $D$ of the size $N\times N$ with the diagonal entries $d_{i}$. We take a matrix $T$, whose entries are random numbers uniformly distributed in $[-1,1]$. Then, we consider the matrix $A=TDT^{-1}$. Clearly, $\sigma(A)$ consists of the numbers $d_{i}$. We interpret $A$ as a random matrix whose spectrum is contained in the rectangle $[-1,0]\times[-i\pi,i\pi]$. In Fig. 1 we show an example of the spectrum of such a matrix.

For $t\in[0,1]$, we take as the sharp matrix $e^{tA}$ the matrix $Te^{tD}T^{-1}$.

We take the following 16 points as interpolation points:

[TABLE]

They are marked in Fig. 1 by the sign $\otimes$. These points are chosen heuristically. We calculate the interpolation polynomial $p$ and the polynomial $\Omega$, which correspond to these points, and substitute the matrix $A$ into them.

Next we calculate $\lVert\Omega(A)e^{tA}\rVert_{2\to 2}=\lVert\Omega(A)Te^{tD}T^{-1}\rVert_{2\to 2}$ for $t=0.01k$, where $k=0,\dots,100$, and take the maximum of these numbers as an approximate value of $\max_{t\in[0,1]}e^{(1-t)\beta}\lVert\Omega(A)e^{tA}\rVert$ (we put $\beta=0$). Finally, we divide the result by $16!$ and, thus, obtain the estimate from Corollary 3; we denote it by $e_{1}$. We also calculate the true accuracy $e_{0}=\lVert e^{A}-p(A)\rVert_{2\to 2}$ and the condition number $\varkappa(T)=\lVert T\rVert_{2\to 2}\cdot\lVert T^{-1}\rVert_{2\to 2}$.

We repeated the described experiment 100 times. After that, we excluded 3 results when $\varkappa(T)>10^{5}$. Finally, we calculated the average values. They are as follows: average $e_{0}$ is $9.36\cdot 10^{-6}$ with the standard deviation $1.37\cdot 10^{-5}$, average $e_{1}$ is $2.57\cdot 10^{-5}$ with the standard deviation $3.33\cdot 10^{-5}$, average $e_{1}/e_{0}$ is $3.03$ with the standard deviation $1.11$, average $\varkappa(T)$ is $1.84\cdot 10^{5}$ with the standard deviation $2.44\cdot 10^{5}$.

The average value $3.03$ of $e_{1}/e_{0}$ shows that the estimate is rather close to the true value. So, we can assume that for not very bad matrices $A$ of the size $N\times N$ with the spectrum in the rectangle $[-1,0]\times[-i\pi,i\pi]$, the interpolation polynomial with the considered interpolation points usually approaches $e^{A}$ with accuracy about $3\cdot 10^{-5}$.

Acknowledgements

The first author was supported by the Ministry of Education and Science of the Russian Federation under state order No. 3.1761.2017/4.6. The second author was supported by the Russian Foundation for Basic Research under research project No. 19-01-00732 .