Matrix factorization I'd like to factorize matrices as follows:
$$ \left(\begin{array}{cc}X_1&X_2\\X_3&X_4\end{array}\right) = \left(\begin{array}{cc}D_1&D_2\\D_3&D_4\end{array}\right)\left(\begin{array}{cc}P_1&\\&P_2\end{array}\right)\left(\begin{array}{cc}D_5&D_6\\D_7&D_8\end{array}\right) $$ 
provided such a factorization exists. 
The left-hand side is known and the right-hand side isn't. The $X_i$ are full $2^t\times 2^t$ ($t\ge 1$), the other blocks have same dimension with the $D_i$ being diagonal and the $P_i$ full. 
(edit: the previous notations I used for the diagonal blocks $D_i$ might have been unclear, these blocks are not necessarily scalar multiples of the identity)
Assuming a simple case where the diagonal matrices are invertible, the above problem amounts to solving the following nonlinear system of matrix equations (I use the same notations to simplify notations but the matrices are not the same as the previous one. The reason for this is to keep notations simple esp. given that what matters here is the structure of the blocks.)
$$ \left\{ \begin{array}{lcl} X_1 &=& P_1 + P_2 \\ X_2 &=& P_1 D_1 + P_2 D_2\\ X_3 &=& D_3 P_1 + D_4 P_2\\ X_4 &=& D_3 P_1 D_1 + D_4 P_2 D_2\end{array} \right. $$ 
I've been looking for a way to solve this system but so far without much success. One thing I tried is to fix $D_1$ and $D_2$ (say), obtain the corresponding $P_1$ and $P_2$ with the first two equations and then find the best $D_3$ and $D_4$ in the Frobenius sense using the last two equations. Then start the other way around using $D_3$ and $D_4$. However this does not seem to converge (looks like projections on non-convex sets.). Also, given that there is $4n^2$ ($n:=2^t)$ equations with $2n^2+4n$ unknowns, maybe that this system can be further simplified.
Any insight is most welcome, thanks!
Edit: any result or idea on the potential infeasibility of finding an efficient/elegant way of solving this is also welcome. I've also been looking into the simplest case of this problem where each block is exactly $2\times 2$. In order not to use notations which would make this question even more confusing that it already is, I'll just use bullets to denote potentially non-zero entries since what really matters here is structure.
$$ \left(\begin{array}{cc}  \left(\begin{array}{cc}\bullet&\bullet\\\bullet&\bullet\end{array}\right)   & \left(\begin{array}{cc}\bullet&\bullet\\\bullet&\bullet\end{array}\right)\\\left(\begin{array}{cc}\bullet&\bullet\\\bullet&\bullet\end{array}\right)&\left(\begin{array}{cc}\bullet&\bullet\\\bullet&\bullet\end{array}\right)\end{array}\right) 
= 
\left(\begin{array}{cc}\left(\begin{array}{cc}\bullet&\\&\bullet\end{array}\right)&\left(\begin{array}{cc}\bullet&\\&\bullet\end{array}\right)\\
\left(\begin{array}{cc}\bullet&\\&\bullet\end{array}\right)&\left(\begin{array}{cc}\bullet&\\&\bullet\end{array}\right)\end{array}\right)
\left(\begin{array}{cc}  \left(\begin{array}{cc}\bullet&\bullet\\\bullet&\bullet\end{array}\right) &\\&\left(\begin{array}{cc}\bullet&\bullet\\\bullet&\bullet\end{array}\right)\end{array}\right)
\left(\begin{array}{cc}\left(\begin{array}{cc}\bullet&\\&\bullet\end{array}\right)&\left(\begin{array}{cc}\bullet&\\&\bullet\end{array}\right)\\
\left(\begin{array}{cc}\bullet&\\&\bullet\end{array}\right)&\left(\begin{array}{cc}\bullet&\\&\bullet\end{array}\right)\end{array}\right) $$ 
(lots of bullets...) when applying a perfect shuffle on this (aka bit-reversal) on this system, we get an equivalent system with the following form:
$$ \left(\begin{array}{cc}  \left(\begin{array}{cc}\bullet&\bullet\\\bullet&\bullet\end{array}\right)   & \left(\begin{array}{cc}\bullet&\bullet\\\bullet&\bullet\end{array}\right)\\\left(\begin{array}{cc}\bullet&\bullet\\\bullet&\bullet\end{array}\right)&\left(\begin{array}{cc}\bullet&\bullet\\\bullet&\bullet\end{array}\right)\end{array}\right) 
= 
\left(\begin{array}{cc}  \left(\begin{array}{cc}\bullet&\bullet\\\bullet&\bullet\end{array}\right) &\\&\left(\begin{array}{cc}\bullet&\bullet\\\bullet&\bullet\end{array}\right)\end{array}\right)
\left(\begin{array}{cc}\left(\begin{array}{cc}\bullet&\\&\bullet\end{array}\right)&\left(\begin{array}{cc}\bullet&\\&\bullet\end{array}\right)\\
\left(\begin{array}{cc}\bullet&\\&\bullet\end{array}\right)&\left(\begin{array}{cc}\bullet&\\&\bullet\end{array}\right)\end{array}\right)
\left(\begin{array}{cc}  \left(\begin{array}{cc}\bullet&\bullet\\\bullet&\bullet\end{array}\right) &\\&\left(\begin{array}{cc}\bullet&\bullet\\\bullet&\bullet\end{array}\right)\end{array}\right)
$$
which maybe can be considered in the context of simultaneous diagonalization of matrices? If solving this particular system can be done, I'm hoping it could give hints for the more general case.
Thanks!
 A: Short answer : your problem can be transformed into an equivalent diagonalization problem, except in a few degenerate and easier cases. 
Longer answer :
I’m sorry but I don’t like your notations at all. I will write
$$  
 \left(\begin{array}{cc} {\mathbf{X_1}}& {\mathbf{X_2}}\\ {\mathbf{X_3}}& {\mathbf{X_4}}\end{array}\right)
=\left(\begin{array}{cc} a{\mathbf I}& b{\mathbf I}\\ c{\mathbf I}& d{\mathbf I}\end{array}\right)\left(\begin{array}{cc}{\mathbf P_1}&\\&{\mathbf P_2}\end{array}\right)\left(\begin{array}{cc} e{\mathbf I}& f{\mathbf I}\\ g{\mathbf I}& h{\mathbf I}\end{array}\right)
$$
instead of
$$
\left(\begin{array}{cc}X_1&X_2\\X_3&X_4\end{array}\right) = \left(\begin{array}{cc}D_a&D_b\\D_c&D_d\end{array}\right)\left(\begin{array}{cc}P_1&\\&P_2\end{array}\right)\left(\begin{array}{cc}D'_a&D'_b\\D'_c&D'_d\end{array}\right)
$$ 
We then can write
$$
\begin{array}{l}
{\mathbf X_1}=(ae){\mathbf P_1}+(bg){\mathbf P_2} \\
{\mathbf X_2}=(af){\mathbf P_1}+(bh){\mathbf P_2} \\
{\mathbf X_3}=(ce){\mathbf P_1}+(dg){\mathbf P_2} \\
{\mathbf X_4}=(cf){\mathbf P_1}+(dh){\mathbf P_2} \\
\end{array}
$$
So an obvious necessary condition for your problem to have a solution is that the vector space $\cal X$ spanned by $\mathbf{X_1},\mathbf{X_2},\mathbf{X_3},\mathbf{X_4}$ has dimension at most 2.
If the dimension of $\cal X$ is $0$, all the $\mathbf X_i$ are zero and we are done. If the dimension of $\cal X$ is $1$, all the $\mathbf{X_i}$ are multiples of a single matrix $\mathbf Z$ : ${\mathbf{X_i}}=x_i{\mathbf{Z}}
$ for all $i \in [1,4]$ . Then we have a solution as follows :
$$  \left(\begin{array}{cc} {\mathbf X_1}& {\mathbf X_2}\\ {\mathbf X_3}& {\mathbf X_4}\end{array}\right)
=\left(\begin{array}{cc} x_1{\mathbf I}& x_2{\mathbf I}\\ x_3{\mathbf I}& x_4{\mathbf I}\end{array}\right)\left(\begin{array}{cc}{\mathbf Z}&\\&{\mathbf Z}\end{array}\right)\left(\begin{array}{cc} {\mathbf I}& \\ & {\mathbf I}\end{array}\right) 
$$ 
So we can assume that the dimension of $\cal X$ is exactly $2$.  Then there are two indices $1 \le i <j \le 4$ such that $\lbrace \mathbf{X_i},\mathbf{X_j}\rbrace$ is a basis of $\cal X$. And $\lbrace P_1,P_2 \rbrace$ is another basis of $\cal X$, so there is a $2 \times 2$ matrix (let’s call it $M=\left(\begin{array}{cc}m_{11}&m_{12}\\m_{21}&m_{22}\end{array}\right)
$) such that 
$$
\left(\begin{array}{c} {\mathbf P_1}\\ {\mathbf P_2}\end{array}\right)=M\left(\begin{array}{c} {\mathbf X_i}\\ {\mathbf X_j}\end{array}\right)
$$
Suppose for example that $i=1$ and $j=2$. Then, there are numbers $x_1,x_2,y_1,y_2$ such that
$$
\begin{array}{l}
{\mathbf {X_3}}=x_1 {\mathbf {X_1}} + x_2 {\mathbf {X_2}} \\
{\mathbf {X_4}}=y_1 {\mathbf {X_1}} + y_2 {\mathbf {X_2}}
\end{array}
$$
Then the equality we want can be rewritten as
$$  
 \left(\begin{array}{cc} {\mathbf{X_1}}& {\mathbf{X_2}}\\ x_1 {\mathbf {X_1}} + x_2 {\mathbf {X_2}}
& y_1 {\mathbf {X_1}} + y_2 {\mathbf {X_2}}\end{array}\right)
=\left(\begin{array}{cc} a{\mathbf I}& b{\mathbf I}\\ c{\mathbf I}& d{\mathbf I}\end{array}\right)\left(\begin{array}{cc}m_{11}{\mathbf X_1}+m_{12}{\mathbf X_2}&\\&m_{21}{\mathbf X_1}+m_{22}{\mathbf X_2}
\end{array}\right)\left(\begin{array}{cc} e{\mathbf I}& f{\mathbf I}\\ g{\mathbf I}& h{\mathbf I}\end{array}\right)
$$
Since $\mathbf X_1$ and $\mathbf X_2$ are linearly independent, we may identify coefficients; this yields :
$$  
 \left(\begin{array}{cc}  1 & 0\\ x_1 
& y_1 \end{array}\right)
=\left(\begin{array}{cc} a& b\\ c& d\end{array}\right)\left(\begin{array}{cc}m_{11}&\\&m_{21}
\end{array}\right)\left(\begin{array}{cc} e& f\\ g& h\end{array}\right)
$$
and
$$  
 \left(\begin{array}{cc}  0 & 1\\ x_2 
& y_2 \end{array}\right)
=\left(\begin{array}{cc} a& b\\ c& d\end{array}\right)\left(\begin{array}{cc}m_{12}&\\&m_{22}
\end{array}\right)\left(\begin{array}{cc} e& f\\ g& h\end{array}\right)
$$
If we put
$$
{\mathbf A}=\left(\begin{array}{cc} a& b\\ c& d\end{array}\right) 
\ \ {\text{and}} \ \ {\mathbf E}=\left(\begin{array}{cc} e& f\\ g& h\end{array}\right), 
$$
we deduce
$$  
 \left(\begin{array}{cc}  1 & 0\\ x_1 
& y_1 \end{array}\right)
={\mathbf A}\left(\begin{array}{cc}m_{11}&\\&m_{21}
\end{array}\right){\mathbf E} \ \ {\text{and}} \ \ 
 \left(\begin{array}{cc}  0 & 1\\ x_2 
& y_2 \end{array}\right)
={\mathbf A}\left(\begin{array}{cc}m_{12}&\\&m_{22}
\end{array}\right){\mathbf E}
$$
If $x_1$ and $y_2$ are both nonzero, all the matrices above are invertible (and in particular, all the $m_{ij}$ are nonzero), and we can write
$$
\frac{1}{y_2}\left(\begin{array}{cc}  -y_2 & 1\\ x_2 
& 0 \end{array}\right)
=\left(\begin{array}{cc}  0 & 1\\ x_2 
& y_2 \end{array}\right)^{-1}
={\mathbf E}^{-1}\left(\begin{array}{cc}\frac{1}{m_{12}}&\\&\frac{1}{m_{22}}
\end{array}\right){\mathbf A}^{-1}
$$
and multiplying,
$$
\left(\begin{array}{cc}  1 & 0\\ x_1 
& y_1 \end{array}\right)\frac{1}{y_2}\left(\begin{array}{cc}  -y_2 & 1\\ x_2 
& 0 \end{array}\right)
={\mathbf A}\left(\begin{array}{cc}m_{11}&\\&m_{21}
\end{array}\right){\mathbf E}
{\mathbf E}^{-1}\left(\begin{array}{cc}\frac{1}{m_{12}}&\\&\frac{1}{m_{22}}
\end{array}\right){\mathbf A}^{-1}
$$
which means that
$$
\frac{1}{y_2}
\left(\begin{array}{cc}  -y_2 & 1\\ -x_1y_2+x_2y_1 
& x_1 \end{array}\right)
={\mathbf A}\left(\begin{array}{cc}\frac{m_{11}}{m_{12}}&\\&\frac{m_{21}}{m_{22}}
\end{array}\right){\mathbf A}^{-1}
$$
So the matrix $N=\left(\begin{array}{cc}  -y_2 & 1\\ -x_1y_2+x_2y_1 & x_1 \end{array}\right)$ must be diagonalizable. Conversely, if $N$ is diagonalizable,  the computations above can be reversed to find a solution to the initial problem.
