Generalized Inverses of Polynomial Matrices

Abstract

A polynomial matrix is a matrix whose entries are polynomials. Equivalently, a polynomial matrix can be expressed as a polynomial with matrix coefficients. Formally speaking, in the univariable case, $(\mathbb {R}[x])^{m \times n}$ and $(\mathbb {R}^{m \times n}) [x]$ are isomorphic. In other words, extending the entries of matrices to polynomials is the same as extending the coefficients of polynomials to matrices.

A polynomial matrix is a matrix whose entries are polynomials. Equivalently, a polynomial matrix can be expressed as a polynomial with matrix coefficients. Formally speaking, in the univariable case, $(\mathbb {R}[x])^{m \times n}$ and $(\mathbb {R}^{m \times n}) [x]$ are isomorphic. In other words, extending the entries of matrices to polynomials is the same as extending the coefficients of polynomials to matrices. An example of a $3 \times 2$ polynomial matrix of degree 2:

$$ \left[ \begin{array}{cc} 1 &{} x^2 \\ x &{} 0 \\ x+1 &{} x^2 - 1 \end{array} \right] = \left[ \begin{array}{cc} 0 &{} \;1 \\ 0 &{} \;0 \\ 0 &{} \;1 \end{array} \right] x^2 + \left[ \begin{array}{cc} 0 &{} \;0 \\ 1 &{} \;0 \\ 1 &{} \;0 \end{array} \right] x + \left[ \begin{array}{cr} 1 &{} 0 \\ 0 &{} 0 \\ 1 &{} -1 \end{array} \right] . $$

In this chapter, we study the Moore-Penrose and Drazin inverses of a polynomial matrix and algorithms for computing the generalized inverses.

10.1 Introduction

We start with the scalar nonsingular case. Let $A \in \mathbb {R}^{n \times n}$ be nonsingular and

$$\begin{aligned} p( \lambda ) = \det (\lambda I - A) = c_0 \lambda ^n + c_1 \lambda ^{n-1} + \cdots + c_{n-1} \lambda + c_n , \end{aligned}$$

(10.1.1)

where $c_0 = 1$, be the characteristic polynomial of A. The Cayley-Hamilton theorem says

$$ c_0 A^n + c_1 A^{n-1} + \cdots + c_{n-1} A + c_n I = \mathrm{O}. $$

Thus, we have

$$ A^{-1} = -c_n^{-1} P_n = (-1)^{n-1} \frac{P_n}{\det (A)} , $$

where

$$ P_n = c_0 A^{n-1} + c_1 A^{n-2} + \cdots c_{n-2}A + c_{n-1}I, $$

which can be efficiently computed using the iterative Horner’s Rule:

$$ P_0 = \mathrm{O}; \ c_0 = 1; \quad P_i = AP_{i-1} + c_{i-1} I, \ i = 1,...,n. $$

How can the coefficients $c_i$, $i=1,...,n$ be obtained efficiently? It turns out that

$$ c_i = - \frac{\mathrm{tr}(AP_i)}{i}, \quad i=1,...n. $$

Putting all things together, the following algorithm, known as the Faddeev-LeVerrier algorithm [1, 2] and presented in Sect. 5.5, efficiently computes the coefficients $c_i$, $i=1,...,n$, of the characteristic polynomial of A.

Algorithm 10.1.1

Given a nonsingular matrix $A \in \mathbb {C}^{n \times n}$, this algorithm computes the coefficients $c_i$, $i=1,...,n$, of its characteristic polynomial (10.1.1) and its inverse $A^{-1}$.

• 1. $P_0 = \mathrm{O}$; $c_0 = 1$;

• 2. for $i=1$ to n

  • $P_i = AP_{i-1} + c_{i-1} I$;

  • $c_i = -\mathrm{tr}(AP_i) / i$;

• 3. $A^{-1} = -c_n^{-1} P_n = (-1)^{n-1} P_n / \det (A)$;

As pointed out in Sect. 5.5, Decell generalized the above algorithm to general scalar matrices and the Moore-Penrose inverse [3]. Let $B = AA^*$ and

$$\begin{aligned} f( \lambda ) = \det (\lambda I - B) = a_0 \lambda ^n + a_1 \lambda ^{n-1} + \cdots + a_{n-1} \lambda + a_n, \quad a_0 = 1, \end{aligned}$$

(10.1.2)

be the characteristic polynomial of B. If $k>0$ is the largest integer such that $a_k \ne 0$, then the Moore-Penrose inverse of A is given by

$$ A^{\dagger } = -a_k^{-1} A^* (B^{k-1} + a_1 B^{k-2} + \cdots + a_{k-1} I). $$

If $k=0$ is the largest integer such that $a_k \ne 0$, that is, $a_0$ is the only nonzero coefficient, then $A^{\dagger } = \mathrm{O}$.

Analogous to Algorithm 10.1.1, the following Decell’s algorithm computes $A^{\dagger }$.

Algorithm 10.1.2

[3] Given $A \in \mathbb {C}^{m \times n}$ and k the largest integer such that $a_k$ in (10.1.2) is nonzero, this algorithm computes the Moore-Penrose inverse of A.

• 1. if $k=0$ return $A^{\dagger } = \mathrm{O}$;

• 2. $B = AA^*$;

• 3. $P_0 = \mathrm{O}$; $A_0 = \mathrm{O}$; $a_0 = 1$;

• 4. for $i=1$ to k

  • $P_i = A_{i-1} + a_{i-1} I$;

  • $A_i = BP_i$;

  • $a_i = -\mathrm{tr}(A_i) / i$;

• 5. $A^{\dagger } = -a_k^{-1} A^* P_k$;

In the following sections, we generalize the above algorithm to polynomial matrices and their generalized inverses.

10.2 Moore-Penrose Inverse of a Polynomial Matrix

The definition of the Moore-Penrose inverse of a polynomial matrix is the same as the scalar case, that is, it is the polynomial matrix satisfying the four Penrose conditions.

The Decell’s Algorithm 10.1.2 for computing the Moore-Penrose inverse is generalized to polynomial matrices [4]. Consider the polynomial matrix

$$ A(x) = A_0 x^n + A_1 x^{n-1} + \cdots + A_{n-1} x + A_n, $$

where $A_i \in \mathbb {R}^{m \times n}$, $i=0,1,...,n$. Let $B(x) = A(x) A(x)^T$ and

$$\begin{aligned} p(\lambda , x)= & {} \det (\lambda I - B(x)) \nonumber \\= & {} a_0(x) \lambda ^n + a_1 \lambda ^{n-1} + \cdots + a_{n-1} (x) \lambda + a_n(x) , \end{aligned}$$

(10.2.1)

where $a_0(x) = 1$, be the characteristic polynomial of B(x). It is shown in [4] that if k is the largest integer such that $a_k (x) \ne 0$ and $\mathbb {Z}$ is the set containing the zeros of $a_k (x)$, then the Moore-Penrose inverse ${A(x)}^{\dagger }$ of A(x) for $x \in \mathbb {R}{\setminus }\mathbb {Z}$ is given by

$$ {A(x)}^{\dagger } = -{a_k (x)}^{-1} {A(x)}^T \left( {B(x)}^{k-1} + a_1 (x) {B(x)}^{k-2} + \cdots + a_{k-1}(x) I \right) . $$

If $k=0$ is the largest integer such that $a_k (x) \ne 0$, then ${A(x)}^{\dagger } = \mathrm{O}$. Moreover, for each $x_i \in \mathbb {Z}$, if $k_i < k$ is the largest integer such that $a_{k_i} (x_i) \ne 0$, then the Moore-Penrose inverse ${A(x_i)}^{\dagger }$ of $A(x_i)$ is given by

$$\begin{aligned}&{A(x_i)}^{\dagger } \\= & {} -{a_{k_i} (x_i)}^{-1} {A(x_i)}^T \left( {B(x_i)}^{k_i-1} + a_1 (x_i) {B(x_i)}^{k_i-2} + \cdots + a_{k_i-1}(x_i) I \right) . \end{aligned}$$

The algorithm for computing the polynomial matrix ${A(x)}^{\dagger }$ is completely analogous to Algorithm 10.1.2, replacing the scalar matrices A, $A_i$, B, and $P_i$ with the polynomial matrices A(x), $A_i(x)$, B(x), and P(x) respectively and the scalars $a_i$ with the polynomials $a_i(x)$. Obviously, the algorithm involves symbolic computation. Also in [4], a two-dimensional algorithm that avoids symbolic computation is presented. From the definition $B(x) = A(x) {A(x)}^T$, the degree of B(x) can be as high as 2n. Consequently, $A_i(x)$, $a_i(x)$, and $P_{i+1} (x)$ are of degrees up to 2in. Let

$$ a_i (x) = \sum _{j=0}^{2in} a_{i,j} x^j, \quad i=1,...,k, $$

where $a_{i,j}$ are scalars, and

$$ P_i (x) = \sum _{j=0}^{2(i-1)n} P_{i,j} x^j , \quad i=1,...,k, $$

where $P_{i,j}$ are scalar matrices, then ${A(x)}^{\dagger }$ can be written as

$$\begin{aligned} {A(x)}^{\dagger }= & {} -\left( \sum _{j=0}^{2kn} a_{k,j} x^j \right) ^{-1} \left( \sum _{j=0}^{n} A_j^T x_j \right) \left( \sum _{j=0}^{2(k-1)n} P_{k,j} x^j \right) \nonumber \\= & {} -\left( \sum _{j=0}^{2kn} a_{k,j} x^j \right) ^{-1} \left( \sum _{j=0}^{(2k-1)n} \sum _{l=0}^j (A_{j-l}^T P_{k,l}) x^j \right) . \end{aligned}$$

(10.2.2)

Now we derive $a_{i,j}$ and $P_{i,j}$. Following Algorithm 10.1.2, first we have

$$\begin{aligned} A_i (x)= & {} B(x) P_i (x) \\= & {} \left( \sum _{j=0}^{n} A_j x_j \right) \left( \sum _{j=0}^{n} A_j^T x_j \right) \left( \sum _{j=0}^{2(i-1)n} P_{i,j} x^j \right) \\= & {} \left( \sum _{j=0}^{2n} \left( \sum _{p=0}^j A_{j-p} A_p^T \right) x^j \right) \left( \sum _{j=0}^{2(i-1)n} P_{i,j} x^j \right) \\= & {} \sum _{j=0}^{2in} \left( \sum _{p=0}^j \left( \sum _{l=0}^{j-p} (A_{j-p-l} A_l^T) \right) P_{i,p} \right) x^j . \end{aligned}$$

It then follows that

$$\begin{aligned} a_i (x)= & {} -i^{-1} \mathrm{tr}(A_i (x) ) \\= & {} - i^{-1} \sum _{j=0}^{2in} \mathrm{tr} \left( \sum _{p=0}^j \left( \sum _{l=0}^{j-p} (A_{j-p-l} A_l^T) \right) P_{i,p} \right) x^j , \end{aligned}$$

which gives

$$\begin{aligned} a_{i,j} = -\frac{1}{i} \mathrm{tr} \left( \sum _{p=0}^j \left( \sum _{l=0}^{j-p} (A_{j-p-l} A_l^T) \right) P_{i,p} \right) , \quad j=0,...,2in . \end{aligned}$$

(10.2.3)

Moreover,

$$\begin{aligned} P_i (x)= & {} A_{i-1} (x) - a_{i-1} (x) I \\= & {} \sum _{j=0}^{2(i-1)n} \left( \sum _{p=0}^j \left( \sum _{l=0}^{j-p} (A_{j-p-l} A_l^T) \right) P_{i-1,p} \right) x^j + \sum _{j=0}^{2(i-1)n} a_{i-1,j} x^j \\= & {} \sum _{j=0}^{2(i-1)n} \left( \left( \sum _{p=0}^j \left( \sum _{l=0}^{j-p} (A_{j-p-l} A_l^T) P_{i-1,p} \right) \right) + a_{i-1,j} I \right) x^j , \end{aligned}$$

which gives

$$\begin{aligned} P_{i,j} = \left( \sum _{p=0}^j \left( \sum _{l=0}^{j-p} (A_{j-p-l} A_l^T) P_{i-1,p} \right) \right) + a_{i-1,j} I, \end{aligned}$$

(10.2.4)

for $j=0,...,2(i-1)n$.

Finally, we have the following two-dimensional algorithm for computing ${A(x)}^{\dagger }$.

Algorithm 10.2.1

[4] Given a polynomial matrix $A(x) \in \mathbb {R}^{m \times n}[x]$ and k, the largest integer such that $a_k (x)$ in (10.2.1) is nonzero, this algorithm computes the Moore-Penrose inverse of A(x).

• 1. if $k=0$ return ${A(x)}^{\dagger } = \mathrm{O}$;

• 2. $P_{0,0} = \mathrm{O}$; $a_{0,0} = 1$;

• 3. for $i=1$ to k

  • compute $P_{i,j}, \quad j=0,...,2(i-1)n$, by (10.2.4);

  • compute $a_{i,j}, \quad j=0,...,2in$, by (10.2.3);

• 4. compute ${A(x)}^{\dagger }$ by (10.2.2).

Note that in the above algorithm it is assumed that $P_{i,j} = \mathrm{O}$ when $j > 2(i-1)n$. This algorithm is called two-dimensional, since it involves the computation of two-dimensional variables $a_{i,j}$ and $P_{i,j}$.

10.3 Drazin Inverse of a Polynomial Matrix

The definition of the Drazin inverse ${A(x)}_d$ of a polynomial matrix $A(x) \in \mathbb {R}^{n \times n}[x]$ is defined as the same as the scalar case, that is, ${A(x)}_d$ is the matrix satisfying the three conditions:

$$\begin{aligned} {A(x)}^{k+1} {A(x)}_d= & {} {A(x)}^k , \\ {A(x)}_d A(x) {A(x)}_d= & {} {A(x)}_d, \\ A(x) {A(x)}_d= & {} {A(x)}_d A(x) , \end{aligned}$$

where $k = \mathrm{ind}(A(x))$, the index of A(x), defined as the smallest integer such that $\mathrm{rank}({A(x)}^k) = \mathrm{rank}({A(x)}^{k+1})$.

Let

$$\begin{aligned} p(\lambda , x)= & {} \det (\lambda I - A(x)) \\= & {} a_0(x) \lambda ^n + a_1(x) \lambda ^{n-1} + \cdots + a_{n-1}(x) \lambda + a_n (x) , \end{aligned}$$

where $a_0 (x) = 1$, be the characteristic polynomial of ${A(x)}^{k+1}$, then

$$ a_{r+1}(x) = ... = a_n(x) = 0, \quad \hbox {and} \quad a_r (x) \ne 0, $$

where $r = \mathrm{rank}({A(x)}^{k+1})$. Let $\mathbb {Z}$ be the set containing the zeros of $a_r(x)$, then the Drazin inverse ${A(x)}_d$ of A(x) for $x \in \mathbb {R}{\setminus }\mathbb {Z}$ is given by

$$\begin{aligned} {A(x)}_d= & {} -a_r(x)^{-1} {A(x)}^k ( ({A(x)}^{k+1})^{r-1} + a_1(x) ({A(x)}^{k+1})^{r-2} + \cdots \\&+\, a_{r-2}(x) {A(x)}^{k+1} + a_{r-1}(x) I ) . \end{aligned}$$

If $r = 0$, then ${A(x)}_d = \mathrm{O}$.

Following the Decell’s Algorithm 10.1.2, we have the following finite algorithm for computing the Drazin inverse ${A(x)}_d$ of A(x) [5, 6].

Algorithm 10.3.1

[5] Given $A(x) \in \mathbb {R}^{n \times n}[x]$, $k = \mathrm{Ind} (A(x))$ and $r = \mathrm{rank} ({A(x)}^{k+1})$, this algorithm computes the Drazin inverse ${A(x)}_d$ of A(x).

• 1. if $r=0$ return ${A(x)}_d = \mathrm{O}$;

• 2. $B(x) = {A(x)}^{k+1}$;

• 3. $P_0(x) = \mathrm{O}$; $a_0(x) = 1$;

• 4. for $i=1$ to r

  • $P_i(x) = B(x)P_{i-1}(x) + a_{i-1}(x) I$;

  • $a_i = -\mathrm{tr}(B(x)P_i(x)) / i$;

• 5. ${A(x)}_d = -a_r(x)^{-1} {A(x)}^k P_r (x)$.

Following derivation of the two-dimensional Algorithm 10.2.1, we can obtain a two-dimensional algorithm for computing the Drazin inverse of a polynomial matrix [5, 6].

Notice that the degrees of $B(x) = {A(x)}^{k+1}$, ${A(x)}^k$, $P_i(x)$, and $a_i(x)$ are respectively $(k+1)n$, kn, $(i-1)(k+1)n$, and $i(k+1)n$. Let

$$ B(x) = \sum _{j=0}^{(k+1)n} B_j x^j, \quad {A(x)}^k = \sum _{j=0}^{kn} \widehat{A_j} x^j, \quad P_i (x) = \sum _{j=0}^{(i-1)(k+1)n} P_{i,j} x^j, $$

where $B_j$, $\widehat{A_j}$, and $P_{i,j}$ are scalar matrices, and

$$ a_i (x) = \sum _{j=0}^{i(k+1)n} a_{i,j} x^j , $$

where $a_{i,j}$ are scalars, then

$$\begin{aligned}&{A(x)}_d \nonumber \\= & {} -\left( \sum _{j=0}^{r(k+1)n} a_{r,j} x^j \right) ^{-1} \left( \sum _{j=0}^{kn} \widehat{A_j} x^j \right) \left( \sum _{j=0}^{(r-1)(k+1)n} P_{r,j} x^j \right) \nonumber \\= & {} -\left( \sum _{j=0}^{r(k+1)n} a_{r,j} x^j \right) ^{-1} \left( \sum _{j=0}^{(kr+r-1)n} \left( \sum _{l=0}^j (\widehat{A_{j-l}} P_{r,l}) \right) x^j \right) . \end{aligned}$$

(10.3.1)

Now we derive $a_{i,j}$ and $P_{i,j}$. Firstly,

$$\begin{aligned} a_i (x)= & {} -i^{-1} \mathrm{tr}(B(x)P_i(x) ) \\= & {} - i^{-1} \mathrm{tr} \left( \left( \sum _{j=0}^{(k+1)n} B_j x^j \right) \left( \sum _{j=0}^{(i-1)(k+1)n} P_{i,j}(x) x^j \right) \right) \\= & {} - i^{-1} \mathrm{tr} \left( \sum _{j=0}^{i(k+1)n} \left( \sum _{l=0}^j B_{j-l} P_{i,l} \right) x^j \right) , \end{aligned}$$

implying that

$$\begin{aligned} a_{i,j} = -\frac{1}{i} \mathrm{tr} \left( \sum _{l=0}^j B_{j-l} P_{i,l} \right) . \quad j=0,...,i(k+1)n . \end{aligned}$$

(10.3.2)

Secondly,

$$\begin{aligned} P_i (x)= & {} B(x) P_{i-1} (x) - a_{i-1} (x) I \\= & {} \left( \sum _{j=0}^{(k+1)n} B_j x^j \right) \left( \sum _{j=0}^{(i-2)(k+1)n} P_{i-1,j} x^j \right) + \sum _{j=0}^{(i-1)(k+1)n} a_{i-1,j} I x^j \\= & {} \sum _{j=0}^{(i-1)(k+1)n} \left( \sum _{l=0}^j B_{j-l} P_{i-1,l} \right) x^j + \sum _{j=0}^{(i-1)(k+1)n} a_{i-1,j} x^j \\= & {} \sum _{j=0}^{(i-1)(k+1)n} \left( \left( \sum _{l=0}^j B_{j-l} P_{i-1,l} \right) a_{i-1,j} I \right) x^j , \end{aligned}$$

implying that

$$\begin{aligned} P_{i,j} = \left( \sum _{l=0}^j B_{j-l} P_{i-1,l} \right) + a_{i-1,j} I, \quad j=0,...,(i-1)(k+1)n. \end{aligned}$$

(10.3.3)

Finally, we have the following two-dimensional algorithm for computing the Drazin inverse ${A(x)}_d$ of A(x). Comparing with Algorithm 10.3.1, this algorithm avoids symbolic computation.

Algorithm 10.3.2

[5] Given $A(x) \in \mathbb {R}^{n \times n}[x]$, $k = \mathrm{Ind} (A(x))$ and $r = \mathrm{rank} ({A(x)}^{k+1})$, this algorithm computes the Drazin inverse ${A(x)}_d$ of A(x).

• 1. if $k=0$ return ${A(x)}_d = \mathrm{O}$;

• 2. $P_{0,0} = \mathrm{O}$; $a_{0,0} = 1$;

• 3. for $i=1$ to k

  • compute $P_{i,j}, \quad j=0,...,(i-1)(k+1)n$, by (10.3.3);

  • compute $a_{i,j}, \quad j=0,...,i(k+1)n$, by (10.3.2);

• 4. compute ${A(x)}_d$ by (10.3.1).

In the algorithm, it is assumed that $B_j = \mathrm{O}$, when $j > (k+1)n$, and $P_{i,j} = \mathrm{O}$, when $j > (i-1)(k+1)n$.

Example 10.3.1

[5] Let

$$ A(x) = \left[ \begin{array}{cccc} x-1 &{} \;1 &{} \;0 &{} \;0 \\ 0 &{} \;1 &{} \;x &{} \;0 \\ 0 &{} \;0 &{} \;0 &{} \;x \\ 0 &{} \;0 &{} \;0 &{} \;0 \end{array} \right] , $$

then $n=1$.

It can be determined

$$ \mathrm{rank} ({A(x)}^2) = \mathrm{rank} ({A(x)}^3) = 2, \quad \hbox {when} \ x \ne 1 . $$

Thus $r = k = 2$. Initially, we have

$$ P_{0,0} = \mathrm{O} \quad \hbox {and} \quad a_{0,0} = 1. $$

When $i=1$, $P_{1,0} = I$ and

$$\begin{aligned} a_{1,0}= & {} -\mathrm{tr} (B_0 P_{1,0}) = 0, \\ a_{1,1}= & {} -\mathrm{tr} (B_1 P_{1,0}) = -3, \\ a_{1,2}= & {} -\mathrm{tr} (B_2 P_{1,0}) = 3, \\ a_{1,3}= & {} -\mathrm{tr} (B_3 P_{1,0}) = -1. \end{aligned}$$

When $i=2$,

$$\begin{aligned} P_{2,0}= & {} B_0 P_{1,0} + a_{1,0} I = \left[ \begin{array}{rrrr} -1 &{} \;1 &{} \;0 &{} \;0 \\ 0 &{} \;1 &{} \;0 &{} \;0 \\ 0 &{} \;0 &{} \;0 &{} \;0 \\ 0 &{} \;0 &{} \;0 &{} \;0 \end{array} \right] \\ P_{2,1}= & {} B_1 P_{1,0} + a_{1,1} I = \left[ \begin{array}{rrrr} 0 &{} -1 &{} 0 &{} 0 \\ 0 &{} -3 &{} 1 &{} 0 \\ 0 &{} 0 &{} -3 &{} 0 \\ 0 &{} 0 &{} 0 &{} -3 \end{array} \right] \\ P_{2,2}= & {} B_2 P_{1,0} + a_{1,2} I = \left[ \begin{array}{rrrr} 0 &{} \;1 &{} \;1 &{} \;1 \\ 0 &{} \;3 &{} \;0 &{} \;1 \\ 0 &{} \;0 &{} \;3 &{} \;0 \\ 0 &{} \;0 &{} \;0 &{} \;3 \end{array} \right] \\ P_{2,3}= & {} B_3 P_{1,0} + a_{1,3} I = \left[ \begin{array}{rrrr} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} -1 &{} 0 &{} 0 \\ 0 &{} 0 &{} -1 &{} 0 \\ 0 &{} 0 &{} 0 &{} -1 \end{array} \right] \end{aligned}$$

and

$$\begin{aligned} a_{2,0}= & {} -\frac{1}{2}\mathrm{tr}(B_0 P_{2,0}) = -1 \\ a_{2,1}= & {} -\frac{1}{2}\mathrm{tr} (B_1 P_{2,0} + B_0 P_{2,1}) = 3 \\ a_{2,2}= & {} -\frac{1}{2}\mathrm{tr} (B_2 P_{2,0} + B_1 P_{2,1} + B_0 P_{2,2}) = -3 \\ a_{2,3}= & {} -\frac{1}{2}\mathrm{tr} (B_3 P_{2,0} + B_2 P_{2,1} + B_1 P_{2,2} + B_0 P_{2,3}) = 1 \\ a_{2,4}= & {} -\frac{1}{2}\mathrm{tr} (B_3 P_{2,1} + B_2 P_{2,2} + B_1 P_{2,3}) = 0 \\ a_{2,5}= & {} -\frac{1}{2}\mathrm{tr} (B_3 P_{2,2} + B_2 P_{2,3}) = 0 \\ a_{2,6}= & {} -\frac{1}{2}\mathrm{tr} (B_3 P_{2,3}) = 0 \end{aligned}$$

Finally, we obtain

$$ {A(x)}_d = -\frac{1}{(x-1)^3} \left[ \begin{array}{cccc} -(x-1)^2 &{} (x-1)^2 &{} x^2 (x-1) &{} x^2(x^2 - x + 1) \\ 0 &{} -(x-1)^3 &{} -x(x-1)^3 &{} -x^2(x-1)^3 \\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \end{array} \right] , $$

for $x \ne 1$. The case when $x=1$ can be dealt with as a special case. $\square$

Remarks

An algorithm for computing the Moore-Penrose inverse of a polynomial matrix with two variables is presented in [4]. In [7], Karampetakis and Tzekis improved Algorithm 10.2.1 for the case when there are big gaps between the powers of x, for example, $A(x) = A_0 x^{80} + A_{79} x + A_{80}$. The above algorithms can be generalized to rational matrices by using the least common denominator of $P_i (x)$ [6].