If the goal would be to solve the (matrix) equation: $${\bf P}^2 = {\bf A}$$ Do you think a numerical scheme alternatingly minimizing for $\bf P_1, P_2$ would be stable: $$\|{\bf P_1P_2-A}\|_2,\hspace{1cm} \|{\bf P_1-P_2}\|_2$$
I do not expect there to be any unique solution. On the contrary I am especially interested in the cases when there aren't, as the approach would be very convenient in providing ways of incorporating constraints on a solution ( if we can get it to work ).
On
Yes this approach proves to have some merit. To take a elementary example, we take maybe the simplest case $\bf A=I$ or in other words an operation which is it's own inverse - an involution. We choose to conduct our experiments over $\mathbb C^{6}$ and starting with random matrices $\bf P_1$ and $\bf P_2$.
After 6 iterations, we calculate eigenvalues and eigenvectors of $\bf P_1$ and we do get:
$${\bf D} = \left[\begin{array}{rrrrrr}-1&0&0&0&0&0\\0&-1&0&0&0&0\\0&0&1&0&0&0\\0&0&0&1&0&0\\0&0&0&0&1&0\\0&0&0&0&0&1\end{array}\right] \pm {\bf I}\cdot10^{-4}$$ which is completely reasonable as the set of solutions to $x^2=1$ over our field of scalars is $x\in \{-1,1\}$
We can then keep going restarting on the solution found ($\bf A = P_1$) and iterating 12 times:
$${\bf D} = \left[\begin{array}{rrrrrr}1&0&0&0&0&0\\0&-i&0&0&0&0\\0&0&i&0&0&0\\0&0&0&-1&0&0\\0&0&0&0&-1&0\\0&0&0&0&0&-1\end{array}\right] \pm {\bf I}\cdot 10^{-4}$$
which is also reasonable since the roots to $x^4 = 1$ are $x \in \{-1,1,-i,i\}$ and if we look carefully we see that $x^2 = -1$ have exactly two roots stemming from the two $-1$ above : $\{i,-i\}$, the other four, $\{1,-1,-1,-1\}$ stem from the 1 eigenvalues.
If this works well there is in theory no stopping us from doing a kind of iterative refinement solving ${\bf P}^{2^k} = {\bf A}$ by first solving ${\bf P}^2 = {\bf A}$ and then setting ${\bf A = P}$ and iterating $k$ times.
You can use iterative methods. Babylonian method, can be an option.