Factorizing an 8×8 unitary matrix into tensor product of three 2×2 unitaries I am an undergraduate physics student currently in the final stages of a BSc project. I am trying to decompose an $8 \times 8$ unitary matrix into a tensor product of three $2 \times 2$ unitaries. This is in the context of quantum information where we need to rotate measurement operators for a 3 qubit system but are limited to performing local (ie $2 \times 2$) unitary rotations.
Concretely, the problem we have is as follows:
We have an 'unphysical' unitary rotation $U_{8 \times 8}^{unphysical}$ that needs to be performed. It is unphysical as it is not a product of single qubit ($2 \times 2$) unitary matrices.
So we have tried to numerically (in python) generate $$U_{8 \times 8}^{physical} = U_{2 \times 2}^{(1)} \otimes U_{2 \times 2}^{(2)} \otimes U_{2 \times 2}^{(3)}$$ 
where we guess the $4 \times 3$ (since $k = 1,2,3$) parameters $\theta^{(k)}, \phi^{(k)}_1, \phi^{(k)}_2, \phi^{(k)}$ in the single qubit unitaries $$U_{2 \times 2}^{(k)} = e^{i \phi^{(k)}} \left( \begin{array}{cc}
  e^{i \phi^{(k)}_1} \cos(\theta^{(k)})& e^{i \phi^{(k)}_2} \sin(\theta^{(k)})\\
  -e^{-i \phi^{(k)}_2}\sin(\theta^{(k)}) & e^{-i \phi^{(k)}_1} \cos(\theta^{(k)})\\
 \end{array} \right)$$
I have a suspicion that it may not be possible to obtain a factorization that produces $U_{8 \times 8}^{physical} = U_{8 \times 8}^{unphysical}$ but I would like to get as close as possible. 
Our current cost function is simply the Frobenius norm of of $U_{8 \times 8}^{physical} - U_{8 \times 8}^{unphysical}$ where every time the solver guesses the parameters, the cost function builds the corresponding single qubit unitaries, tensors (kronecker products) them together as in the construction above and calculates this 'distance'.
The black box solver (scipy.optimize.minimize) stops after a suspiciously short number of iterations (we have used it throughout the project and it has been robust) and finishes with the cost function being about 2 or 3, which is clearly not good enough.
We can either come up with a better cost function and be able to decompose $U_{8 \times 8}^{unphysical}$ numerically, or there may be a mathematical theorem/construct that can help us do this analytically.
I am happy to provide more detail on the context or clarify the problem we have.
Thank you,
Alex
 A: This is not an exact answer to this question. But, it may be beneficial in this context.
In terms of quantum mechanics, a unitary operator is local if it can be expressed as a tensor product of smaller unitary matrices of suitable order. Otherwise, the operator is a global unitary operator, as mentioned in the comment. We use matrix alignment method to classify a unitary operator. An outline of the method may be as follows. Partition the unitary matrix $U$ into blocks. Write all the elements of the block as a row vector of a new matrix $\overline{U}$. If rank of $\overline{U} = 1$ then $U$ is a local unitary operator, and you may write it as a tensor product of smaller matrices. But, how to calculate the smaller matrices is not still clear in the literature. 
Details of this method is discussed in https://arxiv.org/pdf/1310.2042.pdf.
A: Tensor product factorizations of any particular shape, when they exist, are actually really easy to compute. For example,
$$ \begin{pmatrix} a & b \\ c & d \end{pmatrix} \otimes B
= \begin{pmatrix} a B & b B \\ c B & d B \end{pmatrix} $$
So, you just pick some specific coordinate of the blocks on the right — e.g. the top-left entry of each block — and that tells you exactly what the entries are of the left factor. (up to an arbitrary choice of scale)
(this answer does not address the more general question where the shape of the factorization is allowed to vary — e.g. the case where a factorization into Kronecker products only exists after an unknown change of basis of a constrained type)
