Have you tried checking the condition number of the three matrices that you invert in the RHS, instead of the LHS? I am afraid that's the only sure way to tell.
I've heard in many conferences that many experts in matrix computations (led by the late Gene Golub) consider the SMW identity to be "numerically dangerous"; I am not sure that I can find a paper where this is stated black on white, though (sorry for that, I was always curious about that myself).
Alternative idea: switch to solving with the extended matrix $\begin{bmatrix}A & U\\\\ U^* & B\end{bmatrix}$ --- this could be stabler.
The following is just a minor variation of Martin Argerami's proof of the old question. I am even copying his equations and some of his text. If you are +1ing this post, please also +1 his one (if not already done).
Here is a counterexample for $n=2$. Let $\varepsilon $ be a positive real $<
\dfrac{1}{2}$. The matrix
$
a= \begin{bmatrix}
1 & 0 & 0 & 1-\varepsilon \\
0 & \varepsilon & 0 & 0 \\
0 & 0 & \varepsilon & 0 \\
1-\varepsilon & 0 & 0 & 1
\end{bmatrix}
\in \mathrm{M}_{4}\left( \mathbb{C}\right)
$
is positive-definite, but it cannot be written as a sum of tensor products
of nonnegative-semidefinite $2\times 2$-matrices. Here is why:
Assume the contrary. Thus, $a$ is written in the form
$a=\sum_{j} \left[ \begin{matrix}
\alpha _{j} & \overline{\gamma _{j}} \\
\gamma _{j} & \beta _{j}
\end{matrix} \right]
\otimes \left[
\begin{matrix}
\alpha _{j}^{\prime } & \overline{\gamma _{j}^{\prime }} \\
\gamma _{j}^{\prime } & \beta _{j}^{\prime }
\end{matrix} \right]
=\left[
\begin{matrix}
\sum_{j}\alpha _{j}^{\prime }\alpha _{j} & \sum_{j}\alpha _{j}^{\prime }
\overline{\gamma _{j}} & \sum_{j}\overline{\gamma _{j}^{\prime }}\alpha _{j}
& \sum_{j}\overline{\gamma _{j}^{\prime }}\gamma _{j} \\
\sum_{j}\alpha _{j}^{\prime }\gamma _{j} & \sum_{j}\alpha _{j}^{\prime
}\beta _{j} & \ast & \ast \\
\ast & \ast & \sum_{j}\beta _{j}^{\prime }\alpha _{j} & \ast \\
\ast & \ast & \ast & \ast
\end{matrix}\right] $,
where $j$ ranges from $1$ to some positive integer $N$. Since each $
\begin{bmatrix}
\alpha _{j} & \overline{\gamma _{j}} \\
\gamma _{j} & \beta _{j}
\end{bmatrix}
$ is nonnegative-semidefinite, we have $\alpha _{j}\geq 0$, $\beta _{j}\geq 0$,
and $\alpha _{j}\beta _{j}\geq \left\vert \gamma _{j}\right\vert ^{2}$ for
all $j$. Similarly, $\alpha _{j}^{\prime }\geq 0$, $\beta _{j}^{\prime }\geq
0$, and $\alpha _{j}^{\prime }\beta _{j}^{\prime }\geq \left\vert \gamma
_{j}^{\prime }\right\vert ^{2}$ for all $j$.
Now, comparing the entries of $a$ in this equation, we get $\varepsilon
=\sum_{j}\alpha _{j}^{\prime }\beta _{j}$ (from the $\left( 2,2\right) $-th
entry) and $\varepsilon =\sum_{j}\beta _{j}^{\prime }\alpha _{j}$ (from the $
\left( 3,3\right) $-th entry). Taking the arithmetic mean of these two
equations, we get
$\varepsilon =\dfrac{1}{2}\left( \sum_{j}\alpha
_{j}^{\prime }\beta _{j}+\sum_{j}\beta _{j}^{\prime }\alpha _{j}\right)
=\sum_{j}\dfrac{\alpha _{j}^{\prime }\beta _{j}+\beta _{j}^{\prime }\alpha
_{j}}{2}\geq \sum_{j}\sqrt{\alpha _{j}^{\prime }\beta _{j}\beta _{j}^{\prime
}\alpha _{j}}$
(by AM-GM, since we are dealing with nonnegative reals). But
for every $j$, we have
$\sqrt{\alpha _{j}^{\prime }\beta _{j}\beta
_{j}^{\prime }\alpha _{j}}=\sqrt{\alpha _{j}\beta _{j}}\sqrt{\alpha
_{j}^{\prime }\beta _{j}^{\prime }}\geq \left\vert \gamma _{j}\right\vert
\left\vert \gamma _{j}^{\prime }\right\vert $
(since $\alpha _{j}\beta_{j}\geq \left\vert \gamma _{j}\right\vert ^{2}$ and $\alpha _{j}^{\prime }\beta _{j}^{\prime }\geq \left\vert \gamma _{j}^{\prime }\right\vert ^{2}$
), so this becomes
$\varepsilon \geq \sum_{j}\underbrace{\left\vert \gamma
_{j}\right\vert \left\vert \gamma _{j}^{\prime }\right\vert }_{=\left\vert
\overline{\gamma _{j}^{\prime }}\gamma _{j}\right\vert } =
\sum_{j}\left\vert \overline{\gamma _{j}^{\prime }}\gamma _{j}\right\vert
\geq \left\vert \sum_{j}\overline{\gamma _{j}^{\prime }}\gamma
_{j}\right\vert $
(by the triangle inequality). But since $1-\varepsilon
=\sum_{j}\overline{\gamma _{j}^{\prime }}\gamma _{j}$ (by comparing the $
\left( 1,4\right) $-th entry of the matrices in the above equation), this
becomes $\varepsilon \geq \left\vert 1-\varepsilon \right\vert $, what
contradicts the definition of $\varepsilon $.
Best Answer
Solving a linear system with that matrix is equivalent to solving the linear matrix equation $X + B^TXA + D^TXC = E$. As far as I know, it is an open problem how to exploit the structure of that matrix product to solve it (with a direct algorithm) in fewer than $O(n^5)$ operations (which is the complexity of $n^2$ steps of GMRES using the Kronecker structure to compute the products). It would be a big surprise to find a cheaper algorithm.
If one of the terms is small in norm or in rank, you can obtain a preconditioner for an iterative method by dropping it. That is currently your best bet to solve this problem. Pointers to some research in this direction: https://arxiv.org/abs/1704.02167, https://doi.org/10.1016/j.laa.2017.06.027.