---
title: Partial Fractions with Matrix Coefficients
url: https://www.emergentmind.com/topics/partial-fractions-with-matrix-coefficients
type: topic
---

# Partial Fractions with Matrix Coefficients

Searching arXiv for recent and foundational papers on partial fractions with matrix coefficients and related matrix-function PFEs.
Partial fractions with matrix coefficients are decompositions of rational matrix functions—most notably the resolvent $(sI-A)^{-1}$—into sums of pole terms with constant matrix multipliers. In the resolvent case, if $A \in \mathbb{C}^{n\times n}$ has eigenvalues $\lambda_1,\dots,\lambda_s$, the method writes
$$
(sI-A)^{-1}=\sum_{i=1}^s\sum_{j=1}^{r_i} B_{ij}(s-\lambda_i)^{-j},
$$
with uniquely determined matrices $B_{ij}\in\mathbb{C}^{n\times n}$. The method is developed as a unified tool for finding chains of generalized eigenvectors, evaluating matrix exponentials, and solving linear systems of ordinary differential equations with constant coefficients; it also connects naturally to spectral projectors, nilpotent parts, and Laplace-transform methods [2507.09661].

## 1. Formal framework and spectral data

Let $A_k\in\mathbb{C}^{n\times n}$. A matrix polynomial is
$$
P(\lambda)=\sum_{k=0}^m A_k\lambda^k.
$$
A rational matrix function is any function $R(\lambda)$ of the form
$$
R(\lambda)=P(\lambda)Q(\lambda)^{-1},
$$
where $P(\lambda)$ is a matrix polynomial and $Q(\lambda)$ is usually a scalar polynomial or, more generally, a matrix polynomial if specified. The central object is the resolvent
$$
(\lambda I-A)^{-1},
$$
whose scalar denominator is $\det(\lambda I-A)$ and whose numerator is $\operatorname{adj}(\lambda I-A)$.

If $\lambda_1,\dots,\lambda_s$ are the eigenvalues of $A$ and $r_j$ denotes the algebraic multiplicity of $\lambda_j$, then
$$
\det(\lambda I-A)=\prod_{j=1}^s (\lambda-\lambda_j)^{r_j}.
$$
The poles of the resolvent occur at the eigenvalues, with orders equal to the sizes $m_j$ of the largest Jordan blocks for $\lambda_j$, equivalently the multiplicities in the minimal polynomial. For each eigenvalue, let $P_j$ denote the spectral projector onto the generalized eigenspace and define
$$
N_j=(A-\lambda_j I)P_j.
$$
Then $N_j$ is nilpotent on that generalized eigenspace, with $N_j^{m_j}=0$.

A Jordan chain of length $p$ for $\lambda_j$ is a sequence $v^{(1)},\dots,v^{(p)}$ satisfying
$$
(A-\lambda_j I)v^{(1)}=0,\qquad (A-\lambda_j I)v^{(k+1)}=v^{(k)}.
$$
The space spanned by all such chains is the generalized eigenspace onto which $P_j$ projects. In this setting, the matrix coefficients in the partial fraction expansion are not auxiliary parameters; they encode spectral projectors, nilpotent structure, and generalized eigenvector data.

## 2. Decomposition theorem and matrix residues

For $A\in\mathbb{C}^{n\times n}$ with eigenvalues $\lambda_1,\dots,\lambda_s$ and algebraic multiplicities $r_1,\dots,r_s$, the resolvent admits a unique partial fraction decomposition with matrix coefficients,
$$
(sI-A)^{-1}=\sum_{i=1}^s\sum_{j=1}^{r_i} B_{ij}(s-\lambda_i)^{-j},
$$
where the matrices $B_{ij}$ are constant. There is no polynomial regular part in the resolvent case because $\operatorname{adj}(sI-A)$ has degree $n-1$ while $\det(sI-A)$ has degree $n$, so the rational matrix function is proper.

Existence and uniqueness follow entrywise from scalar partial fraction decomposition. Equivalently, one may write
$$
(sI-A)^{-1}=\frac{\operatorname{adj}(sI-A)}{\det(sI-A)}
$$
and expand each entry against the factorization of the scalar denominator. No diagonalizability is required. Distinct eigenvalues produce only simple poles, whereas defective matrices produce higher-order poles, up to the size of the largest Jordan block for the corresponding eigenvalue.

A constructive procedure begins from the ansatz above, multiplies by the common denominator, and obtains the polynomial identity
$$
\sum_{i=1}^s\sum_{j=1}^{r_i}(s-\lambda_i)^{r_i-j}(sI-A)B_{ij}
=
\prod_{i=1}^s (s-\lambda_i)^{r_i}I.
$$
Evaluation at $s=\lambda_i$ and differentiation up to order $r_i-1$ yield linear equations for the unknown matrices. In particular,
$$
(\lambda_i I-A)B_{i,r_i}=0,
$$
so the columns of $B_{i,r_i}$ are eigenvectors for $\lambda_i$, and
$$
(A-\lambda_i I)B_{i,j}=B_{i,j+1},\qquad j=1,\dots,r_i-1.
$$
Hence
$$
B_{i,j}=(A-\lambda_i I)^{j-1}B_{i,1}.
$$

For repeated poles, the coefficients may also be written by the usual repeated-pole residue formula applied entrywise:
$$
B_{i,j}=\frac{1}{(j-1)!}\lim_{s\to\lambda_i}\frac{d^{\,j-1}}{ds^{\,j-1}}\Big[(s-\lambda_i)^j(sI-A)^{-1}\Big].
$$
The method therefore mirrors scalar Heaviside-style partial fractions, but with matrix unknowns rather than scalar coefficients.

## 3. Spectral projectors, nilpotent parts, and generalized eigenvector chains

The coefficient $B_{i,1}$ has a distinguished spectral meaning. If $\Gamma_i$ is a simple positively oriented contour enclosing $\lambda_i$ and no other eigenvalues, then the spectral projector is
$$
P_i=\frac{1}{2\pi i}\oint_{\Gamma_i}(sI-A)^{-1}\,ds.
$$
Substituting the partial fraction expansion and integrating term-by-term gives
$$
P_i=B_{i,1}.
$$
Thus the first matrix coefficient is exactly the spectral projector onto the generalized eigenspace of $\lambda_i$.

Near a pole, the resolvent has the nilpotent expansion
$$
(sI-A)^{-1}
=
\sum_{k=0}^{m_i-1}(N_i^kP_i)(s-\lambda_i)^{-(k+1)},
$$
where $N_i=(A-\lambda_i I)P_i$ and $m_i$ is the size of the largest Jordan block for $\lambda_i$. Comparison with the partial fraction coefficients yields
$$
B_{i,1}=P_i,\qquad B_{i,2}=N_iP_i,\qquad \dots,\qquad B_{i,m_i}=N_i^{m_i-1}P_i.
$$
This identifies the matrix residues with projector and nilpotent data without explicit recourse to Jordan canonical form.

The coefficients also encode generalized eigenvectors. Any nonzero column of $B_{i,r_i}$ is an eigenvector of $A$ for $\lambda_i$. Any nonzero column of $B_{i,j}$ is a generalized eigenvector of rank at most $r_i+1-j$. For fixed $i$ and a fixed column index $m$, the nonzero $m$-th columns of
$$
B_{i,1},B_{i,2},\dots,B_{i,r_i}
$$
form a Jordan chain ending in an eigenvector. If the last nonzero column occurs at index $l_m$, then
$$
(A-\lambda_i I)v_{l_m}=0,\qquad
(A-\lambda_i I)v_{l_m-1}=v_{l_m},\qquad \dots,\qquad
(A-\lambda_i I)v_1=v_2.
$$

A common misconception is that partial fraction methods for the resolvent require diagonalizability or explicit Jordan reduction. The method does not require diagonalizability, and the coefficient-matching construction avoids explicit Jordan form. The paper states that this is often computationally simpler and more stable for moderate-size problems, while still reproducing the same spectral information.

## 4. Matrix exponentials and linear ODE systems

The matrix exponential follows directly from the Bromwich integral,
$$
e^{At}
=
\frac{1}{2\pi i}\int_{\gamma-i\infty}^{\gamma+i\infty} e^{st}(sI-A)^{-1}\,ds,
$$
for $\gamma$ larger than the spectral abscissa of $A$. Substituting the partial fraction expansion and using
$$
\mathcal{L}^{-1}\{(s-\lambda)^{-j}\}
=
e^{\lambda t}\frac{t^{j-1}}{(j-1)!},
$$
one obtains
$$
e^{At}
=
\sum_{i=1}^s e^{\lambda_i t}\sum_{j=1}^{r_i}\frac{t^{j-1}}{(j-1)!}B_{ij}.
$$
In terms of the nilpotent data,
$$
e^{At}
=
\sum_{i=1}^s e^{\lambda_i t}\sum_{k=0}^{m_i-1}\frac{t^k}{k!}N_i^kP_i.
$$
This is the standard Jordan-block expansion, recovered from the partial fractions of the resolvent.

Over the real field, if $\det(sI-A)$ factors into linear and quadratic factors, the decomposition may be written in real partial fractions. For a factor $s^2+a^2$, terms of the form
$$
\frac{s}{s^2+a^2}B_1+\frac{1}{s^2+a^2}B_2
$$
invert to
$$
\cos(at)B_1+\frac{1}{a}\sin(at)B_2.
$$
The method therefore accommodates complex conjugate spectral pairs without leaving real arithmetic.

For the linear system
$$
x'(t)=Ax(t)+b(t),
$$
the homogeneous solution is
$$
x_h(t)=e^{At}x(0).
$$
Variation of parameters gives
$$
x(t)=e^{At}x(0)+\int_0^t e^{A(t-\tau)}b(\tau)\,d\tau.
$$
If $b(t)$ is a combination of exponentials and polynomials, or of sinusoids, the convolution can be evaluated explicitly because $e^{A(t-\tau)}$ decomposes into sums of $e^{\lambda_i(t-\tau)}$ times polynomials in $(t-\tau)$. In Laplace form,
$$
X(s)=(sI-A)^{-1}x(0)+(sI-A)^{-1}B(s),
$$
so the transfer function from input $b$ to state $x$ is $(sI-A)^{-1}$, and its poles and matrix residues determine the explicit time-domain response.

## 5. Computation, worked examples, and numerical issues

A practical computation proceeds without Jordan form. The characteristic polynomial is first factored as
$$
\det(sI-A)=\prod_{i=1}^s (s-\lambda_i)^{r_i},
$$
or, over $\mathbb{R}$, complex conjugate pairs are grouped into irreducible quadratics. One then writes a partial fraction ansatz with undetermined matrices, multiplies by the common denominator to obtain a polynomial identity, and determines the coefficients by substitution at poles, differentiation at repeated poles, or coefficient matching. If $\operatorname{adj}(sI-A)=M(s)$ is computed symbolically, one may also write
$$
(sI-A)^{-1}=\frac{M(s)}{\det(sI-A)}
$$
and solve for the matrix coefficients by polynomial matching.

Once the coefficients are known, the exponential is obtained termwise from inverse Laplace transforms. The same coefficients then provide generalized eigenvectors: columns of $B_{i,r_i}$ are eigenvectors, while columns of $B_{i,j}$ are generalized eigenvectors linked by
$$
(A-\lambda_i I)B_{i,j}=B_{i,j+1}.
$$

The numerical limitations are explicit. Factoring $\det(sI-A)$ is sensitive to roundoff when eigenvalues are clustered or repeated; symbolic or high-precision arithmetic is recommended. Matching coefficients and differentiation at repeated poles amplify conditioning issues; contour-integral-based projectors $P_i$ may be more stable numerically. For large $n$, computing $\operatorname{adj}(sI-A)$ is expensive, so the substitution and matching procedure is preferred. The paper does not provide error bounds, and standard numerical linear algebra considerations apply.

Two examples illustrate the method. In the diagonalizable $2\times2$ case,
$$
A=
\begin{bmatrix}
6 & 4\\
-3 & -1
\end{bmatrix},
$$
one has
$$
\det(sI-A)=(s-2)(s-3),
$$
and
$$
(sI-A)^{-1}
=
\frac{1}{s-2}B_1+\frac{1}{s-3}B_2.
$$
From substitution at $s=2$ and $s=3$,
$$
-B_1=
\begin{bmatrix}
3 & 4\\
-3 & -4
\end{bmatrix},
\qquad
B_2=
\begin{bmatrix}
4 & 4\\
-3 & -3
\end{bmatrix}.
$$
Therefore
$$
e^{tA}
=
e^{2t}
\begin{bmatrix}
-3 & -4\\
3 & 4
\end{bmatrix}
+
e^{3t}
\begin{bmatrix}
4 & 4\\
-3 & -3
\end{bmatrix}.
$$
Because the eigenvalues are distinct, only simple poles occur and there are no polynomial-in-$t$ factors.

In the defective $3\times3$ example,
$$
A=
\begin{bmatrix}
0 & 1 & 2\\
-2 & 4 & 0\\
-1 & 1 & 2
\end{bmatrix},
$$
the characteristic polynomial is
$$
\det(sI-A)=(s-2)^3,
$$
so the only eigenvalue is $2$, with algebraic multiplicity $3$ and a single Jordan chain of length $3$. The decomposition
$$
(sI-A)^{-1}
=
\frac{1}{s-2}B_1+\frac{1}{(s-2)^2}B_2+\frac{1}{(s-2)^3}B_3
$$
has
$$
B_1=I,
$$
$$
B_2=
\begin{bmatrix}
-2 & 1 & 2\\
-2 & 2 & 0\\
-1 & 1 & 0
\end{bmatrix},
\qquad
B_3=
\begin{bmatrix}
0 & 2 & -4\\
0 & 2 & -4\\
0 & 1 & -2
\end{bmatrix}.
$$
Hence
$$
e^{tA}=e^{2t}\left(B_1+tB_2+\frac{t^2}{2}B_3\right),
$$
and the generalized eigenvector chain is read from
$$
(A-2I)B_1=B_2,\qquad (A-2I)B_2=B_3,\qquad (A-2I)B_3=0.
$$

## 6. Relation to rational matrix-function approximation and large sparse computation

The partial fraction method for the resolvent belongs to a broader family of techniques for evaluating matrix functions through rational approximations. For large, sparse, and/or localized matrices, one studies a scalar analytic function $\Psi$ and approximates $\Psi(A)$ or $\Psi(A)\mathbf v$ by a partial fraction expansion such as
$$
f(A)=\sum_{j=1}^N c_j(\xi_j I-A)^{-1},
$$
or, more generally,
$$
r_m(A)=c_0I+\sum_{j=1}^m \alpha_j(A-\sigma_j I)^{-1}.
$$
In this setting, the computation is reduced to families of shifted linear systems that share the same matrix structure, and the shifts can be parallelized [1709.06351].

The conceptual relation is direct. In the sparse-matrix literature, the coefficients are usually scalar residues multiplying matrix resolvents. When $\Psi(A)$ is expressed through spectral or Jordan decomposition, however, projectors and nilpotent blocks play the role of matrix residues, and repeated poles correspond to higher-order resolvents $(A-\lambda_k I)^{-\ell}$ with coefficients involving derivatives of $\Psi$. This suggests that Airapetyan’s matrix-coefficient decomposition of the resolvent and large-scale rational approximations of general matrix functions are two instances of a common rational-calculus viewpoint.

For large-scale problems, the computational emphasis shifts from symbolic coefficient matching to efficient shifted solves and approximate inverses. The cited work studies seed approximate inverse factorizations, sparse updates across shifts, off-diagonal decay in $A^{-1}$, and an error decomposition into rational approximation error and solve or preconditioner error. It reports applications to $\log(A)$, $\exp(A)$, fractional powers, localized matrices, PDE discretizations, Markov generator matrices, and large power systems matrices. The approach is advantageous when $A$ is sparse or localized, $\Psi$ is analytic on a region enclosing the spectrum, multiple right-hand sides are present, and moderate accuracy suffices. Its limitations include non-normality, poorly conditioned shifts, and failure of decay in $A^{-1}$.

The contrast between the two settings is therefore one of scale and objective rather than of principle. For moderate-size matrices, matrix-coefficient partial fractions expose spectral projectors, generalized eigenvectors, and closed-form expressions for $e^{At}$ and ODE solutions. For large sparse matrices, partial fraction expansions serve as a computational reduction to structured shifted systems. In both regimes, the central object is the pole structure of the resolvent and the algebra encoded by its residues.

Source: https://www.emergentmind.com/topics/partial-fractions-with-matrix-coefficients