Sylvester equations and polynomial separation of spectra

TL;DR

This paper explores conditions under which spectra of operators can be separated using polynomial transformations, providing new insights into the solvability of Sylvester equations in Banach spaces.

Contribution

It introduces criteria for polynomial separation of spectra of operators and offers an explicit series expansion for solutions based on these polynomials.

Findings

Spectra of operators can be separated using specific polynomials.

Explicit series expansion for Sylvester equation solutions is derived.

Conditions for spectral separation depend on the inclusion sets $V_p(T)$.

Abstract

Sylvester equations $AX-XB=C$ have unique solutions for all $C$ when the spectra of $A$ and $B$ are disjoint. Here $A$ and $B$ are bounded operators in Banach spaces. We discuss the existence of polynomials $p$ such that the spectra of $p(A)$ and $p(B)$ are well separated, either inside and outside of a circle or separated into different half planes. Much of the discussion is based on the following inclusion sets for the spectrum: $V_p(T)=\{\lambda \in \mathbb C \ : \ |p(\lambda)| \le \|p(T)\| \}$ where $T$ is a bounded operator. We also give an explicit series expansion for the solution in terms of $p(M)$, where $M=\begin{pmatrix} A&C\\ &B\end{pmatrix}$, in the case where the spectra of $A$ and $B$ lie in different components of $V_p(M)$ .

Equations

$$AX-XB=C$$

$$p(A)Y-Yp(B)=C$$

$$X=q(A,B)(Y).$$

$$p(A)X-Xp(B)=q(A,B)(C),$$

$$\phi (z)=\sum _{j=1}^d{\delta }_j(z)f_j(p(z))$$

$$\sigma (T)=\{\lambda \in \mathbb{C}:\lambda -T\text{ is not invertible}\}.$$

$$\sigma (A)\cap \sigma (B)=\varnothing ,$$

$$(A-\lambda )X-X(B-\lambda )=C$$

$$X=\frac{1}{2\pi i}{\int }_{\gamma }(\lambda -A{)}^{-1}C(\lambda -B{)}^{-1}d\lambda .$$

$$X=-\sum _{n=0}^{\infty }{A}^{n}C{B}^{-n-1}$$

$$\Vert A^n{\Vert }^{1/n}\Vert C{B}^{-1}{\Vert }^{1/n}\Vert B^{-n}{\Vert }^{1/n}\to \rho (A)\rho (B^{-1})<1.$$

$$(\lambda -A{)}^{-1}=\sum _{n=0}^{\infty }{\lambda }^{-n-1}{A}^n\text{ and }(\lambda -B{)}^{-1}=-\sum _{n=0}^{\infty }{\lambda }^n{B}^{-n-1}.$$

$$\rho (p(A))\rho (p(B^{-1}))<1$$

$$Y=-\sum _{n=0}^{\infty }p(A{)}^{n}Cp(B{)}^{-n-1}.$$

$$\hat{K}=\{z\in \mathbb{C}:|p(z)|\le \Vert p{\Vert }_K\text{ for all polynomials }p\}$$

$$\widehat{\sigma (A)}\cap \sigma (B)=\varnothing .$$

$$X={\int }_0^{\infty }{e}^{-tA}C{e}^{tB}dt.$$

$$\Vert e^{-tA}\hat{E}\Vert \le K{e}^{-\epsilon t}\text{ and }\Vert e^{tB}\hat{E}\Vert \le K{e}^{-\epsilon t}.$$

$$\sigma (p(A))\subset {\mathbb{C}}_{+}\text{ and }\sigma (p(B))\subset {\mathbb{C}}_{-}$$

$$\widehat{\sigma (A)}\cap \widehat{\sigma (B)}=\varnothing$$

$${\Sigma }_{\epsilon }(T)=\{\lambda \in \mathbb{C}:\text{ either }\lambda \in \sigma (T)\text{ or }\Vert (\lambda -T{)}^{-1}\Vert \ge \frac{1}{\epsilon }\}$$

$$\widehat{{\Sigma }_{\epsilon }(A)}\cap {\Sigma }_{\epsilon }(B)=\varnothing \text{ and }\widehat{{\Sigma }_{\epsilon }(A)}\cap \widehat{{\Sigma }_{\epsilon }(B)}=\varnothing$$

$$V_p(T)=\{\lambda \in \mathbb{C}:|p(\lambda )|\le \Vert p(T)\Vert \}$$

$$V_p(A)\cap \sigma (B)=\varnothing$$

$$\widehat{\sigma (A)}=\bigcap _{k\ge 1}{K}_k.$$

$$q(\lambda ,\mu )=\frac{p(\lambda )-p(\mu )}{\lambda -\mu }.$$

$$q(\lambda ,\mu )=\sum _{j=1}^d{\alpha }_j{q}_{j-1}(\lambda ,\mu ).$$

$$f(A,B):C\mapsto f(A,B)(C)=\sum _{i,j}{\alpha }_{ij}{A}^{i}C{B}^j.$$

$$g(A,B)(f(A,B)(C))=A^m(A^{i}C{B}^j)B^n=A^{m+i}C{B}^{n+j}=h(A,B)(C).$$

$$X=q(A,B)(Y)$$

Full text

*Aalto University

Department of Mathematics and Systems Analysis

email: Olavi.Nevanlinna@aalto.fi

Key words: Sylvester equation, multicentric calculus, preconditioning, spectral separation

2010 Mathematics Subject Classification: 15A24, 47A10, 47A60, 47A62, 65F08, 65F10, 65J10

1 Introduction

We discuss the solution of the Sylvester equation

[TABLE]

by solving first a related equation

[TABLE]

which is assumed to be easier to solve and then recover the solution of (1.1) as

[TABLE]

Here the operator $q(A,B)$ is obtained by the bivariate polynomial functional calculus from the divided difference of $p$, see Section 2, below. Alternatively, one can first form a new right hand side and consider solving

[TABLE]

see Propositions 2.2 and 2.4.

We consider the equations in the generality of bounded operators in Banach spaces. Given Banach spaces $\mathcal{X},\mathcal{Y}$ we assume that $A$ is bounded in $\mathcal{X}$, $B$ in $\mathcal{Y}$ and while $C$ and the unknowns $X$ and $Y$ are bounded operators from $\mathcal{Y}$ to $\mathcal{X}$. We discuss solution methods which can be formulated in infinite dimensional cases but which should be useful in matrix problems, in particular when the dimensions are large so that direct methods may not be practical. In this introduction we mention two basic representations for the solution, and then provide the spectral conditions under which a polynomial $p$ exist so that these methods can be used.

In a series of papers [10,11,12] we have studied the possibility of taking a polynomial as a new global variable. As polynomials are not injective we represent scalar functions $\varphi:z\mapsto\varphi(z)\in\mathbb{C}$ by vector valued functions $f:w\mapsto f(w)\in\mathbb{C}^{d}$ where $w=p(z)$ and $p$ is a polynomial of degree $d$ with simple roots $\lambda_{j}$. Then $\varphi$ is represented in the multicentric form

[TABLE]

where $\delta_{j}$ is the Lagrange polynomial $\delta_{j}(z)=\prod_{k\not=j}\frac{z-\lambda_{k}}{\lambda_{j}-\lambda_{k}}$. In this representation $\delta_{j}(A)$ is always well defined for any bounded operator and if $p(A)$ is "simpler" than $A$, small in norm, diagonalizable, normal, etc, an efficient functional calculus may be available for defining and computing $f_{j}(p(A))$.

Here the idea is again to replace the operators $A$ and $B$ by $p(A)$ and $p(B)$ but part of our dicussion is independent of the multicentric calculus. However, we discuss an application of the multicentric calculus which can be viewed as a modification of the *sign-*function approach, leading to a series expansion given in powers of $p(M)$ where $M=\begin{pmatrix}A&C\\ &B\end{pmatrix}$.

We shall now summarize the key results on the Sylvester equation, needed in the following. If $T$ is a bounded operator in a Banach space, then we denote by $\sigma(T)$ the spectrum:

[TABLE]

Bhatia and Rosenthal have written a readable survey of (1.1), [1]. They call the following as Sylvester-Rosenblum Theorem.

Theorem 1.1**.**

Let $\mathcal{X}$ and $\mathcal{Y}$ be Banach spaces and $A$, $B$ bounded operators in $\mathcal{X}$ and $\mathcal{Y}$, respectively. If

[TABLE]

then the equation (1.1) has a unique solution $X\in\mathcal{B}(\mathcal{Y},\mathcal{X})$ for every $C\in\mathcal{B}(\mathcal{Y},\mathcal{X})$.

We shall only consider the cases where (1.5) holds. Thus at least one of the operators $A$ and $B$ can be assumed to be nonsingular, and we shall assume that $B$ is. This is no restriction of generality as we could "transpose" the equation. Further, if $\lambda$ is a regular point for both $A$ and $B$ we could consider the equivalent equation

[TABLE]

instead and then both operators are invertible. This leads to the following representation of the solution.

Theorem 1.2**.**

([14]) If $\gamma$ is a union of closed contours with total winding numbers 1 around $\sigma(A)$ and 0 around $\sigma(B)$, then the solution of (1.1) can be expressed as

[TABLE]

Proof.

Operate (1.6) by $(\lambda-A)^{-1}$ from left and with $(\lambda-B)^{-1}$ from right. Integrating over $\gamma$ yields the claim.∎

Denote by $\rho(T)$ the spectral radius of $T$: $\rho(T)=\sup\{|\lambda|\ :\lambda\in\sigma(T)\}$.

Proposition 1.3**.**

Assume that $B$ is invertible and that $\rho(A)\rho(B^{-1})<1$. Then the series $\sum_{n=0}^{\infty}A^{n}CB^{-n-1}$ converges and setting

[TABLE]

we have a representation for the solution.

Proof.

The series converges as

[TABLE]

Multiplying the series by $A$ from left and subtracting the result of multiplying the series by $B$ from right then yields the claim.

Notice that this also follows from Theorem 1.2 since by assumption there exists an $r>0$ such that $\rho(A)<r$ and $\rho(B^{-1})<1/r$. Then we can integrate along $|\lambda|=r$ substituting

[TABLE]

∎

Our first aim is to discuss whether for given $A$ and $B$ there is a polynomial $p$ such that

[TABLE]

so that (1.2) could be solved as

[TABLE]

Recall, that the polynomially convex hull $\widehat{K}$ of a compact set $K\subset\mathbb{C}$ is defined as

[TABLE]

where $\|p\|_{K}=\sup_{z\in K}|p(z)|$. Thus $\widehat{K}$ is obtained by "filling the holes" of $K$. We have the following.

Theorem 1.4**.**

There exists a polynomial $p$ such that $p(B)$ is invertible and (1.9) holds if and only if

[TABLE]

The proof is in Section 3 where we also show how small the product in (1.9), when properly normalized, can be.

The second aim concerns another sufficient condition, based on the separation of the spectra of $A$ and $B$ by a vertical line. Again, by subtracting a suitable constant from the operators we may assume that the line is the imaginary axis. We shall denote by $\mathbb{C}_{+}$ the open right half plane and by $\mathbb{C}_{-}$ the open left half plane.

Theorem 1.5**.**

([5]) Suppose that the operators $A$, $B$ and $C$ are all bounded and that $\sigma(A)\subset\mathbb{C}_{+}$ and $\sigma(B)\subset\mathbb{C}_{-}$. Then the solution of (1.1) can be represented as

[TABLE]

Proof.

For a small enough $\varepsilon>0$ and large enough $K$ we have for $t>0$

[TABLE]

Thus, the integral converges and the claim follows by operating with $A$ from left and integrating by parts. ∎

Recall that under the assumptions of Theorem 1.5 the sign-function of the block operator $M$ is well defined and can be used to solve the Sylvester equation, see (4.2). On the possibility of separation into half planes we have the following result with proof in Section 4.

Theorem 1.6**.**

There exists a polynomial $p$ such that

[TABLE]

if and only if

[TABLE]

holds.

While (1.9) and (1.14) give the conditions under which these separating polynomials exist, one should expect that replacing the spectra by $\varepsilon$-pseudospectra should give useful information on the difficulty of computing these polynomials. Denoting the

[TABLE]

we could ask for how large $\varepsilon$ the conditions

[TABLE]

would hold. However, it seems that a more useful concept in this connection is the following inclusion set

[TABLE]

where $p$ is a polynomial. For (1.9) we would look for a polynomial $p$ such that

[TABLE]

while for (1.14) we would look for a polynomial such that $V_{p}(A\oplus B)$ separates into different components, containing $\sigma(A)$ and $\sigma(B)$, respectively.

In the practical search for separating polynomials, Krylov methods can be uselful, but one cannot in general guarantee that they would always produce separating polynomials when the necessary and sufficent spectral conditions hold. However, an idealized procedure exists with guaranteed performance. It assumes that one can perform minimizations of norms at polynomials of the operator and the key point is that one need not to know about the spectrum in advance. The following is Theorem 1.3 in [9], see also [4].

Theorem 1.7**.**

There exists a procedure which, given $A\in\mathcal{B}(\mathcal{X})$, produces a sequence of compact sets $K_{k}\subset\mathbb{C}$ and polynomials $p_{k}$ satisfying the following: $K_{k+1}\subset K_{k}$, $V_{p_{k}}(A)\subset K_{k}$, and

[TABLE]

In Section 2 we show how the post-processing is done. Sections 3 and 4 contain proofs of Theorems 1.4 and 1.6 and refinement of these.

At the end in Section 5 we take a somewhat different approach. We assume that we have a polynomial $p$ such that $V_{p}(M)$ separates into two components in which we define a piecewise constant holomorphic function. Using multicentric representation of this function we obtain a series expansion in terms of $p(M)$ from which the solution for the Sylvester equation can be read out in the same way as from sgn$(M)$. The coefficients of the series expansions can be computed with an explicit recursion depending on the polynomial $p$.

2 Post-processing

Assume that one has in one way or another solved the modified equation (1.2). We assume that we know the operators $A\in\mathcal{B}(\mathcal{X}),B\in\mathcal{B}(\mathcal{Y})$ and $Y\in\mathcal{B}(\mathcal{Y},\mathcal{X})$ and the (scalar) polynomial $p$. We shall use the bivariate polynomial calculus to write down the solution $X$ satisfying (1.1). To that end we associate with $p$ the bivariate polynomial $q$ as the divided difference of $p$:

[TABLE]

Denote $q_{k-1}(\lambda,\mu)=\lambda^{k-1}+\lambda^{k-2}\^{E}\mu+\cdots+\mu^{k-1}$ with $q_{0}=1$. Since $\lambda^{k}-\mu^{k}=(\lambda-\mu)q_{k-1}(\lambda,\mu)$ we then have with $p(\lambda)=\sum_{j=0}^{d}\alpha_{j}\lambda^{j}$

[TABLE]

On bivariate holomorphic functional calculus we recommend [7]. Since we deal here only with polynomials we can give the calculus without reference to integral representations. In the notation of [7], $q\{A,B^{T}\}(C)$ stands for our $q(A,B)(C)$.

Definition 2.1**.**

Let the operators $A\in\mathcal{B}(\mathcal{X}),B\in\mathcal{B}(\mathcal{Y})$ and $C\in\mathcal{B}(\mathcal{Y},\mathcal{X})$ and the polynomial $f(\lambda,\mu)=\sum_{i,j}\alpha_{ij}\lambda^{i}\mu^{j}$ be given. Then we denote by $f(A,B)$ the bounded linear operator in $\mathcal{B}(\mathcal{Y},\mathcal{X})$:

[TABLE]

When $f$ is holomorphic in two variables one defines $f(A,B)$ using a double integral and based on that one can prove that if $h(\lambda,\mu)=g(\lambda,\mu)f(\lambda,\mu)$ one gets $h(A,B)(C)=g(A,B)(f(A,B)(C))$. For polynomials this is obvious from (2.3) as we may work termwise. If $g(\lambda,\mu)=\lambda^{m}\mu^{n}$, $f(\lambda,\mu)=\lambda^{i}\mu^{j}$ then $g(\lambda,\mu)f(\lambda,\mu)=\lambda^{i+m}\mu^{j+n}=h(\lambda,\mu)$ and we have

[TABLE]

Taking linear combinations we see that $h(A,B)=g(A,B)\circ f(A,B)$ holds for polynomials $f$, $g$ where $h=gf$.

Consider now the post-processing step which is contained in the following simple result.

Proposition 2.2**.**

Let $A\in\mathcal{B}(\mathcal{X}),B\in\mathcal{B}(\mathcal{Y})$ and $Y\in\mathcal{B}(\mathcal{Y},\mathcal{X})$ be given and a polynomial $p$, such that (1.2) holds. Then

[TABLE]

satisfies the original Sylvester equation (1.1).

Proof.

We have

[TABLE]

Taking the left hand side as a polynomial of two variables and applying the polynomial functional calculus yields, by (1.2), $p(A)Y-Yp(B)=C$. Now the right hand side gives $Aq(A,B)(Y)-q(A,B)(Y)B=AX-XB$ , completing the proof.

∎

Example 2.3**.**

Let $A$ be a nonsingular real symmetric matrix, $B$ a real skew symmetric one. Then $A^{2}$ is positive definite while $B^{2}$ is negative semidefinite and

[TABLE]

solves the modified equation. Now $q(\lambda,\mu)=\lambda+\mu$ and we have the solution of the original Sylvester equation as

[TABLE]

The simple choice, $p(\lambda)=\lambda^{2}$ works naturally in a somewhat lager set of matrices. In fact, if there exists $\theta<1$ such that if $\alpha+i\beta\in\sigma(A)$ then $|\beta|\leq\theta|\alpha|$ while with $\gamma+i\delta\in\sigma(B)$ we ask for $|\gamma|\leq\theta|\delta|$. If at least one of $A$ or $B$ is nonsingular, then again the integral in (2.5) converges.

Denoting $S(\lambda,\mu)=\lambda-\mu$ the solution operator is the inverse of $S(A,B)$ satisfying

[TABLE]

Extending the bivarite polynomial calculus to holomorphic calculus one can show that if $f,g$ are holomorphic in two variables near the spectra and $h=gf$, then

[TABLE]

see e.g. Lemma 4.2 in [7]. Assuming this allows us to commute the terms in (2.6) and we conclude that rather than post-processing with $q(A,B)$ we may equally well begin with processing $C$. Clearly the order of computation is not the same but the operations needed to be excecuted essentially are. To summarise:

Proposition 2.4**.**

Let $A\in\mathcal{B}(\mathcal{X}),B\in\mathcal{B}(\mathcal{Y})$ and $C\in\mathcal{B}(\mathcal{Y},\mathcal{X})$ be given and a polynomial $p$ such that $\sigma(p(A))\cap\sigma(p(B))=\emptyset$. Then

[TABLE]

has a unique solution $X$ which also satisfies (1.1).

3 Disc separation

As before, $A\in\mathcal{B}(\mathcal{X})$ and $B\in\mathcal{B}(\mathcal{Y})$ and here we consider the convergence condition $\rho(p(A))\ \rho(p(B)^{-1})<1.$ Theorem 1.4 covers the existence of such polynomials and we give the proof here. We also derive an expression for the normalized infimum of the product of spectral radii. At the end of this section we discuss a more quantitative result.

If the spaces are finite dimensional, or more generally, if $A$ is an algebraic operator, then there exists a minimal polynomial $m_{A}$ such that $m_{A}(A)=0$, and assuming $\sigma(A)\cap\sigma(B)=\emptyset$, then trivially $\rho(m_{A}(A))\ \rho(m_{A}(B)^{-1})=0$. However, the degree of $m_{A}$ may be impractically high and computation of $m_{A}$ unstable.

Proof of Theorem 1.4

Suppose first that $\lambda_{0}\in\widehat{\sigma(A)}\cap\sigma(B)$ and let $p$ be a polynomial such that $p(B)$ is invertible. Then

[TABLE]

Since $\rho(p(A))\geq|p(\lambda_{0})|$ we have

[TABLE]

and we see that the condition (1.12) is necessary.

Assume then that (1.12) holds. As $\widehat{\sigma(A)}$ and $\sigma(B)$ are both compact, there exists an open $U$ such that $\widehat{\sigma(A)}\subset U$ while $\sigma(B)\cap U=\emptyset.$ By Hilbert Lemniscate Theorem, see e.g. Theorem 5.5.8 in [13], there exists a polynomial $p$ such that

[TABLE]

Thus, in particular

[TABLE]

and so $\rho(p(A))\ \rho(p(B)^{-1})<1,$ completing the proof. $\hfill\Box$

In practical computation, the spectral radius $\rho(p(A))$ should rather be replaced by $\|p(A)\|$ and scaled properly. To that end put

[TABLE]

where the infimum is over all polynomials $p$.

Lemma 3.1**.**

We have

[TABLE]

Proof.

The claim follows from the spectral radius formula. In fact, given $\varepsilon>0$ there exists a polynomial $q$ of degree $k$ such that

[TABLE]

But we have as $n\rightarrow\infty$

[TABLE]

so that $\eta(A,B)$ cannot be larger than $\inf(\rho(p(A))\ \rho(p(B)^{-1}))^{1/{deg(p)}}$. As it trivially cannot be smaller, (3.3) holds.

∎

It is of interest to know how small $\eta(A,B)$ can be. Given a polynomially convex compact set $K$ with positive logarithmic capacity, denote by $g$ the Green’s function of the complement of $K$, with singularity at $\infty$. That is, $g$ is harmonic in $\mathbb{C}\setminus K$,

[TABLE]

and such that for nearly everywhere on $\partial K$ $g(\zeta)\rightarrow 0$ as $\zeta$ tends to $\partial K$ from $\mathbb{C}\setminus K$, e.g. [13].

Theorem 3.2**.**

Assume (1.12) holds and $A$ is such that $\text{cap}(\widehat{\sigma(A)})>0.$ Denote by $g$ the Green’s function of $\mathbb{C}\setminus\widehat{\sigma(A)}$. Set $\alpha=\min_{\mu\in\sigma(B)}g(\mu)$. Then we have $0<\alpha<\infty$ and

[TABLE]

Proof.

Here we use Bernstein’s Lemma, as formulated in Theorem 5.5.7 of [13]. Since $\sigma(B)$ and $\widehat{\sigma(A)}$ are both compact, there is a positive distance between them and since $g$ is continuous and postive, we conclude $0<\alpha<\infty$. Then Bernstein’s Lemma yields for any polynomial $p$ of degree $d$

[TABLE]

which means

[TABLE]

Thus

[TABLE]

To get $\eta(A,B)$ bounded from above we use the following part of Theorem 5.5.7, [13]: if $p$ is a Fekete polynomial for $\widehat{\sigma(A)}$ of degree $d>1$, then

[TABLE]

Here $h$ is as follows:

[TABLE]

where $\tau$ is the Harnack distance for $\mathbb{C}\setminus\widehat{\sigma(A)}$. For us it suffices to know that $\tau$ is continuous and that

[TABLE]

Thus, for any $\varepsilon>0$ there exists a Fekete polynomial $p$ of degree $d$ such that

[TABLE]

But then

[TABLE]

Multiplying this with $\rho(p(A))^{1/d}$ gives

[TABLE]

which implies the bound from above. ∎

Recall, that operators $A\in\mathcal{B}(\mathcal{X})$ are called quasialgebraic if there exists a sequence $\{p_{j}\}$ of monic polynomials such that

[TABLE]

Halmos [3] has shown that a bounded operator is quasialgberaic if and only if the capacity of its spectrum vanishes. So, quasinilpotent, compact, polynomially compact, Riesz operators ect, are all quasialgebraic.

Theorem 3.3**.**

Let $A\in\mathcal{B}(\mathcal{X})$, $B\in\mathcal{B}(\mathcal{Y})$ satisfy $\widehat{\sigma(A)}\cap\sigma(B)=\emptyset$. Then $\eta(A,B)=0$ if and only if $A$ is quasialgebraic and, in particular $\widehat{\sigma(A)}=\sigma(A)$.

Proof.

That $A$ being quasialgebraic is necessary, follows immediately from Theorem 3.2. To obtain the other direction one needs to conclude that the superlinear decay quaranteed for $A$ can be obtained with a sequence of polynomials with roots staying away from the spectrum of $B$. This can be done for example by taking a nested sequence of compact sets $K_{n}$ such that $\cap K_{n}=\sigma(A)$, using Hilbert Lemniscate Theorem to get polynomials such that the associated lemniscates include $K_{n+1}$ but stay inside $K_{n}$. The related Green’s functions shall blow up at $\sigma(B)$.

∎

Remark 3.4**.**

In [8] we studied the polynomial acceleration speeds for the equation $x=Lx+f$ with $L$ a bounded operator in a Banach space $\mathcal{X}$. We formulated the equation in the fixed point form, rather than the usual $Ax=b$, to make the relationship between fixed point iteration and e.g. Krylov methods more apparent. Notice that viewing $x$ and $f$ as bounded operators $\mathbb{C}\rightarrow\mathcal{X}$, the fixed point equation can be viewed as a very special case of (1.1) with $A=L$ and $B=1$. The optimal asymptotic convergence rate is, in agreement with the results above,

[TABLE]

provided $1\notin\widehat{\sigma(A)}$, see Theorem 3.4.9 in [8]. Here $g$ denotes the Green’s function when the capacity is positive and can be thought as $+\infty$ when the capacity vanishes. We also discussed the superlinear behavior when the capacity vanishes and modelling the early behavior of iterations by assuming $1\in\partial{\widehat{\sigma(A)}}$ when the speed is sublinear.

We now derive a quantitative version of Theorem 1.4. Denote by $S(A,B)$ again the mapping $X\mapsto AX-XB$. Then the norm of $S(A,B)^{-1}$ can be used to bound the perturbation sensitivity. Since $S(A,B)^{-1}=q(A,B)\circ S(p(A),p(B))^{-1}$ we have

[TABLE]

When separating the operators using a polynomial $p$ the inversion should become easier but one would pay the prize of $q(A,B)$ typically having a large norm. However, as $q(A,B)$ is written out explicitly it can be thought of be applied exactly while the inversion part - when the dimensions are large or infinite - would typically be done only approximatively, e.g. by truncating an iteration.

It is tempting to replace the separation condition $\widehat{\sigma(A)}\cap\sigma(B)=\emptyset$ by the corresponding one on pseudospectra:

[TABLE]

in particular, as one of the the early applications of pseudospectrum was related to measuring the separation between matrices [15], [2]. However, we shall rather use the following condition

[TABLE]

which connects the polynomial $p$ directly into the estimates. In practice, one could calculate $\Sigma_{\varepsilon}(B)$ with moderate $\varepsilon$ and search for a polynomial $p$ e.g. by running an Arnoldi type Krylov process for a while and testing whether (3.8) is satisfied. This, or even the "ideal Arnodi" method, may not always produce polynomials with level set staying close to the spectrum. In fact, already the minimizing $\|p\|_{K}$ of monic polynomials of odd degree over $K=[-2,-1]\cup[1,2]$ necessarily has a zero at origin, staying far away from $K$. For that reason the process behind the proof of Theorem 1.7 is based on minimizing $\|p(A)\|$ over monic polynomials of given degree but includes a "cleaning" process - which most likely would not usually be needed. Notice also, that if $\widehat{\Sigma_{\varepsilon}(A)}$ is known and such that (3.7) holds, then one could compute Fekete points on $\widehat{\Sigma_{\varepsilon}(A)}$ to get a polynomial for which (3.8) could hold.

Assume now that $\varepsilon$ and $p$ are such that (3.8) holds. Then there exists $\delta>0$ and a contour $\gamma_{B}$ surrounding $\Sigma_{\varepsilon}(B)$, having vanishing total winding around $V_{p}(A)$, and such that along $\gamma_{B}$ we have $|p(\mu)|>\|p(A)\|+\delta$. Let $\ell_{B}$ be the length of $\gamma_{B}$. Then

[TABLE]

which implies

[TABLE]

so that

[TABLE]

Summing up we have the following.

Proposition 3.5**.**

Assume that there is a polynomial $p$ and $\varepsilon>0$ so that (3.8) holds. Then with $\delta$, $\ell_{B}$ as above we have

[TABLE]

Remark 3.6**.**

If $X=S(A,B)^{-1}(C)$ is wanted within some tolerance, notice that (3.6) and (3.10) allow one to calculate a safe truncation of the series expansion

[TABLE]

In fact, truncating

[TABLE]

and denoting $\widetilde{X}=q(A,B)(\widetilde{Y})$ we obtain $\|\widetilde{X}-X\|<tol$, providing $N$ is large enough so that

[TABLE]

holds, where $r=\|p(A)\|/(\|p(A)\|+\delta)$.

4 Half plane separation

The Theorem 1.6 deals with the question of existence of $p$ such that the spectra are separated into different half planes, allowing one to solve the modified equation using the integral representation (1.13) or the *sign-*function.

Observe that

[TABLE]

is satisfied exactly when $AX-XB=C$. If $\sigma(A)\subset\mathbb{C}_{+}$ and $\sigma(B)\subset\mathbb{C}_{-}$, the sign-function is well defined at $M$ and we have

[TABLE]

Thus, $X$ can be obtained if sgn$(M)$ can be computed. This is a rather popular route to compute the solution to Sylverster equation, see e.g. [1], [6].

We first prove the qualitative result of Theorem 1.6, then discuss how the lemniscate set $V_{p}(A\oplus B)$ can be used to obtain a quantitative result.

Proof of Theorem 1.6. The condition $\widehat{\sigma(A)}\cap\widehat{\sigma(B)}=\emptyset$ is necessary. In fact, assuming (1.14) holds, then we also have $\widehat{\sigma(p(A))}\subset\mathbb{C}_{+}$ and $\widehat{\sigma(p(B))}\subset\mathbb{C}_{-}$ and hence

[TABLE]

If $\lambda_{0}\in\widehat{\sigma(A)}\cap\widehat{\sigma(B)}$ we get a contradiction as

[TABLE]

Here the last step follows from the general fact that if $z\in\widehat{K}$ and q is any polynomial, then $|(q\circ p)(z)|\leq\|q\circ p\|_{K}=\|q\|_{p(K)}$ and so, $p(z)\in\widehat{p(K)}$.

Assume therefore that (1.15) holds and denote dist$(\sigma(A),\sigma(B))=\delta$. Put $U_{1}=\{\lambda\ :{\rm dist}(\lambda,\widehat{\sigma(A)})<\delta/3\}$ and $U_{2}=\{\mu\ :{\rm dist}(\mu,\widehat{(\sigma(B)})<\delta/3\}$. Then denote by $K$ the union of the closures of $U_{1}$ and $U_{2}$. Recall that $A(K)$ stands for continuous functions in $K$ which are holomorphic in the interior of $K$. Denote $c=$max$\{\|A\|,\|B\|\}+1$. Then we define a function $\varphi\in A(K)$ as follows

[TABLE]

Since $\mathbb{C}\setminus K$ is connected we may by Mergelyan’s Theorem approximate $\varphi$ arbitrarily accurately on $K$ by polynomials, say $\|\varphi-p\|_{K}<\varepsilon$. If $\gamma_{1}$ is a contour such that $\gamma_{1}$ surrounds $\widehat{\sigma(A)}$ inside $U_{1}$, then we have

[TABLE]

and in particular if $\varepsilon$ is small enough, $\sigma(p(A))\subset\mathbb{C}_{+}.$ Defining $\gamma_{2}$ in the similar way and integrating we get $p(B)$ with spectrum in the left half plane. $\hfill\Box$

We may replace the Mergelyan’s Theorem in the proof of Theorem 1.6 by the use of multicentric representation of $\varphi$. To that end, assume we have found polynomials $p_{1}$, $p_{2}$ such that

[TABLE]

e.g. based on Theorem 1.7. Let then $U_{i}$ be open, $V_{p_{1}}(A)\subset U_{1}$ and $V_{p_{2}}(B)\subset U_{2}$ and such that $\overline{U}_{1}\cap\overline{U}_{2}=\emptyset$. Then, again by Theorem 1.7, we may assume that, applied to the block diagonal operator $A\oplus B$, we have a polynomial $p$ such that

[TABLE]

Without loss of generality we may assume that $p$ is of degree $d$ and has simple roots $\lambda_{j}$. Let $t>0$ be small enough so that

[TABLE]

Define $\varphi$ on $\overline{U}_{1}\cup\overline{U}_{2}$ as in (4.3). We now use the multicentric representation (1.4) of $\varphi$ to approximate $\varphi(A)$ and $\varphi(B)$ by polynomials. When $|w|<|p(\lambda)|$ we have

[TABLE]

and the functions $f_{j}$ in

[TABLE]

satisfy

[TABLE]

see [10]. We put

[TABLE]

where we truncate the series expansion for the integral kernel after the index $N$

[TABLE]

so that in particular $P$ is a polynomial of degree $(N+1)d-1$ at most.

Let $\gamma=\gamma_{1}\cup\gamma_{2}$ with $\gamma_{i}\subset U_{i}$. The roots of $p$ are divided into two parts, say $\lambda_{j}\in U_{1}$ for $j\leq m$ and $\lambda_{k}\in U_{2}$ for $m<k\leq d$. Since $\sigma(A)\subset U_{1}$, the integral over $\gamma_{2}$ does not contribute into $\varphi(A)$ and we may estimate as follows. Denote

[TABLE]

Now

[TABLE]

and thus

[TABLE]

where we set $r=\frac{\|p(A)\|}{\|p(A\oplus B)\|+t}$. Likewise we obtain

[TABLE]

with $s=\frac{\|p(B)\|}{\|p(A\oplus B)\|+t}$. By the choice of $\varphi$ the spectrum of $\varphi(A)$ is in the half plane $Re\ \lambda>1$ while that of $\varphi(B)$ is likewise in the half plane $Re\ \mu<-1$. Choosing $N$ large enough so that

[TABLE]

we have $\sigma(P(A))\subset\mathbb{C}_{+}$ and $\sigma(P(B))\subset\mathbb{C}_{-}$.M

To summarize:

Proposition 4.1**.**

Assume that we have a polynomial $p$ such that (4.5) holds. Then we can estimate a truncation index $N$ such that

[TABLE]

holds with the polynomial $P$ in (4.6).

5 Explicit series expansion using multicentric calculus

In the previous section we demonstrated the existence polynomials for half plane separation. One could then compute the sign-function of

[TABLE]

and obtain the solution $X$ to the Sylvester equation from (4.2). This can be done for example using Newton’s iteration. We shall here *bypass * the mapping into different half planes. We use piecewise holomorphic functions to define the formal solution as a Cauchy-integral and then show how using multicentric calculus we get an explicit series expression for it. In the following we again assume all the time that $A\in\mathcal{B}(\mathcal{X})$, $B\in\mathcal{B}(\mathcal{Y})$ and $C\in\mathcal{B}(\mathcal{Y},\mathcal{X})$

Suppose we have open sets $U_{1},U_{2}$ such that $\widehat{\sigma(A)}\subset U_{1}$ and $\widehat{\sigma(B)}\subset U_{2}$ and $U_{1}\cap U_{2}=\emptyset$. Let $\varphi$ be the locally constant holomorphic function taking value 1 in $U_{1}$ and value $-1$ in $U_{2}$. If $\gamma_{1}$ is a contour inside $U_{1}$ surrounding $\sigma(A)$ we set

[TABLE]

Then the following holds.

Proposition 5.1**.**

In the notation above

[TABLE]

where $X$ is the solution of $AX-XB=C$.

Proof.

From (5.1) and

[TABLE]

we have

[TABLE]

∎

Our aim is now to compute $Q$. To that end let $\gamma=\gamma_{1}\cup\gamma_{2}$ where $\gamma_{2}$ is a contour surrounding $\sigma(B)$ inside $U_{2}$ so that, as $\gamma$ surrounds $\sigma(M)$, we have

[TABLE]

But then adding this to both sides of

[TABLE]

yields $\varphi(M)=2Q-I$ and $Q=\frac{1}{2}(\varphi(M)+I).$ Suppose we have a polynomial $p$ such that

[TABLE]

and $t>0$ small enough so that $\gamma=\{\lambda\ :|p(\lambda)|=\|p(M)\|+t\}\subset U_{1}\cup U_{2}$. Then $\gamma$ splitts into $\gamma_{1}$ and $\gamma_{2}$ in a natural way. We now write down the series expansion of $\varphi$ which converge inside $\gamma$, uniformly in compact subsets.

On the polynomial $p$ we assume that it has simple roots and is monic and of degree $d$. We write $\varphi$ in the multicentric form

[TABLE]

where the Taylor coefficients $\alpha_{j,k}$ in

[TABLE]

can be computed by an explicit recursion. The recursion is derived in [10]. Let $p$ have roots $\lambda_{j}$ and $\delta_{j}(\lambda)$ denote the polynomials taking value 1 at $\lambda_{j}$ and vanishing at the other roots. We may assume that $\lambda_{j}\in U_{1}$ for $j\leq s$ and $\lambda_{j}\in U_{2}$ for $s+1\leq j\leq d$. We first compute recursively polynomials $b_{n,m}$ as follows:

Put $b_{0,0}=1$, $b_{1,1}=p^{\prime}$ , $b_{n,0}=0$ for $n>0$ and for $m>n$ $b_{n,m}=0$. Then

[TABLE]

Then given the values $\varphi^{(n)}(\lambda_{j})$ we can compute $f_{j}^{(n)}(0)$ from the following

[TABLE]

This is Proposition 4.3 in [10]111where the last line (5.11) had dropped out. We can summarize:

Proposition 5.2**.**

Let $\varphi=1$ in $U_{1}$ and $\varphi=-1$ in $U_{2}$, and assume $p$ is such that (5.7) holds. Then we have $Q=\frac{1}{2}(\varphi(M)+I)$ where

[TABLE]

The Taylor coefficients of $f_{j}$ satisfy, see Proposition 4.4 in [10],

[TABLE]

Denote $L_{j}=\frac{1}{2\pi}\int_{\gamma}\frac{|d\lambda|}{|\lambda-\lambda_{j}|}$ then, we have $|\frac{f_{j}^{(n)}(0)}{n!}|\leq L_{j}(\|p(M)\|+t)^{-n},$ which allows us to truncate the series. Put

[TABLE]

so that

[TABLE]

where

[TABLE]

Let $tol>0$ be given and compute $N$ such that

[TABLE]

Proposition 5.3**.**

In the notation above, if $N$ is large enough so that (5.14) holds, then we have an approximation $\widetilde{X}$ to $X$ solving $AX-XB=C$ such that $\|\widetilde{X}-X\|<tol$, where $\widetilde{X}$ is the right upper corner element of $\widetilde{Q}=\frac{1}{2}(\widetilde{\varphi}(M)+I)$.

Remark 5.4**.**

We may assume without loss of generality that $p$ has simple rational roots, as conditions such as (5.7) allow small perturbations if needed. This means that the Taylor coefficients $\alpha_{j,n}$ are rational as well.

Remark 5.5**.**

Observe that we have an explicit formula for $p(M)^{k}$. In fact

[TABLE]

and so

[TABLE]

where $q_{k-1}(\lambda,\mu)=(\lambda^{k}-\mu^{k})/(\lambda-\mu)$.

Example 5.6**.**

We shall again demonstrate the approach using the special case as in Example 2.3. Let $A$ and $B$ be nonsigular bounded operators in a Hilbert space, such that $A$ and $iB$ are self adjoint, normalized e.g. so that both have norms bounded by 1. In particular then $A^{2}$ and $-B^{2}$ are both positive definite with spectra in some interval $[\alpha,1]$, with $\alpha>0$. We can proceed in two slightly different ways.

We could start by setting $\zeta=\lambda^{2}$ and solve

[TABLE]

using sign-function expansion in the polynomial $p(\zeta)=\zeta^{2}-1$. Or, you could solve

[TABLE]

with $p(\lambda)=\lambda^{4}-1$ so that $q(\lambda,\mu)=\lambda^{3}+\lambda^{2}\mu+\lambda\mu^{2}+\mu^{3}$. Here you should define $\varphi=1$ in the open sectors where $\arg(\lambda^{4})>0$ and $\varphi=-1$ where $\arg(\lambda^{4})<0$. Both approaches lead to an expansion in terms of powers of $M^{4}-I$ which is easy to derive directly. Consider the sign-function, defined for $Re\ \zeta\not=0$ as

[TABLE]

where $Re(\zeta^{2})^{1/2}>0$. With $w=\zeta^{2}-1$ and assuming that $|w|=|\zeta^{2}-1|<1$ we may expand $(1+w)^{-1/2}$ to get

[TABLE]

Now, we can apply this to the operator $M^{2}$. In fact, we have

[TABLE]

which converges as the spectral radius $\rho(M^{4}-I)=\|(A^{2}\oplus B^{2})^{2}-I\|<1$. The solution to the original equation is then the right upper corner element of $\frac{1}{2}{\rm sgn}(M^{2})$.

References

[1] R. Bhatia, P. Rosenthal, How and why to solve the operator equation AX - XB = Y, Bull. London Math. Soc., 29 (1997)1 - 21

[2] J. W. Demmel: The Condition Number of Equivalence Transformations that Block Diagonalize Matrix Pencils, SIAM J. Numer. Anal. 20, No. 3, (1983)

[3] P.R.Halmos, Capacity in Banach Algebras, IndianaUniv. Math., 20, (1971), pp 855-863,

[4] Anders C. Hansen, Olavi Nevanlinna, Complexity issues in computing spectra, pseudospectra and resolvents, Études Opératorielles, Banach Center Publ. Volume 112, (2017), pp 171 - 194, DOI: 10.4064/bc112-0-10

[5] E. Heinz, Beiträge zur Störungstheorie der Spektralzerlegung, Math. Ann. 123 (1951) 415 - 438.

[6] N. J. Higham. Functions of matrices. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, (2008)

[7] Daniel Kressner, Bivariate Matrix Functions, Operators and Matrices, 8,2, pp 449-466, (2014), DOI : 10.7153/oam-08-23

[8] Olavi Nevanlinna, Convergence of Iterations for Linear Equations, Birkhäuser, (1993)

[9] Olavi Nevanlinna, Computing the spectrum and representing the resolvent, Numer. Funct. Anal. Optim. 30 no.9 - 10 (2009), 1025 - 1047.

[10] O. Nevanlinna, Multicentric Holomorphic Calculus, Computational Methods and Function Theory, June 2012, Vol. 12, Issue 1, 45 - 65.

[11] O. Nevanlinna, Lemniscates and K -spectral sets, J. Funct. Anal. 262, (2012), 1728 - 1741.

[12] O. Nevanlinna, Polynomial as a new variable - a Banach algebra with functional calculus, Oper. and Matrices 10 (3) (2016) 567 - 592

[13] T. Ransford, Potential Theory in the Complex Plane, London Math. Soc. Student Texts 28, Cambridge Univ. Press, (1995)

[14] M. Rosenblum, On the operator equation BX-XA=Q, Duke Math.J. 23 (1956) 263 - 270

[15] J.M.Varah, On the Separation of Two Matrices, SIAM J. Numer. Anal., Vol. 16(2) pp. 216 - 222, (1979)