Partial Fractions with Matrix Coefficients

• Partial fractions with matrix coefficients decompose the resolvent (sI-A)⁻¹ into pole terms that encode spectral projectors, nilpotency, and generalized eigenvector chains.
• The method provides a unified framework for evaluating matrix exponentials and solving ODE systems by linking matrix residues to eigen-structure without requiring explicit Jordan reduction.
• This approach offers computational advantages for moderate-size matrices and is applicable to both real and complex systems, balancing numerical stability and analytic insight.

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 (Airapetyan, 13 Jul 2025).

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 $A \in \mathbb{C}^{n\times n}$0 is usually a scalar polynomial or, more generally, a matrix polynomial if specified. The central object is the resolvent

$A \in \mathbb{C}^{n\times n}$1

whose scalar denominator is $A \in \mathbb{C}^{n\times n}$2 and whose numerator is $A \in \mathbb{C}^{n\times n}$3.

If $A \in \mathbb{C}^{n\times n}$4 are the eigenvalues of $A \in \mathbb{C}^{n\times n}$5 and $A \in \mathbb{C}^{n\times n}$6 denotes the algebraic multiplicity of $A \in \mathbb{C}^{n\times n}$7, then

$A \in \mathbb{C}^{n\times n}$8

The poles of the resolvent occur at the eigenvalues, with orders equal to the sizes $A \in \mathbb{C}^{n\times n}$9 of the largest Jordan blocks for $\lambda_1,\dots,\lambda_s$0, equivalently the multiplicities in the minimal polynomial. For each eigenvalue, let $\lambda_1,\dots,\lambda_s$1 denote the spectral projector onto the generalized eigenspace and define

$\lambda_1,\dots,\lambda_s$2

Then $\lambda_1,\dots,\lambda_s$3 is nilpotent on that generalized eigenspace, with $\lambda_1,\dots,\lambda_s$4.

A Jordan chain of length $\lambda_1,\dots,\lambda_s$5 for $\lambda_1,\dots,\lambda_s$6 is a sequence $\lambda_1,\dots,\lambda_s$7 satisfying

$\lambda_1,\dots,\lambda_s$8

The space spanned by all such chains is the generalized eigenspace onto which $\lambda_1,\dots,\lambda_s$9 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 $$(sI-A)^{-1}=\sum_{i=1}^s\sum_{j=1}^{r_i} B_{ij}(s-\lambda_i)^{-j},$$0 with eigenvalues $$(sI-A)^{-1}=\sum_{i=1}^s\sum_{j=1}^{r_i} B_{ij}(s-\lambda_i)^{-j},$$1 and algebraic multiplicities $$(sI-A)^{-1}=\sum_{i=1}^s\sum_{j=1}^{r_i} B_{ij}(s-\lambda_i)^{-j},$$2, 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},$$3

where the matrices $$(sI-A)^{-1}=\sum_{i=1}^s\sum_{j=1}^{r_i} B_{ij}(s-\lambda_i)^{-j},$$4 are constant. There is no polynomial regular part in the resolvent case because $$(sI-A)^{-1}=\sum_{i=1}^s\sum_{j=1}^{r_i} B_{ij}(s-\lambda_i)^{-j},$$5 has degree $$(sI-A)^{-1}=\sum_{i=1}^s\sum_{j=1}^{r_i} B_{ij}(s-\lambda_i)^{-j},$$6 while $$(sI-A)^{-1}=\sum_{i=1}^s\sum_{j=1}^{r_i} B_{ij}(s-\lambda_i)^{-j},$$7 has degree $$(sI-A)^{-1}=\sum_{i=1}^s\sum_{j=1}^{r_i} B_{ij}(s-\lambda_i)^{-j},$$8, so the rational matrix function is proper.

Existence and uniqueness follow entrywise from scalar partial fraction decomposition. Equivalently, one may write

$$(sI-A)^{-1}=\sum_{i=1}^s\sum_{j=1}^{r_i} B_{ij}(s-\lambda_i)^{-j},$$9

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

$B_{ij}\in\mathbb{C}^{n\times n}$0

Evaluation at $B_{ij}\in\mathbb{C}^{n\times n}$1 and differentiation up to order $B_{ij}\in\mathbb{C}^{n\times n}$2 yield linear equations for the unknown matrices. In particular,

$B_{ij}\in\mathbb{C}^{n\times n}$3

so the columns of $B_{ij}\in\mathbb{C}^{n\times n}$4 are eigenvectors for $B_{ij}\in\mathbb{C}^{n\times n}$5, and

$B_{ij}\in\mathbb{C}^{n\times n}$6

Hence

$B_{ij}\in\mathbb{C}^{n\times n}$7

For repeated poles, the coefficients may also be written by the usual repeated-pole residue formula applied entrywise:

$B_{ij}\in\mathbb{C}^{n\times n}$8

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_{ij}\in\mathbb{C}^{n\times n}$9 has a distinguished spectral meaning. If $A_k\in\mathbb{C}^{n\times n}$0 is a simple positively oriented contour enclosing $A_k\in\mathbb{C}^{n\times n}$1 and no other eigenvalues, then the spectral projector is

$A_k\in\mathbb{C}^{n\times n}$2

Substituting the partial fraction expansion and integrating term-by-term gives

$A_k\in\mathbb{C}^{n\times n}$3

Thus the first matrix coefficient is exactly the spectral projector onto the generalized eigenspace of $A_k\in\mathbb{C}^{n\times n}$4.

Near a pole, the resolvent has the nilpotent expansion

$A_k\in\mathbb{C}^{n\times n}$5

where $A_k\in\mathbb{C}^{n\times n}$6 and $A_k\in\mathbb{C}^{n\times n}$7 is the size of the largest Jordan block for $A_k\in\mathbb{C}^{n\times n}$8. Comparison with the partial fraction coefficients yields

$A_k\in\mathbb{C}^{n\times n}$9

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 $P(\lambda)=\sum_{k=0}^m A_k\lambda^k.$0 is an eigenvector of $P(\lambda)=\sum_{k=0}^m A_k\lambda^k.$1 for $P(\lambda)=\sum_{k=0}^m A_k\lambda^k.$2. Any nonzero column of $P(\lambda)=\sum_{k=0}^m A_k\lambda^k.$3 is a generalized eigenvector of rank at most $P(\lambda)=\sum_{k=0}^m A_k\lambda^k.$4. For fixed $P(\lambda)=\sum_{k=0}^m A_k\lambda^k.$5 and a fixed column index $P(\lambda)=\sum_{k=0}^m A_k\lambda^k.$6, the nonzero $P(\lambda)=\sum_{k=0}^m A_k\lambda^k.$7-th columns of

$P(\lambda)=\sum_{k=0}^m A_k\lambda^k.$8

form a Jordan chain ending in an eigenvector. If the last nonzero column occurs at index $P(\lambda)=\sum_{k=0}^m A_k\lambda^k.$9, then

$R(\lambda)$0

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,

$R(\lambda)$1

for $R(\lambda)$2 larger than the spectral abscissa of $R(\lambda)$3. Substituting the partial fraction expansion and using

$R(\lambda)$4

one obtains

$R(\lambda)$5

In terms of the nilpotent data,

$R(\lambda)$6

This is the standard Jordan-block expansion, recovered from the partial fractions of the resolvent.

Over the real field, if $R(\lambda)$7 factors into linear and quadratic factors, the decomposition may be written in real partial fractions. For a factor $R(\lambda)$8, terms of the form

$R(\lambda)$9

invert to

$R(\lambda)=P(\lambda)Q(\lambda)^{-1},$0

The method therefore accommodates complex conjugate spectral pairs without leaving real arithmetic.

For the linear system

$R(\lambda)=P(\lambda)Q(\lambda)^{-1},$1

the homogeneous solution is

$R(\lambda)=P(\lambda)Q(\lambda)^{-1},$2

Variation of parameters gives

$R(\lambda)=P(\lambda)Q(\lambda)^{-1},$3

If $R(\lambda)=P(\lambda)Q(\lambda)^{-1},$4 is a combination of exponentials and polynomials, or of sinusoids, the convolution can be evaluated explicitly because $R(\lambda)=P(\lambda)Q(\lambda)^{-1},$5 decomposes into sums of $R(\lambda)=P(\lambda)Q(\lambda)^{-1},$6 times polynomials in $R(\lambda)=P(\lambda)Q(\lambda)^{-1},$7. In Laplace form,

$R(\lambda)=P(\lambda)Q(\lambda)^{-1},$8

so the transfer function from input $R(\lambda)=P(\lambda)Q(\lambda)^{-1},$9 to state $P(\lambda)$0 is $P(\lambda)$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

$P(\lambda)$2

or, over $P(\lambda)$3, 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 $P(\lambda)$4 is computed symbolically, one may also write

$P(\lambda)$5

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 $P(\lambda)$6 are eigenvectors, while columns of $P(\lambda)$7 are generalized eigenvectors linked by

$P(\lambda)$8

The numerical limitations are explicit. Factoring $P(\lambda)$9 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 $A \in \mathbb{C}^{n\times n}$00 may be more stable numerically. For large $A \in \mathbb{C}^{n\times n}$01, computing $A \in \mathbb{C}^{n\times n}$02 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 $A \in \mathbb{C}^{n\times n}$03 case,

$A \in \mathbb{C}^{n\times n}$04

one has

$A \in \mathbb{C}^{n\times n}$05

and

$A \in \mathbb{C}^{n\times n}$06

From substitution at $A \in \mathbb{C}^{n\times n}$07 and $A \in \mathbb{C}^{n\times n}$08,

$A \in \mathbb{C}^{n\times n}$09

Therefore

$A \in \mathbb{C}^{n\times n}$10

Because the eigenvalues are distinct, only simple poles occur and there are no polynomial-in-$A \in \mathbb{C}^{n\times n}$11 factors.

In the defective $A \in \mathbb{C}^{n\times n}$12 example,

$A \in \mathbb{C}^{n\times n}$13

the characteristic polynomial is

$A \in \mathbb{C}^{n\times n}$14

so the only eigenvalue is $A \in \mathbb{C}^{n\times n}$15, with algebraic multiplicity $A \in \mathbb{C}^{n\times n}$16 and a single Jordan chain of length $A \in \mathbb{C}^{n\times n}$17. The decomposition

$A \in \mathbb{C}^{n\times n}$18

has

$A \in \mathbb{C}^{n\times n}$19

$A \in \mathbb{C}^{n\times n}$20

Hence

$A \in \mathbb{C}^{n\times n}$21

and the generalized eigenvector chain is read from

$A \in \mathbb{C}^{n\times n}$22

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 $A \in \mathbb{C}^{n\times n}$23 and approximates $A \in \mathbb{C}^{n\times n}$24 or $A \in \mathbb{C}^{n\times n}$25 by a partial fraction expansion such as

$A \in \mathbb{C}^{n\times n}$26

or, more generally,

$A \in \mathbb{C}^{n\times n}$27

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 (Bertaccini et al., 2017).

The conceptual relation is direct. In the sparse-matrix literature, the coefficients are usually scalar residues multiplying matrix resolvents. When $A \in \mathbb{C}^{n\times n}$28 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 \in \mathbb{C}^{n\times n}$29 with coefficients involving derivatives of $A \in \mathbb{C}^{n\times n}$30. 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 \in \mathbb{C}^{n\times n}$31, and an error decomposition into rational approximation error and solve or preconditioner error. It reports applications to $A \in \mathbb{C}^{n\times n}$32, $A \in \mathbb{C}^{n\times n}$33, fractional powers, localized matrices, PDE discretizations, Markov generator matrices, and large power systems matrices. The approach is advantageous when $A \in \mathbb{C}^{n\times n}$34 is sparse or localized, $A \in \mathbb{C}^{n\times n}$35 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 \in \mathbb{C}^{n\times n}$36.

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 $A \in \mathbb{C}^{n\times n}$37 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.