How to compute the operator Schmidt decomposition of a matrix? Given a matrix $A \in \mathcal{L}(\mathcal{H}_1\otimes\mathcal{H_2})$ is there always an operator Schmidt decomposition like
$$
A = \sum_{i=0}^n \lambda_i B_i \otimes C_i,
$$
with $\lambda_i \geq 0$ and $B_i \in \mathcal{H}_1, C_i \in \mathcal{H}_2$ ?
If yes, how can one (numerically) compute the decomposition with $A$ given? What role does the singular value decomposition (SVD) play here?
Could you illustrate the method for the matrix
$$
A = \sigma_x \otimes \sigma_z = \begin{pmatrix} 0&0&1&0\\0&0&0&-1\\1&0&0&0\\0&-1&0&0 \end{pmatrix}?
$$
 A: Concretely, given an object/tensor with any number of indices, you can choose any bipartition of the indices and perform the corresponding SVD.
More precisely, this means to think of the object as an operator acting between two linear spaces, and perform its SVD as you would normally.
For example, your $A$ is an operator acting in a bipartite space, and therefore has $4$ indices. Write these as $A=\sum_{ijk\ell}A_{ijk\ell} |ij\rangle\!\langle k\ell|$. The "standard SVD" corresponds to thinking of this object as the linear operator
$$|k\ell\rangle\mapsto \sum_{ij} A_{ijk\ell} |ij\rangle,$$
and doing the SVD of the corresponding matrix. However, you could also consider the linear operator
$$|j\ell\rangle\mapsto \sum_{ik} A_{ijk\ell} |ik\rangle.$$
If you perform the SVD of the corresponding matrix (which is obtained from the other one transposing some layers), you get a decomposition of the form
$$A_{ijk\ell} = \sum_m s_m u^m_{ik} \,\bar v^m_{j\ell},$$
for some $s_m\ge0$ and collections of orthonormal vectors $\{u^m\}_m,\{v^m\}_m$. If you use this decomposition in the original formula for $A$, you get
$$A = \sum_{ijk\ell} \sum_m s_m u^m_{ik}\,\bar v^m_{j\ell} |ij\rangle\!\langle k\ell|
= \sum_m s_m B_m\otimes C_m,$$
where
$$B_m\equiv \sum_{ik} u^m_{ik} |i\rangle\!\langle k|,
\qquad
C_m\equiv \sum_{ik} \bar v^m_{j\ell} |j\rangle\!\langle \ell|.$$
Note that the orthonormality of the vectors $u^m,v^m$, given by the SVD, translate into the orthonormality of the operators $B_m,C_m$ with respect to the Hilbert-Schmidt inner product:
$$\operatorname{Tr}(B_m^\dagger B_n)=\operatorname{Tr}(C^\dagger_m C_n)=\delta_{mn}.$$
Now, if you want to do this operator decomposition in practice, given an operator expressed as a matrix, you need to "reshape" the matrix to obtain a four-index tensor object, then convert this into a new matrix, and compute the SVD of that. If you have access to Mathematica, here's a quick snippet to take a matrix and return the corresponding operator Schmidt decomposition:
reshapeAndComputeSVD[matrix_, dims_List] := ArrayReshape[matrix, dims] // Flatten[#,
    {{1, 3}, {2, 4}}
] & // SingularValueDecomposition;

visualizeOperatorSchmidtDecomposition[matrix_, dims_List] := reshapeAndComputeSVD[matrix, dims] // Thread[
        {Transpose @ #[[1]], Diagonal @ #[[2]], Transpose @ #[[3]]}
    ]& // Map[
        #[[2]] * Inactive[CircleTimes][
            ArrayReshape[#[[1]], dims[[{1, 3}]]] // MatrixForm,
            ArrayReshape[#[[3]], dims[[{2, 4}]]] // MatrixForm
        ]
    &] // Total;

and here's an example of what you get in a few simple cases:

