Is there a way to solve a linear system of equations where the coefficients do not commute? (Warning: my two-days of looking so far suggest this may be an unsolved math problem, but I figured I could start here.)
For a minimal-working example:
$$ \begin{array}{rcl} 0 & = & M_1 x_1 + M_2 x_2 + M_3 x_3 , \\ 0 & = & M_4 x_1 + M_5 x_2 + M_6 x_3 \end{array} $$ where $$ M_i M_j \ne M_j M_i. $$ Is there an algorithmic way to solve this system?
In this case, a solution by hand is straightforward: $x_3 = f(M, M^{-1}) x_1$, and $x_2 = g(M, M^{-1}) x_1$, where $f$ and $g$ are some products and sums of the $M$'s and their inverses.
Can Mathematica do this for more complicated systems?
Further comments:
I've looked a little bit at the Groebner basis features of NCAlgebra (also mentioned in this question), but that doesn't look like it can go the whole way.
It looks to me like mathematicians have only solved such expressions when the $M$'s are the same, maybe? See this paper.
I have in mind specifically the case where the $M's$ are not matrices but rather some arbitrary operators acting on the $x$'s. In particular, in the equations I'm dealing with the $M$'s represent some relabeling of the (suppressed) arguments of the $x$'s. Something like $M_1 x_1 \equiv M_1[x_1({1,2,3,4})] = x_1[{3,1,2,4}]$; so the equations I care about are actually functional equations relating functions with different arguments.
Edit:
Including MWE of using NCAlgebra:
(* Set up equations *)
(* I used 'o' here instead of 'M'*)
(* and 'f' instead of 'x' *)
eqns = {
o[1] ** f[1] + o[2] ** f[2] + o[3] ** f[3]
,
o[4] ** f[1] + o[5] ** f[2] + o[6] ** f[3]
};
(* Useful for defining variables *)
ClearAll[or, os, fs]
or = Range[1, 6];
os = Join[o /@ or, oi /@ or];
fs = f /@ Range[1, 3];
(* Defining the variables *)
vars = Join[os, fs];
(* For the NCAlgebra package *)
SetNonCommutative[vars];
SetMonomialOrder[vars];
(* Defining the operator inverse equations *)
invEqns = Flatten[{o[#] ** oi[#] - 1, oi[#] ** o[#] - 1} & /@ or];
(* Find Groebner basis *)
NCMakeGB[Join[eqns, invEqns], 1]
(*
{o[1] ** f[1] + o[2] ** f[2] + o[3] ** f[3],
o[4] ** f[1] + o[5] ** f[2] + o[6] ** f[3], -1 + o[1] ** oi[1], -1 +
oi[1] ** o[1], -1 + o[2] ** oi[2], -1 + oi[2] ** o[2], -1 +
o[3] ** oi[3], -1 + oi[3] ** o[3], -1 + o[4] ** oi[4], -1 +
oi[4] ** o[4], -1 + o[5] ** oi[5], -1 + oi[5] ** o[5], -1 +
o[6] ** oi[6], -1 + oi[6] ** o[6],
f[3] + oi[3] ** o[1] ** f[1] + oi[3] ** o[2] ** f[2],
f[3] + oi[6] ** o[4] ** f[1] + oi[6] ** o[5] ** f[2]}
*)
So the last line of the output gets f[3]
in the form I want, but it's not clear how to proceed from there.