According to the work of David Eisenbud cited below, it is possible to factorize almost any polynomial with matrices (exclusions ~ no linear terms), and this has apparently been of great use to string theorists, amongst others, and I understood him to say (Numberphile YouTube Video) that the result is not merely an existence theorem, there is an algorithm for doing it.

(The paper is a proof of a generalisation of the technique invented by Paul Dirac that resulted in Dirac Matrices)

A cursory examination of the original paper suggests that even though the result may be simply stated, the algorithm is not necessarily simple (I have not even found a description of it) and it would probably require both advanced mathematical and advanced Mathematica knowledge.

Does anyone know of an existing implementation of the algorithm in Mathematica?

A quick check of the current online documentation for Factor does not indicate that matrix factorisation is an option.

Eisenbud, David. ‘Homological Algebra on a Complete Intersection, with an Application to Group Representations’. Transactions of the American Mathematical Society 260, no. 1 (1980): 35–64.

  • 1
    $\begingroup$ Have you already seen this? $\endgroup$ May 17, 2020 at 8:34
  • $\begingroup$ @J.M. - No, I had not seen that; thank you. The inefficiency of the "brute force" approach leading to $2n \times 2n$ matrices is perhaps to be expected; apart from the special case considered, I wonder whether there are any other approaches... will try playing with the simple approach first... thanks again. $\endgroup$ May 17, 2020 at 13:33
  • $\begingroup$ @J.M. - thanks to your direction to that paper I now have an implementation and have improved upon the algorithm to stop exponential growth of the factor matrices, i.e. you decide how big you want the matrices to be and it works to that fixed size. Description, code, examples etc. is only a few pages (it would still be long answer though), alas M2MD for Mathematica to Markdown by kuba doesn't quite do enough... any ideas how I might easily format it to include as an answer here? NB I don't use github and have no plans to :) $\endgroup$ May 23, 2020 at 14:32
  • $\begingroup$ @kuba for comment above, just in case you have other suggestions... M2MD nice but left a lot of work to do... $\endgroup$ May 23, 2020 at 14:33
  • $\begingroup$ I'm fussy about formatting myself, so I'm the wrong one to ask about "easy" formatting. If it will take time to write something presentable, then by all means, take the time you need. $\endgroup$ May 24, 2020 at 8:04

1 Answer 1



Thanks to the paper identified by @J.M., an implementation of the standard method is presented, below - together with (what I think may be) a novel refinement that allows matrix growth to be avoided.

Code was viable before transformation to markdown; apologies for any errors that may have crept in as a result.


Crisler, David, and Diveris, Kosmas. 'Matrix Factorizations of Sums of Squares Polynomials', 21 October 2016, 8

Definitions & Notation

We write matrix multiplication as $A\cdot B$.

An $n \times n$ matrix factorisation of a polynomial $f \in S$ , where $S$ is the ring $\mathbb{R}\left[x_1, x_2, \ldots , x_m\right]$, is a pair of $n \times n$ matrices $A$ and $B$ such that $A\cdot B$ = $f I_n$, where $I_n$ is the $n \times n$ identify matrix, i.e. each non-zero element of $A\cdot B$ is a copy of $f$.


In the paper cited, the authors show how to inductively construct a matrix factorisation of polynomials of the form $f_k = g_1 h_1 + g_2 h_2 + ... g_k h_k$ using the technique of Knörrer. The disadvantage of the general method is that for a polynomial of $n$ terms, resultant matrices are $2^{n - 1} \times 2^{n - 1}$, i.e. matrix size increases exponentially with the number of terms (although a method for generating smaller factorisations of "sums of squares polynomials" $f_n= g_1^2 + g_2^2 + ... g_n^2$ for $4 \leq n \leq 8$ was given, it is not of more general applicability).

A refinement is presented here - in addition to a simple Mathematica implementation of the standard method - that allows polynomials to be processed in groups of $m$ terms, where $m \lt n$, so that the resulting matrices are of constant size $2^{m - 1} \times 2^{m - 1}$. In dealing with Mathematica representations of the polynomials to be factorised we note that:

  • Polynomials nominally comprise monomials of integral powers - however the code below does permit real value exponents; "monomial" and "(additive) term" are therefore considered synonymous

  • Mathematica expressions may need to be coerced to specific forms in order to handle $a \times a$ effectively as Times[a, a], whose separable parts are of the same kind, and not as Power[a,2], whose parts are not of the same kind and therefore not separable as required.

  • If complex numbers are used, it is not always possible to recover the original form of the polynomial where I * I terms occur and become -1, however Simplification should demonstrate equality where it exists

Code has been written for clarity rather than elegance or efficiency, and only limited error handling is provided.

There are various ways in which the code could be developed further:

  • User control over the splitting of powers in makeMonomialMultiplicative

  • User control over the extraction of parts of monomials in makeMonomialMultiplicative (currently just First/Rest)

To reproduce the inline examples, the functions provided later must first be defined; in a mathematica notebook the code would typically be in Initialization cells.

The refinement: Pairwise Factorisation and Summation

Consider a polynomial $p$ of four terms expressed a sum of two, two-term polynomials $p1$ and $p2$

p1 = x1 y1 z1 + x2 y2 z2 ;
p2 = x3 y3 z3 + x4 y4 z4;
p = p1 + p2;
AA = mxfactor[p1][[1]]; BB = mxfactor[p1][[2]]; 
CC = mxfactor[p2][[1]]; DD = mxfactor[p2][[2]];

The pairs {AA, BB} and {CC, DD} are the matrix factorisations of p1 and p2 respectively, hence,

$$p I_2 = AA.BB + CC.DD$$

However, here we have the sum of a pair of factorisations where a single pair is required so that we can iterate over a long polynomial, processing sub-polynomials individually and accumulating the results into matrices of a constant size.

Fortunately, we can obtain $p I_2 = EE\cdot FF$ simply: with the aid of elementary matrix operations, we eliminate one pair of terms by pre- and post-multiplying by the inverses of elements in a pair (e.g. CC, DD) and absorbing the identify matrix by representing it as e.g. Inverse[DD].DD and using the distributivity of matrix multiplication.

In fact, by simple permutation we can obtain $EE\cdot FF = AA\cdot BB + CC\cdot DD$ in four distinct ways, and allow the user to choose among them by an Option:

Simplify[(AA.BB + CC.DD).Inverse[BB].BB == ((AA.BB + CC.DD).Inverse[BB]).BB  == (AA + CC.DD.Inverse[BB]).BB] == 
Simplify[(BB.AA + CC.DD).Inverse[AA].AA == ((BB.AA + CC.DD).Inverse[AA]).AA  == (BB + CC.DD.Inverse[AA]).AA] == 
Simplify[(AA.BB + CC.DD).Inverse[DD].DD == ((AA.BB + CC.DD).Inverse[DD]).DD  == (CC + AA.BB.Inverse[DD]).DD] == 
Simplify[(AA.BB + DD.CC).Inverse[CC].CC == ((AA.BB + DD.CC).Inverse[CC]).CC  == (DD + AA.BB.Inverse[CC]).CC]


Tested forms

The following polynomials were factorised and recovered from the factorisation (with the caveat for complex numbers noted above).

p = x1   + x2 y2  + x3 y3 z3 + x4 y4 z4 a4;
p = Exp[Sin[y1]] + Sin[x2] Exp[y2 z2] + x3 Sin[y3] z3 + Sin[x4 y4 z4] + x5 y5;
p = g1^3.6 + h1 ^Pi   + i1 i2 + j1 j2;
p = (x1 + y1 I) (x2 + y2 I) + Sin[x3 + y3 I] + z^3

(* evaluate and compare, varying TermsPerFactorisation, Method as desired using *)
mfp = matrixFactorisePolynomial[p, "TermsPerFactorisation" -> 2, "Method" -> 3]
recoverPoly[mfp] == p



(* makeMonomialMultiplicative coerces form into two multiplied terms 
   so that the terms are suitable for use by mxfactor *)
makeMonomialMultiplicative[monomial_] := 
  (* The best way to produce arbitrary "monomials" from non-standard 
  "polynomials" that may have non-integral exponents is to parse 
  the main expression as a List and take the parts as "monomials" *)
    Module[{\[Alpha], \[Beta]},
              ToString@Head@monomial == "Times", \[Alpha] = First[monomial]; \[Beta] = Rest[monomial];
            , ToString@Head@monomial == "Power", 
                If[IntegerQ[monomial[[2]]] (* this works well if monomial was obtained as a part of a List *)
                    , \[Alpha] = Power[monomial[[1]], Floor[monomial[[2]]/2]]; \[Beta] = Power[monomial[[1]], Ceiling[monomial[[2]]/2]]; (* split powers as close to evenly as possible in integers*)
                    , \[Alpha] = Power[monomial[[1]], monomial[[2]]/2]; \[Beta] = Power[monomial[[1]], monomial[[2]]/2];(* if powers non-integral, then divide by 2 *)
            , True, \[Alpha] = 1; \[Beta] = monomial;
        Return[{\[Alpha], \[Beta]}];


(* mxfactor performs a matrix factorisation Based on Corollary 7 of 
   Crisler & Diveris; this will produce matrices that grow in size 
   exponentially with the number of terms in f *)
mxfactor[f_] := 
    Module[{A, B, Anew, Bnew, monomials = List @@ f, monomial, mCnt, 
   iMd = 1, unity, mPair},
        mCnt = Length@monomials;
              mCnt < 1, Return[{Null, Null}];
            , mCnt == 1, Return[Flatten@{First@(List @@ f), Rest@(List @@ f)}];
            , mCnt > 1,             
                        mPair = makeMonomialMultiplicative[monomials[[i]]];
                              i == 1, A = mPair[[1]]; B = mPair[[2]];
                            , True,
                                    If[SquareMatrixQ@A, iMd = IdentityMatrix@Last@Dimensions@A];
                                    Anew = ArrayFlatten[{{A, -mPair[[2]] iMd}, {mPair[[1]] iMd, B}}]; (* Convert block matrix to flat matrix *) 
                                    Bnew = ArrayFlatten[{{B, mPair[[2]] iMd}, {-mPair[[1]] iMd, A}}];
                                    A = Anew; B = Bnew;
                        , {i, 1, mCnt}
                    Return[{A, B}];


(* Let p1 = term1 + term2 and p2 = term3 + term 4 be two polynomials 
   (of possibly non-integral coefficients) and let {AA, BB}, {CC, DD}
   be their respective matrix factorisations, then
   combinePolymomialMatrixFactors[AA, BB, CC, DD] returns a new pair 
   of matrices, say {EE, FF}, such that {EE, FF} is a matrix factorisation 
   of p = p1 + p2.
   This method of combining solutions allows a polynomial in an 
   arbitrary number of terms n to be expressed as a matrix factorisation
   in terms a pair of 2^(m-1)\[Times]2^(m-1) matrices, where m is the 
   number of terms processed at once, rather than a pair of 
   2^(Length[p]-1)\[Times]2^(Length[p]-1) matrices *)

combinePolynomialMatrixFactors::inconsistentDims = "Error: the matrices are not all the same size."; (* The matrices must be 2D square, and commute pairwise, i.e. such that AA.BB = BB.AA, CC.DD = DD.CC *)
combinePolymomialMatrixFactors::invalidMethod    = "The option \"Method\" value must be in {1, 2, 3, 4}";
Options[combinePolymomialMatrixFactors]          = {"Method" -> 1}; 

combinePolymomialMatrixFactors[AA_?SquareMatrixQ, BB_?SquareMatrixQ, CC_?SquareMatrixQ, DD_?SquareMatrixQ, OptionsPattern[]] :=
        aDim = Last /@ Dimensions /@ {AA, BB, CC, DD};
        If[AnyTrue[aDim, # != aDim[[1]] &], Message[combinePolynomialMatrixFactors::inconsistentDims]; Abort[]];
              OptionValue["Method"] == 1, {AA + CC.DD.Inverse[BB], BB}
            , OptionValue["Method"] == 2, {BB + CC.DD.Inverse[AA], AA}
            , OptionValue["Method"] == 3, {CC + AA.BB.Inverse[DD], DD}
            , OptionValue["Method"] == 4, {DD + AA.BB.Inverse[CC], CC}
            , True, Message[combinePolymomialMatrixFactors::invalidMethod]; {Null, Null}


(* matrixFactorisePolynomial factorises a polynomial of abritrary \
   length into fixed matrices according to the number of terms to be \
   processed at once *)

matrixFactorisePolynomial::invalidOptionValue = "Invalid option value supplied; TermsPerFactorisation must be >= 2";
matrixFactorisePolynomial::invalidExprHead = "The expression to be factorised must be a sum of terms, i.e. with Head = Plus; the head was `1`.";
(* Option AutoExpand converts 2(a+b) into 2a + 2b, A simple fix for
   some common forms that can easily be made compliant with the need for a sum of terms; 
   Option Method is defined for and passed to combinePolymomialMatrixFactors *)
Options[matrixFactorisePolynomial] = {"AutoExpand" -> True, "TermsPerFactorisation" -> 2, "Method" -> 1}; 

matrixFactorisePolynomial[p_, OptionsPattern[]] :=
    Module[{poly = p, subpolys, zero, 
    tpf = OptionValue["TermsPerFactorisation"], mxf, cmpmf, AA, BB},
        If[OptionValue["AutoExpand"], poly = Expand[poly]]; (* A simple fix for some common forms that can easily be made compliant *)
        If[ToString@Head@poly != "Plus", Message[matrixFactorisePolynomial::invalidExprHead, Head@poly]; Abort[]];
        If[tpf < 2, Message[matrixFactorisePolynomial::invalidOptionValue]; Abort[]];
        subpolys = Partition[List @@ poly, tpf, tpf, 1, zero];
        mxf = mxfactor[subpolys[[1]]];
        AA = mxf[[1]]; BB = mxf[[2]];
            mxf = mxfactor[subpolys[[i]]];
            cmpmf = combinePolymomialMatrixFactors[AA, BB, mxf[[1]], mxf[[2]], Method -> OptionValue["Method"]];
            AA = cmpmf[[1]]; BB = cmpmf[[2]];
            , {i, 2, Length[subpolys]}
        Return[{AA /. zero -> 0, BB /. zero -> 0}];


(* recoverPoly recovers the polynomial from the product of the
   factorisation matrices, taking into account all entries for
   verification purposes, i.e. it sums the elements etc. rather 
   than just extract a single diagonal element of what should 
   be a multiple of the identity matrix *)

recoverPoly[factorPair_] := Total@Flatten@Simplify[factorPair[[1]].factorPair[[2]]/Last@Dimensions@factorPair[[1]]];

printMatrixFactoriation[factorPair_] := Print[MatrixForm[factorPair[[1]]], ".", MatrixForm[factorPair[[2]]], " = ", MatrixForm@Simplify[factorPair[[1]].factorPair[[2]]]];


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.