Polynomials and Matrices | KEVOS® Mathematics

Q: Why does the constant term become a multiple of the identity?

Because every term in the sum must be a matrix of the same size for the addition to be defined. The degree-zero term is $a_0A^{0}$, and $A^{0} = I_n$ by the same convention that makes $x^{0} = 1$ for scalars, so the term is $a_0I_n$. Writing it as a bare scalar is a widely used abbreviation, not a literal instruction.

Q: Does substituting into the factored form really give the same matrix?

Yes, and the reason is commutativity. Every factor $A - b_iI_n$ is a polynomial in $A$, and any two polynomials in the same matrix commute with one another. Expanding the product therefore behaves exactly as the scalar expansion does, and the factors may be applied in any order without changing the result.

Q: Can I substitute a matrix into a polynomial that involves division?

Not in general. Polynomials are built only from addition, scalar multiplication and powers, all of which square matrices support. Rational functions require inverses, which exist only for non-singular matrices, and even then $A^{-1}B$ and $BA^{-1}$ differ. Rational matrix functions are studied, but they are a substantially more delicate object.

Q: What does it mean for a polynomial to annihilate a matrix?

It means $p(A) = O$, the zero matrix. Such polynomials always exist — the Cayley-Hamilton theorem shows the characteristic polynomial is one — and every eigenvalue of $A$ must be a root of any annihilating polynomial, because $p(A)\vec{x} = p(\lambda)\vec{x}$ for an eigenvector $\vec{x}$. That is how factoring an annihilating polynomial exposes eigenvalues.

Q: Why is the factored form often preferred computationally?

Forming each factor $A - b_iI_n$ requires only subtractions on the diagonal, so no large intermediate coefficients ever appear, and the product needs $m-1$ multiplications, the same as Horner's scheme. The factored form also makes the structure visible: each factor annihilates one eigenspace, which is why a partially applied product has a recognisable rank and column space.

Q: How do I evaluate a very high degree polynomial in a matrix efficiently?

First reduce the degree. By Cayley-Hamilton, $A^{n}$ is a linear combination of lower powers, so dividing the polynomial by the characteristic polynomial and keeping the remainder leaves an equivalent polynomial of degree less than $n$. Then evaluate the remainder by Horner's scheme, or by the Paterson-Stockmeyer grouping if the degree is still large.

Q: What happens to eigenvalues when a polynomial is applied?

They are mapped through the same polynomial. If $A\vec{x} = \lambda\vec{x}$ then $p(A)\vec{x} = p(\lambda)\vec{x}$, so the eigenvectors are unchanged and each eigenvalue $\lambda$ becomes $p(\lambda)$. This is why $p(A) = O$ forces $p(\lambda) = 0$ for every eigenvalue, and why a factor $(A - bI_n)$ annihilates exactly the eigenspace for $b$.

Overview

A polynomial is built from three ingredients: powers of the variable, multiplication by scalar coefficients, and addition. Division never appears. Square matrices support exactly the same three operations — they can be multiplied by scalars, added to one another, and raised to positive integer powers by repeated matrix multiplication — and, crucially, all of those operations preserve the size of a square matrix. Substituting a matrix for the variable of a polynomial is therefore a well-defined operation with a matrix result.

The only genuine subtlety is the constant term. In the scalar polynomial $a_0 + a_1x + \cdots + a_mx^{m}$ the leading term is a number; in the matrix version it must be a matrix of the right size, and the only sensible candidate is $a_0I_n$ . This follows from the convention $A^{0} = I_n$ , which makes the constant term the degree-zero term like any other. Writing $p(A) = a_0 + a_1A + \cdots$ with a bare scalar is a common abbreviation, but it is an abbreviation, and taking it literally in code is one of the reliable ways to compute the wrong answer.

The second observation is that the factored form of a polynomial may be substituted just as freely as the expanded form, and gives the same matrix. That is not automatic: matrix multiplication does not commute in general, so expanding a product of matrix factors normally requires care about ordering. Here it is safe, because every factor is a polynomial in the single matrix $A$ , and any two polynomials in the same matrix commute with one another. The factored form is often both cheaper to evaluate and far more informative.

Matrix polynomials are not an idle generalisation. They are the machinery behind the proof that every square matrix has an eigenvalue, behind the Cayley-Hamilton theorem, behind Krylov subspace iterative solvers, behind pole-placement formulas in control, and behind closed-form rotation formulas in kinematics. Wherever a matrix appears inside an expression built from sums, scalar multiples and powers, a matrix polynomial is what is actually being evaluated.

Definition

Polynomial Evaluated at a Matrix

Let $p(x) = a_0 + a_1x + a_2x^{2} + \cdots + a_mx^{m}$ be a polynomial with coefficients in $\mathbb{C}$ , and let $A$ be a square matrix of size $n$ . Then

$p(A) = a_0I_n + a_1A + a_2A^{2} + \cdots + a_mA^{m}$ ,

a square matrix of size

n

. Each power

A^{k}

is the

k

-fold matrix product of

A

with itself, and the degree-zero term uses

A^{0} = I_n

The requirement that $A$ be square is essential twice over: a non-square matrix cannot be multiplied by itself, and the sum $a_0I_n + a_1A$ requires the identity to have the same shape as $A$ .

Matrix Power

For a square matrix $A$ and a non-negative integer $k$ , the power $A^{k}$ is defined recursively by $A^{0} = I_n$ and $A^{k} = A\,A^{k-1}$ for $k \ge 1$ . The convention $A^{0} = I_n$ is what makes the index laws $A^{j}A^{k} = A^{j+k}$ hold without exception, and it is the reason the constant term of a polynomial becomes a multiple of the identity.

This is not the entrywise power. Raising each entry of $A$ to the power $k$ gives a completely different matrix, and several programming languages make the entrywise operation the default for the ** operator.

Annihilating Polynomial

A non-zero polynomial $p$ with $p(A) = O$ , the zero matrix, is said to annihilate $A$ . A weaker and equally useful notion fixes a single vector: $p$ annihilates $A$ at $\vec{x}$ when $p(A)\vec{x} = \vec{0}$ . Such polynomials always exist, and factoring them is the route by which eigenvalues can be extracted without ever computing a determinant.

Concepts

The constant term is a matrix

Every term of $p(A)$ must be an $n \times n$ matrix, because they are being added together. The terms $a_1A$ , $a_2A^{2}$ and so on are matrices automatically. The degree-zero term is $a_0A^{0}$ , and since $A^{0} = I_n$ this is the matrix $a_0I_n$ : the scalar $a_0$ on the diagonal and zeros everywhere else. It is not the scalar $a_0$ , and it is emphatically not the matrix with $a_0$ in every position. This distinction is invisible in the standard abbreviated notation and highly visible in a wrong answer.

Why the factored form may be substituted

Matrix multiplication is not commutative, so expanding $(A - bI_n)(A - cI_n)$ appears to demand care. It does not, because $A$ commutes with itself and with every multiple of the identity. Consequently any two polynomials in the same matrix $A$ commute: $p(A)q(A) = q(A)p(A)$ , and both equal $(pq)(A)$ , the polynomial product evaluated at $A$ . Substituting $A$ into a factored polynomial therefore gives exactly the same matrix as substituting into its expanded form, and the factors may be applied in any order.

The factored form is usually the better one

Beyond convenience, the factored form carries structure. Each factor $A - b_iI_n$ is a shifted matrix, and its null space is the eigenspace of $A$ for the eigenvalue $b_i$ . A product of such factors annihilates any vector lying in the eigenspace of any of the roots. That observation drives the constructive proof that every square matrix has an eigenvalue: build an annihilating polynomial for a chosen vector, factor it, and apply the factors one at a time until the result becomes the zero vector. The last factor applied identifies an eigenvalue.

Action on an eigenvector

If $A\vec{x} = \lambda\vec{x}$ then $A^{k}\vec{x} = \lambda^{k}\vec{x}$ by repeated application, and therefore $p(A)\vec{x} = p(\lambda)\vec{x}$ for every polynomial $p$ . An eigenvector of $A$ is an eigenvector of every polynomial in $A$ , with the eigenvalue transformed by the same polynomial. This is the spectral mapping property in its most elementary form, and it explains why $p(A) = O$ forces every eigenvalue of $A$ to be a root of $p$ .

Degree can always be reduced

The powers $I_n, A, A^{2}, \ldots$ cannot be independent forever: the space of $n \times n$ matrices has dimension $n^{2}$ , so at most $n^{2}$ of them can be linearly independent, and the Cayley-Hamilton theorem sharpens this dramatically by showing that $A^{n}$ is a combination of the lower powers. Every polynomial in $A$ therefore equals a polynomial of degree less than $n$ . In practice this means that evaluating a high-degree matrix polynomial should begin by reducing the polynomial modulo the characteristic polynomial, not by computing high powers.

Powers are where the numerical trouble lives

The entries of $A^{k}$ grow or decay roughly like $\rho(A)^{k}$ , where $\rho$ is the spectral radius. For $\rho > 1$ the powers overflow; for $\rho < 1$ they underflow to zero; and for a non-normal matrix the intermediate powers can be enormous even when both $A$ and $A^{k}$ are modest, producing catastrophic cancellation when the terms are summed. Matrix polynomials are also destructive of sparsity: each multiplication by $A$ propagates non-zeros further from the diagonal, so $A^{k}$ of a sparse matrix is generally dense.

Evaluating a matrix polynomial

Confirm the matrix is squarePowers and the addition of

a_0I_n

both require it. A rectangular matrix cannot be substituted into a polynomial at all.

Reduce the degree if it is largeAny polynomial in

A

equals one of degree less than

n

. Divide by the characteristic polynomial and keep the remainder before computing anything.

Factor if the roots are knownA product of shifted matrices

A - b_iI_n

is usually cheaper to form and reveals the eigenstructure. The factors commute, so the order is free.

Otherwise use Horner's schemeEvaluate as

(\cdots((a_mA + a_{m-1}I_n)A + a_{m-2}I_n)A + \cdots)

, which needs

m-1

matrix multiplications and never forms a high power explicitly.

Write the constant term as $a_0I_n$ Never add a bare scalar to a matrix. In array-based languages that operation silently adds the scalar to every entry.

VerifyCheck against the factored form, or apply the result to a known eigenvector and confirm the output equals

p(\lambda)\vec{x}

Equations

Polynomial evaluated at a matrix

EQ-PM-01

p (A) = a_{0} I_{n} + a_{1} A + a_{2} A^{2} + \dots + a_{m} A^{m} = \sum_{k}^{=} a_{k} A^{k}

Each term is an $n \times n$ matrix. The degree-zero term is $a_0 I_n$ , not the scalar $a_0$ .

Matrix powers

EQ-PM-02

A^{0} = I_{n}, A^{k} = A A^{k - 1}, A^{j} A^{k} = A^{j + k}

The recursive definition of a matrix power and the index law it satisfies. The convention $A^{0} = I_n$ is what makes the index law hold for $j = 0$ .

Factored form substituted

EQ-PM-03

p (x) = \prod_{i}^{=} (x - b_{i}) \Rightarrow p (A) = (A - b_{1} I_{n}) (A - b_{2} I_{n}) \dots (A - b_{m} I_{n})

Over $\mathbb{C}$ every polynomial factors into linear factors. The matrix factors commute, so the ordering is immaterial and the product equals the expanded form.

Polynomials in one matrix commute

EQ-PM-04

p (A) q (A) = q (A) p (A) = (p q) (A)

The reason the factored substitution is legitimate. It fails for polynomials in two different matrices unless those matrices themselves commute.

Action on an eigenvector

EQ-PM-05

A \vec{x} = λ \vec{x} \Rightarrow p (A) \vec{x} = p (λ) \vec{x}

An eigenvector of $A$ is an eigenvector of every polynomial in $A$ , with the eigenvalue mapped through the same polynomial.

Horner's scheme for a matrix argument

EQ-PM-06

p (A) = (\dots ((a_{m} A + a_{m - 1} I_{n}) A + a_{m - 2} I_{n}) A + \dots) + a_{0} I_{n}

Evaluates a degree- $m$ polynomial in $m-1$ matrix multiplications without ever forming a high power explicitly, which limits both cost and overflow.

Annihilating polynomial at a vector

EQ-PM-07

p (A) \vec{x} = \vec{0}, \vec{x} \neq \vec{0}

The relation that a linear dependence among $\vec{x}, A\vec{x}, A^{2}\vec{x}, \ldots$ produces. Factoring such a $p$ and applying the factors in turn extracts an eigenvalue of $A$ .

Variable Definitions

Symbols used on this page
Symbol	Name	Meaning	Domain / type
$A$	Matrix argument	The square matrix substituted for the polynomial variable	n x n matrix over C
$n$	Size	Common row and column count of $A$ and of $I_n$	positive integer
$p (x)$	Polynomial	The scalar polynomial whose variable is replaced by a matrix	polynomial with complex coefficients
$m$	Degree	Degree of the polynomial, the highest power appearing	non-negative integer
$a_{k}$	Coefficient	Scalar multiplier of the $k$ -th power term	complex scalar
$b_{i}$	Root	A root of $p$ , appearing in the factor $A - b_iI_n$	complex scalar
$I_{n}$	Identity matrix	The size- $n$ identity, equal to $A^{0}$ and carrier of the constant term	n x n matrix
$A^{k}$	Matrix power	The $k$ -fold matrix product of $A$ with itself; not the entrywise power	n x n matrix
$O$	Zero matrix	The matrix with every entry zero, the result when $p$ annihilates $A$	n x n matrix

Worked Numerical Example

Problem statement

A discrete-time state transition matrix in companion form governs a third-order digital filter. Evaluate the quadratic $q(x) = x^{2} - 3x + 2$ at that matrix, first from the expanded form and then from the factored form $(x-1)(x-2)$ , and interpret the structure of the result.

State the matrix and the polynomial
The companion matrix $D$ shifts the state and forms the new entry from the feedback coefficients in its last row. The polynomial to be evaluated is $q(x) = x^{2} - 3x + 2$ , whose constant term will contribute $2I_3$ .
$D = [\begin{matrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 6 & - 11 & 6 \end{matrix}], q (x) = x^{2} - 3 x + 2$
Compute the required powers
Only $D^{0}$ , $D^{1}$ and $D^{2}$ are needed. Because the first two rows of $D$ are shifted unit vectors, the first two rows of $D^{2}$ are simply the second and third rows of $D$ ; only the last row requires arithmetic, as $6(0,1,0) - 11(0,0,1) + 6(6,-11,6) = (36, -60, 25)$ .
$D^{0} = I_{3} = [\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}], D^{2} = [\begin{matrix} 0 & 0 & 1 \\ 6 & - 11 & 6 \\ 36 & - 60 & 25 \end{matrix}]$
Assemble the expanded form
Combine the three matrices with the coefficients $1$ , $-3$ and $2$ . The constant term enters as $2I_3$ , contributing $2$ to each diagonal entry and nothing elsewhere — not $2$ to every entry.
$q (D) = D^{2} - 3 D + 2 I_{3} = [\begin{matrix} 0 & 0 & 1 \\ 6 & - 11 & 6 \\ 36 & - 60 & 25 \end{matrix}] - [\begin{matrix} 0 & 3 & 0 \\ 0 & 0 & 3 \\ 18 & - 33 & 18 \end{matrix}] + [\begin{matrix} 2 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 2 \end{matrix}]$
Complete the entrywise arithmetic
Adding the three matrices entry by entry gives the result. For instance the $(3,3)$ entry is $25 - 18 + 2 = 9$ , and the $(2,2)$ entry is $-11 - 0 + 2 = -9$ .
$q (D) = [\begin{matrix} 2 & - 3 & 1 \\ 6 & - 9 & 3 \\ 18 & - 27 & 9 \end{matrix}]$
Form the factored version
The polynomial factors as $q(x) = (x-1)(x-2)$ , so $q(D) = (D - I_3)(D - 2I_3)$ . Each factor is $D$ with a constant subtracted from the diagonal only.
$D - I_{3} = [\begin{matrix} - 1 & 1 & 0 \\ 0 & - 1 & 1 \\ 6 & - 11 & 5 \end{matrix}], D - 2 I_{3} = [\begin{matrix} - 2 & 1 & 0 \\ 0 & - 2 & 1 \\ 6 & - 11 & 4 \end{matrix}]$
Multiply the factors and compare
The first row of the product is $(-1)(-2,1,0) + (1)(0,-2,1) + (0)(6,-11,4) = (2,-3,1)$ , matching the expanded result. The remaining rows agree likewise. Because polynomials in $D$ commute, reversing the two factors gives the same matrix.
$(D - I_{3}) (D - 2 I_{3}) = [\begin{matrix} 2 & - 3 & 1 \\ 6 & - 9 & 3 \\ 18 & - 27 & 9 \end{matrix}] = q (D)$
Read the structure of the result
Every row of $q(D)$ is a multiple of $(2, -3, 1)$ and every column is a multiple of $(1, 3, 9)$ , so $q(D)$ has rank one. That column direction is not accidental: $D(1,3,9) = (3, 9, 6 - 33 + 54) = (3, 9, 27) = 3(1,3,9)$ , so it is an eigenvector of $D$ for the eigenvalue $3$ — the one root of the cubic $x^{3} - 6x^{2} + 11x - 6$ that $q$ does not remove.
$D [\begin{matrix} 1 \\ 3 \\ 9 \end{matrix}] = [\begin{matrix} 3 \\ 9 \\ 27 \end{matrix}] = 3 [\begin{matrix} 1 \\ 3 \\ 9 \end{matrix}]$
Extend to the annihilating polynomial
Multiplying by the remaining factor gives the zero matrix, since $q(D)(D - 3I_3)$ has first row $2(-3,1,0) - 3(0,-3,1) + 1(6,-11,3) = (0,0,0)$ and the other rows are multiples of the first. Hence $p(D) = O$ for $p(x) = (x-1)(x-2)(x-3) = x^{3} - 6x^{2} + 11x - 6$ , which is exactly the characteristic polynomial encoded in the last row of the companion matrix.
$(D - I_{3}) (D - 2 I_{3}) (D - 3 I_{3}) = O$

Result

The quadratic evaluates to a rank-one matrix whose column space is the eigenspace of $D$ for the eigenvalue $3$ : the factors $(D - I_3)$ and $(D - 2I_3)$ have annihilated the other two eigendirections and left only the third. Applying the final factor annihilates that one too, so the cubic is an annihilating polynomial for $D$ . For the filter this is the statement that its state transition matrix satisfies its own difference equation — the feedback coefficients in the last row of $D$ are precisely the coefficients of the polynomial that kills it.

Applications & Industry Use

Control systems

Pole placement and state transition

Ackermann's formula computes a state-feedback gain as a row vector times the desired characteristic polynomial evaluated at the plant matrix, $p_{d}(A)$ . Separately, the discrete-time state transition over $k$ steps is $A^{k}$ , and Cayley-Hamilton lets that power be written as a polynomial of degree less than $n$ , which is how transition matrices are evaluated without repeated multiplication.

Robotics and kinematics

Closed-form rotation from a skew matrix

Rodrigues' formula expresses a rotation about a unit axis as $R = I_3 + \sin\theta\,K + (1 - \cos\theta)K^{2}$ , where $K$ is the skew-symmetric cross-product matrix of the axis. This is a quadratic polynomial in $K$ ; higher powers are unnecessary because $K^{3} = -K$ collapses them back. Evaluating it is a matrix polynomial evaluation performed millions of times per second in motion planning code.

Numerical linear algebra

Krylov subspace iterative solvers

Conjugate gradients, GMRES and their relatives construct approximate solutions of the form $\vec{x}_k = \vec{x}_0 + q(A)\vec{r}_0$ for a polynomial $q$ of degree $k-1$ . The residual after $k$ steps is $p(A)\vec{r}_0$ with $p(0) = 1$ , and the entire convergence theory of these methods is the question of how small such a polynomial can be made on the spectrum of $A$ .

Network analysis

Walk counting and graph polynomials

For an adjacency matrix $A$ , the entry $[A^{k}]_{ij}$ counts walks of length exactly $k$ from node $i$ to node $j$ . A polynomial in $A$ therefore counts weighted combinations of walks, which is how communicability, Katz centrality and network resilience measures are defined and computed in transport and communication network engineering.

Signal processing

Filters as polynomials in a shift operator

A finite impulse response filter applied to a finite signal is a polynomial in the shift matrix, with the filter taps as coefficients. Cascading two filters corresponds to multiplying their polynomials, and because polynomials in one matrix commute, the order in which two filters are applied does not change the result — a fact filter designers rely on when reordering a processing chain.

Structural and thermal analysis

Polynomial approximation of matrix functions

Transient heat conduction and structural response require the matrix exponential $e^{At}$ , which is approximated in practice by a truncated or Chebyshev polynomial in $A$ , or by a rational approximation whose numerator and denominator are both matrix polynomials. The evaluation strategy — Horner, Paterson-Stockmeyer, scaling and squaring — is chosen entirely on matrix-polynomial grounds.

Design Considerations

Write the constant term as $a_0I_n$ in code as well as on paper

In array-based languages, adding a scalar to a matrix broadcasts: A + 2 adds $2$ to every entry, not to the diagonal. The correct expression is A + 2*np.eye(n). The bug produces a plausible matrix of the right shape and is easily missed in review, so it should be caught by an explicit unit test on a diagonal input.

Reduce the degree before evaluating

By Cayley-Hamilton, any polynomial in $A$ equals one of degree less than $n$ . For a high-degree polynomial, divide by the characteristic polynomial and evaluate only the remainder. This can turn an evaluation requiring hundreds of matrix multiplications into one requiring a handful, and it also avoids forming the high powers where overflow lives.

Never form high powers explicitly

Computing $A^{20}$ by nineteen multiplications is both slow and numerically poor. Horner's scheme evaluates a degree- $m$ polynomial in $m-1$ multiplications with no explicit high powers; the Paterson-Stockmeyer scheme reduces that to about $2\sqrt{m}$ . Where only $A^{k}$ itself is wanted, binary exponentiation needs about $\log_2 k$ multiplications.

Expect sparsity to be destroyed

Each multiplication by a sparse $A$ spreads non-zeros further, so $A^{k}$ is generally dense even when $A$ is not. Matrix polynomials of sparse matrices should therefore be applied to vectors rather than formed as matrices: computing $p(A)\vec{v}$ by repeated sparse matrix-vector products costs $m$ sparse products and never materialises a dense intermediate.

Prefer the factored form when the roots are known and real

Forming $\prod(A - b_iI_n)$ requires only diagonal subtractions plus $m-1$ multiplications, avoids computing coefficients that may be enormous, and exposes the eigenstructure. If the roots are complex and the matrix real, pair conjugate roots into real quadratic factors to keep the arithmetic real.

Watch conditioning for non-normal matrices

For a matrix that is far from normal, intermediate powers can be many orders of magnitude larger than either the matrix or the final answer, so summing the terms of an expanded polynomial cancels catastrophically. Where accuracy matters, evaluate through a Schur decomposition or a scaling-and-squaring scheme rather than by direct summation.

Standards & Codes

Notation, interchange and numerical standards that govern how this material is written down, stored and computed in production systems.

Applicable standards, conventions and reference implementations
Reference	Title	Relevance to this topic
`ISO 80000-2`	Quantities and units — Part 2: Mathematics	Governs the notation for matrix powers, the identity matrix $I_n$ and polynomial coefficients, and the typographic distinction between the scalar variable $x$ and the matrix argument $A$ .
`BLAS Level 3 (xGEMM)`	Basic Linear Algebra Subprograms — matrix-matrix operations	Every matrix multiplication in a polynomial evaluation is a `GEMM` call. Its accumulate form $C \leftarrow \alpha AB + \beta C$ implements one Horner step directly, which is why Horner's scheme maps so cleanly onto optimised libraries.
`IEEE 754-2019`	IEEE Standard for Floating-Point Arithmetic	Defines the overflow and underflow limits that matrix powers reach quickly when the spectral radius departs from one, and the rounding model behind the cancellation seen in expanded-form evaluation.
`ISO/IEC 14882`	Programming languages — C++	Relevant to the operator-overloading conventions that determine whether `A * A` means a matrix product or an entrywise product; the Eigen library distinguishes them explicitly with `A * A` and `A.array() * A.array()`.
`IEC 61131-3`	Programmable controllers — Programming languages	Applies where a small matrix polynomial — a Rodrigues rotation or a fixed-gain state observer — is evaluated inside deterministic control logic, constraining the arithmetic types and the permissible loop structure.

Material Selection

For a mathematical topic, "material" is the numeric representation: the scalar field, storage format and precision the computation is built from.

Representation and precision selection
Representation	Select when	Trade-off
Exact integer arithmetic	Companion matrices, adjacency matrices and combinatorial walk counting where entries are integers and the answer must be exact.	Exactly correct with no rounding, but the entries of $A^{k}$ grow geometrically and a fixed-width type will overflow for modest $k$ .
Exact rational or symbolic coefficients	Deriving an annihilating or characteristic polynomial, or evaluating a polynomial containing a design parameter.	Yields a factorable exact result, but expression size grows rapidly with the degree and the matrix size.
IEEE 754 binary64, dense storage	General numerical evaluation of moderate-degree polynomials at matrices of moderate size.	Fast, well supported by BLAS Level 3, but vulnerable to overflow at high powers and to cancellation for non-normal matrices.
Sparse storage with matrix-vector application	Very large sparse matrices where only $p(A)\vec{v}$ is required, as in Krylov solvers and network centrality measures.	Keeps memory linear in the number of non-zeros and never densifies, but cannot deliver the matrix $p(A)$ itself.
IEEE 754 binary32 on accelerators	Graphics and robotics kernels such as Rodrigues rotations, evaluated at very high rates and at low degree.	Halves bandwidth and maps well to hardware, but the reduced exponent range makes even modest powers risky and accumulated drift requires periodic re-normalisation.
Finite field arithmetic	Coding theory and cryptography, where companion matrices over $\mathbb{F}_q$ implement polynomial arithmetic.	No overflow and exact results with fixed operand size, but the factorisation of a polynomial over a finite field differs from its factorisation over $\mathbb{C}$ .

Manufacturing Notes

Implementation notes — how the result is actually produced by hand, by algorithm and by library, including cost and numerical behaviour.

Operation counts

Naive evaluation of a degree- $m$ polynomial by forming each power separately costs $m-1$ matrix multiplications, each about $2n^{3}$ operations. Horner's scheme costs the same number of multiplications but avoids storing the powers. The Paterson-Stockmeyer scheme reduces the count to roughly $2\sqrt{m}$ multiplications by grouping terms, which is the standard technique inside matrix-exponential implementations.

Hand procedure

Compute the powers one at a time, each from the previous, and write them out fully before combining. Scale each by its coefficient, then add entry by entry, taking the diagonal contribution of $a_0I_n$ last so that it is not forgotten. Where the polynomial factors, form the shifted matrices instead — subtracting a constant from the diagonal is far less error-prone than computing a cube.

Library behaviour

numpy.linalg.matrix_power performs true matrix exponentiation by binary powering, whereas A ** 2 on a NumPy array is entrywise. numpy.polyvalfromroots and numpy.polyval are scalar routines and must not be applied to a matrix argument. SymPy evaluates matrix polynomials symbolically and Matrix.charpoly() supplies the characteristic polynomial for degree reduction. scipy.linalg.funm and expm implement matrix functions using polynomial and rational approximations internally.

Verification

Three checks are worth running. Compare the expanded and factored evaluations, which exercise entirely different arithmetic. Apply the result to a known eigenvector and confirm the output equals $p(\lambda)\vec{x}$ . Finally, check the trace: the trace of $p(A)$ equals $\sum_i p(\lambda_i)$ over the eigenvalues with multiplicity, which is a cheap scalar test of the whole computation.

Numerical stability

Direct summation of an expanded matrix polynomial is stable only when the terms do not vary wildly in magnitude. For non-normal matrices, or for polynomials of high degree, use a Schur decomposition to reduce $A$ to triangular form and evaluate there, or apply scaling and squaring so that every polynomial evaluation is performed on a matrix of small norm.

Failure Modes & Common Mistakes

Failure modes, root causes and prevention
Failure mode / mistake	Impact	Root cause	Prevention & detection
Adding the constant term as a scalar	high	Writing `A + a0` in an array language, which broadcasts the scalar to every entry instead of forming $a_0I_n$ .	Always write $a_0I_n$ explicitly. Unit-test with a diagonal matrix, where the wrong result is immediately visible in the off-diagonal entries.
Using entrywise powers instead of matrix powers	high	The `**` or `.^` operator meaning elementwise exponentiation in NumPy, MATLAB and similar environments.	Use the explicit matrix-power function (`numpy.linalg.matrix_power`, MATLAB's `^`) and verify on a matrix where the two differ, such as one with a zero on the diagonal.
Taking $A^{0}$ to be the zero matrix	medium	Carrying over the intuition that a zero exponent contributes nothing, which deletes the entire constant term.	Fix $A^{0} = I_n$ as a stated convention and check that a constant polynomial $p(x) = c$ evaluates to $cI_n$ .
Substituting into a non-square matrix	medium	Applying a polynomial to a rectangular matrix, where neither $A^{2}$ nor $A + I$ is defined.	Assert squareness at the entry point of any routine that evaluates a matrix polynomial.
Overflow when forming high powers	high	Computing $A^{k}$ directly for large $k$ when the spectral radius exceeds one, so intermediate entries leave the representable range.	Reduce the polynomial degree using the characteristic polynomial, use Horner's scheme, or scale the matrix before evaluation and undo the scaling afterwards.
Catastrophic cancellation in the expanded form	medium	Summing terms of vastly different magnitude for a non-normal matrix, so the significant digits of the answer are lost.	Evaluate through the factored form or a Schur decomposition, and compare the two routes on a representative case.
Densifying a sparse matrix	medium	Forming $p(A)$ explicitly for a large sparse $A$ , where each power fills in more entries until memory is exhausted.	Apply the polynomial to a vector by repeated sparse matrix-vector products rather than forming the polynomial as a matrix.
Assuming factors of polynomials in two different matrices commute	low	Over-generalising the commutation property, which holds for polynomials in a single matrix but not for $p(A)$ and $q(B)$ with $AB \ne BA$ .	Check that every factor is a polynomial in the same matrix before reordering. If two matrices are involved, verify $AB = BA$ explicitly.

FAQs

Why does the constant term become a multiple of the identity?

Because every term in the sum must be a matrix of the same size for the addition to be defined. The degree-zero term is $a_0A^{0}$ , and $A^{0} = I_n$ by the same convention that makes $x^{0} = 1$ for scalars, so the term is $a_0I_n$ . Writing it as a bare scalar is a widely used abbreviation, not a literal instruction.

Does substituting into the factored form really give the same matrix?

Yes, and the reason is commutativity. Every factor $A - b_iI_n$ is a polynomial in $A$ , and any two polynomials in the same matrix commute with one another. Expanding the product therefore behaves exactly as the scalar expansion does, and the factors may be applied in any order without changing the result.

Can I substitute a matrix into a polynomial that involves division?

Not in general. Polynomials are built only from addition, scalar multiplication and powers, all of which square matrices support. Rational functions require inverses, which exist only for non-singular matrices, and even then $A^{-1}B$ and $BA^{-1}$ differ. Rational matrix functions are studied, but they are a substantially more delicate object.

What does it mean for a polynomial to annihilate a matrix?

It means $p(A) = O$ , the zero matrix. Such polynomials always exist — the Cayley-Hamilton theorem shows the characteristic polynomial is one — and every eigenvalue of $A$ must be a root of any annihilating polynomial, because $p(A)\vec{x} = p(\lambda)\vec{x}$ for an eigenvector $\vec{x}$ . That is how factoring an annihilating polynomial exposes eigenvalues.

Why is the factored form often preferred computationally?

Forming each factor $A - b_iI_n$ requires only subtractions on the diagonal, so no large intermediate coefficients ever appear, and the product needs $m-1$ multiplications, the same as Horner's scheme. The factored form also makes the structure visible: each factor annihilates one eigenspace, which is why a partially applied product has a recognisable rank and column space.

How do I evaluate a very high degree polynomial in a matrix efficiently?

First reduce the degree. By Cayley-Hamilton, $A^{n}$ is a linear combination of lower powers, so dividing the polynomial by the characteristic polynomial and keeping the remainder leaves an equivalent polynomial of degree less than $n$ . Then evaluate the remainder by Horner's scheme, or by the Paterson-Stockmeyer grouping if the degree is still large.

What happens to eigenvalues when a polynomial is applied?

They are mapped through the same polynomial. If $A\vec{x} = \lambda\vec{x}$ then $p(A)\vec{x} = p(\lambda)\vec{x}$ , so the eigenvectors are unchanged and each eigenvalue $\lambda$ becomes $p(\lambda)$ . This is why $p(A) = O$ forces $p(\lambda) = 0$ for every eigenvalue, and why a factor $(A - bI_n)$ annihilates exactly the eigenspace for $b$ .

References

Beezer, R. A. A First Course in Linear Algebra, Version 0.70. University of Puget Sound, 2006. Section EE, Subsection PM. Licensed under the GNU Free Documentation License v1.2.
ISO 80000-2:2019, Quantities and units — Part 2: Mathematics. International Organization for Standardization.
Higham, N. J. Functions of Matrices: Theory and Computation. Society for Industrial and Applied Mathematics, 2008.
Golub, G. H. and Van Loan, C. F. Matrix Computations, 4th edition. Johns Hopkins University Press, 2013.
Paterson, M. S. and Stockmeyer, L. J. On the Number of Nonscalar Multiplications Necessary to Evaluate Polynomials. SIAM Journal on Computing, 1973.

AI Suggested Questions

Show how the Paterson-Stockmeyer scheme evaluates a degree-16 matrix polynomial in about eight matrix multiplications instead of fifteen.
Given a $4 \times 4$ matrix, reduce $A^{20}$ to a polynomial of degree three using the characteristic polynomial, and verify the reduction numerically.
Why does $K^{3} = -K$ hold for the skew-symmetric cross-product matrix of a unit vector, and how does that truncate Rodrigues' rotation formula?
Construct a non-normal matrix for which the expanded and factored evaluations of the same polynomial differ noticeably in floating point.
Explain how the residual of GMRES after $k$ steps is a matrix polynomial applied to the initial residual, and what that implies about convergence.
For a sparse adjacency matrix, compare the cost of forming $p(A)$ explicitly with applying $p(A)$ to a single vector.

Related Calculators

Matrix Polynomial Evaluator

Evaluate any polynomial at a square matrix from coefficients or roots, showing every power and partial sum.

Matrix Power Calculator

Raise a square matrix to an integer power by binary exponentiation, with growth and conditioning diagnostics.

Characteristic Polynomial Calculator

Compute the characteristic polynomial of a matrix and use it to reduce the degree of any polynomial in that matrix.

Overview

Definition

Polynomial Evaluated at a Matrix

Matrix Power

Annihilating Polynomial

Concepts

The constant term is a matrix

Why the factored form may be substituted

The factored form is usually the better one

Action on an eigenvector

Degree can always be reduced

Powers are where the numerical trouble lives

Evaluating a matrix polynomial

Equations

Variable Definitions

Worked Numerical Example

Problem statement

State the matrix and the polynomial

Compute the required powers

Assemble the expanded form

Complete the entrywise arithmetic

Form the factored version

Multiply the factors and compare

Read the structure of the result

Extend to the annihilating polynomial

Applications &amp; Industry Use

Pole placement and state transition

Closed-form rotation from a skew matrix

Krylov subspace iterative solvers

Walk counting and graph polynomials

Filters as polynomials in a shift operator

Polynomial approximation of matrix functions

Design Considerations

Write the constant term as a0In in code as well as on paper

Reduce the degree before evaluating

Never form high powers explicitly

Expect sparsity to be destroyed

Prefer the factored form when the roots are known and real

Watch conditioning for non-normal matrices

Standards &amp; Codes

Material Selection

Manufacturing Notes

Operation counts

Hand procedure

Library behaviour

Verification

Numerical stability

Failure Modes &amp; Common Mistakes

FAQs

Related Engineering Topics

References

AI Suggested Questions

Related Calculators

Continue learning

Applications & Industry Use

Write the constant term as $a_0I_n$ in code as well as on paper

Standards & Codes

Failure Modes & Common Mistakes