Solving linear system with non-commuting coeffcients

Question

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.

Answer 1

The reason why the Groebner basis calculation does not produce the answer you are seeking is because you need an ever growing set of assumptions on what can be inverted and not only the invertibility of the o's in order to move forward.

If you are fine with assuming all letters are invertible than you could try a symbolic matrix factorization. Using NCAlgebra version 5.0.0 you can set up the problem as follows:

M = {{o[1], o[2], o[3]}, {o[4], o[5], o[6]}} /. f_[i_] -> Subscript[f, i]
x = {f[1], f[2], f[3]} /. f_[i_] -> Subscript[f, i]

The rule in the end is to work with Subscripts, which automatically treat all lowercase o's and f's as noncommutative. The problem you're trying to solve is the linear algebra problem $M x = 0$, a problem that can be solved by factoring $M$. A factorization suitable for symbolic computation is the $LU$ decomposition with complete pivoting:

{lu, p, q, rank} = NCLUDecompositionWithCompletePivoting[M];
{l, u} = GetLUMatrices[lu];

from which a solution can be calculated as:

sol = Thread[x[[q[[1 ;; rank]]]] -> -NCDot[NCInverse[u[[All, 1 ;; rank]]], u[[All, rank + 1 ;;]], x[[q[[rank + 1 ;;]]]]]]

The solution obtained this way requires the assumption that $o_5 - o_4 o_1^{-1} o_2$ is invertible, which was never part of the input during the Groebner basis calculation.