# Solve won't find solution to simultaneous matrix equations

I'm trying to solve a pair of simultaneous matrix equations of the form,

$$AX + XA^\dagger + BJB^\dagger = 0,\ \ X C^\dagger + B J D^\dagger = 0,$$ where, $$J = \begin{bmatrix}1&0\\0&-1\end{bmatrix}$$

However sometimes Mathematica cannot find a solution:

\$Assumptions = {s0 \[Element] Reals, s0 > 0};
Module[{a, b, c, d,
n = 2, no = 2, X, ji = {{1, 0}, {0, -1}}},
{a, b, c,
d} = {{{(-2 + s0)/(2 + s0), Sqrt[s0 (4 + s0^2)]/(2 + s0)}, {Sqrt[
s0 (4 + s0^2)]/(2 + s0), -((2 + 3 s0)/(2 + s0))}}, {{0, 1/Sqrt[
2 + s0]}, {Sqrt[(2 + s0)/(
4 + s0^2)], -(2/Sqrt[(((2 + s0) (4 + s0^2))/s0)])}}, {{-4 Sqrt[
s0/(2 + s0)], -((2 Sqrt[8 + 4 s0 + 2 s0^2 + s0^3])/(
2 + s0))}, {-2 Sqrt[2 + s0], 0}}, {{1, 0}, {0, 1}}};
X = Array[x, {n, n}];

Print[Solve[
Simplify[
a.X + X.a\[ConjugateTranspose] + b.ji.b\[ConjugateTranspose]] ==
ConstantArray[0, {n, n}], Flatten[X]] // Simplify];
Print[Solve[
Simplify[X.c\[ConjugateTranspose] + b.ji.d\[ConjugateTranspose]] ==
ConstantArray[0, {n, no}], Flatten[X]] // Simplify];

(* Checking that solution solves both equations *)
Print[Simplify[(X /.
Solve[Simplify[
a.X + X.a\[ConjugateTranspose] +
b.ji.b\[ConjugateTranspose]] == ConstantArray[0, {n, n}],
Flatten[X]]) == (X /.
Solve[Simplify[
X.c\[ConjugateTranspose] + b.ji.d\[ConjugateTranspose]] ==
ConstantArray[0, {n, no}], Flatten[X]])]];

sols = Solve[
Simplify[
a.X + X.a\[ConjugateTranspose] + b.ji.b\[ConjugateTranspose]] ==
ConstantArray[0, {n, n}]
&& Simplify[
X.c\[ConjugateTranspose] + b.ji.d\[ConjugateTranspose]] ==
ConstantArray[0, {n, no}], Flatten[X]]]


The output is:

{{x[1,1]->-(1/(4+2 s0)),x[1,2]->Sqrt[s0/(4+s0^2)]/(2+s0),x[2,1]->s0/((2+s0) Sqrt[s0 (4+s0^2)]),x[2,2]->1/(4+2 s0)}}
{{x[1,1]->-(1/(4+2 s0)),x[1,2]->Sqrt[s0/(4+s0^2)]/(2+s0),x[2,1]->1/Sqrt[((2+s0)^2 (4+s0^2))/s0],x[2,2]->1/(4+2 s0)}}
True
{}


Note that the returned results for the first two Print statements (solving them separately) solve both equations, so why isn't Solve working?

• If you use Reduce instead of Solve you get this: The answer found by Reduce contains unsolved equation(s) ... A likely reason for this is that the solution set depends on branch cuts of Wolfram Language functions. Jul 1 '20 at 15:33
• There is a generic solution: $$X= (AB J D^\dagger-BJB^\dagger C^\dagger)(A^\dagger C^\dagger)^{-1}$$. Jul 1 '20 at 19:28

With Method->Reduce you get the result as ConditionalExpressions.

Solve[Simplify[
a.X + X.a\[ConjugateTranspose] + b.ji.b\[ConjugateTranspose]] ==
ConstantArray[0, {n, n}] &&
Simplify[X.c\[ConjugateTranspose] + b.ji.d\[ConjugateTranspose]] ==
ConstantArray[0, {n, no}], Flatten[X], Method -> Reduce]


Alternatively, you can reformulate the problem as a single equation

$$AXC^\dagger + XA^\dagger C^\dagger + BJB^\dagger C^\dagger = 0,\ \ AX C^\dagger + AB J D^\dagger = 0$$

$$XA^\dagger C^\dagger + BJB^\dagger C^\dagger -AB J D^\dagger= 0$$

FullSimplify[
Solve[X.a\[ConjugateTranspose].c\[ConjugateTranspose] +
b.ji.b\[ConjugateTranspose].c\[ConjugateTranspose] -
a.b.ji.d\[ConjugateTranspose] == ConstantArray[0, {n, n}],
Flatten[X]]]


• Your reformulation as a single equation is the way I will go! Not sure why I didn't think of that. Thank you! I still wonder why this is an issue in the first place though Jul 1 '20 at 18:45
• According to the documentation: "Solve gives generic solutions only. Solutions that are valid only when continuous parameters satisfy equations are removed." I did not look deep to work out why the conditions are needed. Jul 1 '20 at 18:53
• But why do you use Solve for a linear equation? Jul 1 '20 at 19:25