Let $A_{n,d}$ be the space of polynomials of degree $d$ or less whose coefficients are real $n\times n$ matrices --- or, if you prefer, the space of matrices whose entries are degree-$d$ polynomials. Given $A(z),B(z)\in A_{n,d}$, there are $P(z),Q(z)\in A_{n,d}$ so that $P(z)A(z)=Q(z)B(z)$ (proof: consider the $2dn\times (2d-1)n $ Sylvester resultant matrix of the polynomials; its left nullspace has dimension at least $n$ for dimension reasons, and one can build $P(z)$ and $Q(z)$ out of the entries of a $n\times 2nd$ matrix in this nullspace $[P_0\;P_1\;\dots\;P_d\;Q_0\;Q_1\;\dots\;Q_d]$).
As Gerry Myerson noted in the comments, $P=Q=0$ is ok if we do not impose any additional constraint. The proof above shows how to compute two polynomial such that $[P_0\;P_1\;\dots\;P_d\;Q_0\;Q_1\;\dots\;Q_d]$ is full-rank. We could ask for that as a normalizing condition. In the generic case, the left kernel of the Sylvester matrix has rank $n$, so with this constraint $P$ and $Q$ are well-defined up to left multiplication by a nonsingular matrix.
1) is there a simple and efficient way to compute $P(z)$ and $Q(z)$ by exploiting the structure of the matrix?
Currently my best shot would be Fourier transforms+a fast Cauchy-like solver like the GKO algorithm, for a total cost of $O(n^3d^2)$ and a rather cumbersome algorithm (I know about superfast solvers, but they seems to be unstable and convenient only for large input sizes). But this seems excessively complicated and costly. For the scalar case, one can block-triangularize the resultant matrix by running the Euclidean algorithm, but in the matrix case this does not work, as a nonzero leading coefficient need not be invertible. Maybe there is an easy FFT-based solution (with cost $O(n^3d\log^{something}d)$) that I am missing?
2) can we say more about the properties of $P(z)$ and $Q(z)$, for example that their determinants do not vanish for any $z$ whenever those of $A(z)$ and $B(z)$ don't?
