How to find an algebraic representation

I want to find matrix representation of some algebra (Clifford algebra Cl(3,1) in this case). Here is an example, which I would like to get and then extend to higher dimensional matrices.

Suppose we look for 4x4 matrix representation and want to find 4 matrices (matrix dimensions are not directly related to number of matrices)

em1 = Array[a1, {4, 4}]; em2 = Array[a2, {4, 4}];
em3 =  Array[a3, {4, 4}]; em4 = Array[a4, {4, 4}];


I know that the algebra representation can be found using SolveAlways as follows

SolveAlways[
{Det[c[1]*em1 + c[2]*em2 + c[3]*em3 +
c[4]*em4] == (c[1]^2 + c[2]^2 + c[3]^2 - c[4]^2)^2,
Det[em1] == 1, Tr[em1] == 0, Det[em2] == 1, Tr[em2] == 0,
Det[em3] == 1, Tr[em3] == 0, Det[em4] == 1, Tr[em4] == 0},
{c[1],c[2], c[3], c[4]}]


Unfortunatelly, this runs forever (>2 days on my laptop). Below is the test that at least one of representation, which satisfy the the above conditions is

em1test = {{0, 0, 1, 0},
{0, 0, 0, 1},
{1, 0, 0, 0},
{0, 1, 0, 0}};
em2test = {{0, 0, -I, 0},
{0, 0, 0, -I},
{I, 0, 0, 0},
{0, I, 0, 0}};
em3test = {{0, -I, 0, 0},
{I, 0, 0, 0},
{0, 0, 0, I},
{0, 0, -I, 0}};
em4test = {{I, 0, 0, 0},
{0, -I, 0, 0},
{0, 0, -I, 0},
{0, 0, 0, I}};


Indeed,

testRules =


and

Det[c[1]*em1 + c[2]*em2 + c[3]*em3 + c[4]*em4] /. testRules // Simplify
(\* (c[1]^2 + c[2]^2 + c[3]^2 - c[4]^2)^2 \*)


Unfortunatelly separating real and imaginary parts seems do not reduce complexity of the problem (same number of unknowns).

From documentation I know that SolveAlways can be rewriten in terms of, say Solve and Eliminate. So, I hope, that it can be rewriten in terms of GroebnerBasis itself. How can this be done? How far I can hope to go for looking for representations of higher algebras? Can I hope say, 8x8 dimensional representation for 7 matrices to be reachable within this approach? Any other ideas?

Edit 1

Answering to the comment I can reformulate the problem in terms of hermitian Transpose[Conjugate[m]]=m and antihermitian Transpose[Conjugate[m]]=-m matrices as follows:

(reProto={{p[1,1],p[1,2],p[1,3],p[1,4]},{p[1,2],p[2,2],p[2,3],p[2,4]},{p[1,3],p[2,3],p[3,3],p[3,4]},{p[1,4],p[2,4],p[3,4],p[4,4]}})//MatrixForm

(imProto={{0,i[1,2],i[1,3],i[1,4]},{-i[1,2],0,i[2,3],i[2,4]},{-i[1,3],-i[2,3],0,i[3,4]},{-i[1,4],-i[2,4],-i[3,4],0}})//MatrixForm


and then

em1=(reProto+I*imProto)/.{p->ar1,i->ai1};
em2=(reProto+I*imProto)/.{p->ar2,i->ai2};
em3=(reProto+I*imProto)/.{p->ar3,i->ai3};
em4=(I*reProto+imProto)/.{p->ar4,i->ai4};


Unfortunatelly, counting variables we see that number of unknowns did't change, i.e. 64+4, except that now solutions should be real numbers. Unfortunatelly I don't know how to use this restriction in SolveAlways.

Edit2

Answering the comment about restructions I reformulated the problem to include all known constraints

 SolveAlways[Flatten[{Det[c[1]*em1+c[2]*em2+c[3]*em3+c[4]*em4]==(c[1]^2+c[2]^2+c[3]^2-c[4]^2)^2,
Tr[em1]==0,Tr[em2]==0,Tr[em3]==0,Tr[em4]==0,

• I think this might be intractable as posed. You have 9 equations and 68 variables. Are there more relations to consider? For example, are the matrices intended to be unitary (as is the case for the xxxtest matrices above)? – Daniel Lichtblau Nov 10 '15 at 17:12
• To impose the Hermitian-by-signature-parity I think you really do need to double your variables into real and imaginary components. Anything that can be done to break the symmetries will also help. After all this is done, I expect the problem will be intractable. If you set it up as a GroebnerBasis invocation, I'd consider using settings MonomialOrder -> DegreeReverseLexicographic, CoefficientDomain -> \ RationalFunctions where the c[j] for j in Range[4] are in effect parameters (so the "variables" argument in GB would be Variables[{em1, em2, em3, em4}]). – Daniel Lichtblau Nov 11 '15 at 15:19