Consider the matrix equation $${\bf X}^k\bf {A = B}$$ Which we want to solve for $\bf X$
We can put A and B in a "block-vector": $v = [{\bf A}^T,{\bf 0},\cdots,{\bf 0},{\bf B}^T]^T$, assume there exists a matrix: $$M = \begin{bmatrix}\bf I&\bf 0&\bf 0&\cdots&\bf 0\\\bf X_1&\bf 0&\bf 0&\cdots&\bf 0\\\bf 0&\bf X_2&\bf 0&\cdots&\bf 0\\\vdots&\ddots&\ddots&\cdots&\bf 0\\\bf 0&\bf 0&\bf 0&\bf X_k&\bf 0\end{bmatrix}$$ With $k$ nr of $X$s in it. $$\|{\bf C(Mv-d)}\|_F^2$$ with $\bf C$ (like a boundary condition) encoding that fit to first and last element is important and disregard what happens in between.
Finally solve wrt $\bf X$, for example by vectorization, Kronecker products and a new term regularizing to get the $\bf 0$s to $\bf 0$ and a set of equation systems tying the first $\bf X$ to the second, et.c. $+\|{\bf X_i-X_{i+1}}\|_F^2$. When solved, if successful the ${\bf X_i}$ should be so close to each other that we can take for example the mean value of them as the solution for ${\bf X}$.
Would this be a reasonable approach of linearizing the equation or am I missing something important?
EDIT I found a potential brain fart. If expressed as equations, I think we need the diagonal strip of $\bf -I$ below, but I am not sure yet.
$$M = \begin{bmatrix}\bf I&\bf 0&\bf 0&\cdots&\bf 0\\\bf X_1&\bf -I&\bf 0&\cdots&\bf 0\\\bf 0&\bf X_2&\bf -I&\cdots&\bf 0\\\vdots&\ddots&\ddots&\cdots&\bf 0\\\bf 0&\bf 0&\bf 0&\bf X_k&\bf -I\end{bmatrix}$$
NEW EDIT None of the above two seem to work very well. I do get some matrix as a solution so the $$+\displaystyle \sum_{i=1}^{k-1} \|{\bf X_i - X_{i+1}}\|_2^2$$ seems to work as intended, but any of the $\bf X_l$ is several orders of magnitude too small to make any sense. Still I can't find any theoretical reason for why it should not work so my current assumption is that I have made a mistake in the implementation. Another possibility is that maybe it is in general quite hard to even find one solution to such polynomial matrix equations so that some kind of regularization that I have not thought about would be required.
Newest edit! I realized I need one more line in the matrices and one more $\bf 0$ block in the right hand side vector, like this: $$M = \begin{bmatrix}\bf I&\bf 0&\bf 0&\cdots&\bf 0\\\bf X_1&\bf -I&\bf 0&\cdots&\bf 0\\\bf 0&\bf X_2&\bf -I&\cdots&\bf 0\\\vdots&\ddots&\ddots&\cdots&\bf 0\\\bf 0&\bf 0&\bf 0&\bf X_k&\bf -I\\\bf 0 &\cdots& \bf 0& \bf 0 & \bf I\end{bmatrix}, v = \begin{bmatrix}\bf A\\\bf 0\\\bf 0\\ \vdots\\\bf 0\\\bf B\end{bmatrix}$$
This way both $\bf A,B$ get their very own equations for fixation a.k.a. "boundary conditions". Without any similar fixation for $\bf B$ as we had for $\bf A$ before there is no guarantee what properties our $\bf X_k$ should be getting when solving for $M$.