Introduction
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.
Citation
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$.
Algorithms
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:
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]
(*True*)
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]
printMatrixFactoriation[mfp]
recoverPoly[mfp] == p
Implementation
ClearAll[makeMonomialMultiplicative];
(* 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]},
Which[
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]}];
];
ClearAll[mxfactor];
(* 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;
Which[
mCnt < 1, Return[{Null, Null}];
, mCnt == 1, Return[Flatten@{First@(List @@ f), Rest@(List @@ f)}];
, mCnt > 1,
Do[
mPair = makeMonomialMultiplicative[monomials[[i]]];
Which[
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}];
]
]
ClearAll[combinePolynomialMatrixFactors];
(* 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[]] :=
Module[{aDim},
aDim = Last /@ Dimensions /@ {AA, BB, CC, DD};
If[AnyTrue[aDim, # != aDim[[1]] &], Message[combinePolynomialMatrixFactors::inconsistentDims]; Abort[]];
Which[
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}
]
]
ClearAll[matrixFactorisePolynomial];
(* 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]];
Do[
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}];
];
ClearAll[recoverPoly];
(* 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]]];
ClearAll[printMatrixFactoriation];
printMatrixFactoriation[factorPair_] := Print[MatrixForm[factorPair[[1]]], ".", MatrixForm[factorPair[[2]]], " = ", MatrixForm@Simplify[factorPair[[1]].factorPair[[2]]]];
