Closed form for the inverse of a special transition matrix Let $b\in\mathbf{R}^n$ be fixed, and let $A\in\mathcal{M}_{n\times n}(\mathbf{R})$ be given but otherwise arbitrary. Let $a_1,\ldots,a_n$ denote the coefficients of the characteristic polynomial of $A$
\begin{equation}
\chi_A(x) = x^n + a_1 x^{n-1} + \ldots + a_{n-1} x + a_n.
\end{equation}
We shall assume that $\text{span}\big\{b, Ab, \ldots, A^{n-1}b\big\}=\mathbf{R}^n$, and let us consider the invertible matrix
\begin{equation}
P = [f_1 \,\, \ldots \,\, f_n]
\end{equation}
where the columns $f_k$ are of the form
\begin{equation}
f_n = b \quad \text{ and } \quad f_k = \left(A^{n-k}+\sum_{j=1}^{n-k}a_j A^{n-k-j}\right)b \quad \text{ for } 1 \leqslant k \leqslant n-1.
\end{equation}
My question is:

Is there a closed form for $P^{-1}$ where the dependence on $b$ is explicitly seen/tracked?

When $A$ is nilpotent, then we can write the inverse as a finite Neumann series, so I would like to avoid this case and consider more general matrices $A$, possibly invertible as well (e.g. tridiagonal Toeplitz, and so on).
Note that $P=P(b)$ can be rewritten as the product of the two rectangular matrices
\begin{equation}
P = [p_1(A) \,\, \ldots\,\, p_n(A)] \begin{bmatrix} b& & \\ & \ddots & \\ & & b\end{bmatrix},
\end{equation}
where $p_j(A):=A^{n-k}+\sum_{j=1}^{n-k}a_j A^{n-k-j}$ for $j\leqslant n-1$ and $p_n=\text{Id}$, but I am not sure if this could be of any use in further computations.
 A: Let $B$ be a matrix formed by the columns $[b|Ab|\dots|A^{n-1}b]$. Since $A$ is a zero of $\chi_A(x)$, the matrix $B$ satisfies the identity:
$$AB = BC,$$
where
$$C:=\begin{bmatrix} 
0 & 0 & 0 & \dots & 0 & -a_n\\
1 & 0 & 0 & \dots & 0 & -a_{n-1}\\
0 & 1 & 0 & \dots & 0 & -a_{n-2}\\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & 0 & \dots & 0 & -a_2\\
0 & 0 & 0 & \dots & 1 & -a_1
\end{bmatrix}
$$
is the companion matrix of $\chi_A(x)$.
From the definition of $P$, it follows that
$$P = BD$$
with a Hankel matrix
$$D:=\begin{bmatrix} 
a_{n-1} & a_{n-2} & a_{n-3} & \dots & a_1 & 1\\
a_{n-2} & a_{n-3} & a_{n-4} & \dots & 1 & 0\\
a_{n-3} & a_{n-4} & a_{n-5} & \dots & 0 & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
a_1 & 1 & 0 & \dots & 0 & 0\\
1 & 0 & 0 & \dots & 0 & 0
\end{bmatrix}.
$$
It can be easily seen that both $D$ and $CD$ are symmetric matrices, implying that
$D^{-1}CD = C^T$.
Combining all together, we have $AP = PC^T$ and thus
$$P^{-1}A = C^TP^{-1},$$
which is a Sylvester equation with respect to $P^{-1}$.
The last equation implies that $P^{-1}$ is formed by rows $v^TA^i$ for $i=0,1,\dots,n-1$, where vector $v^T$ (the first row of $P^{-1}$) is also the last row of $B^{-1}$.
