Revision history [back]

Well... This problem (decomposing a matrix as the sum of $p$ matrices with your constraints) can be turned in the search of the ideal of some field defined by a polynomial system :

• each element of M will be translated as a sum of $p+1$ terms (the corresponding elements of the $p$ matrices minus the corresponding element of M) ;

• each $(0, 1)$ $r\times c$ matrix (generic term $m$) can be created by nullifying $rc$ $m(m-1)$ polynomials ;

• each $(0, -1)$ $r\times c$ matrix (generic term $m$) can be created by multiplying a $(0, 1)$ matrix by -1 ;

• therefore, all the matrices can be given the same form, by multiplying each of them by $s_k,\ k\in1\dots p$ with the additional constraints of nullifying $s_k^2-1$.

This can be implemented as follows :

def solve_system_for(M, Nmat=6):
    """
    Decompose M as a sum of Nmat matrices of the same shape, each of them
    being composed either of 1 and 0 OR -1 and 0
    """
    Nr = M.nrows()
    Nc = M.ncols()
    Nel = Nr*Nc
    # Names for signs of matrices
    sn = ["s%d"%u for u in (1..Nmat)]
    # Names for matrix elements
    mn = [[["m%d%d%d"%(t,u,v)for v in (1..Nc)] for u in (1..Nr)]
          for t in (1..Nmat)]
    # Polynomial ring : must be a field containing -1, 0, 1,
    # with arithmetic compatible with ZZ's ==> QQ
    R1 = PolynomialRing(QQ, sn+flatten(mn))
    R1.inject_variables()
    # Vector of inteterminates related to matrice's signs
    s = R1.gens()[0:Nmat]
    # vector of indeterminates for matrices elements
    a = R1.gens()[Nmat:]
    # List of matrices of indeterminates
    A = [s[u-1]*matrix(R1.gens()[Nmat+u*Nel:Nmat+(u+1)*Nel],
                       nrows=Nr, ncols=Nc)
         for u in range(Nmat)]
    # Matrix of equations to solve
    Eq = sum(A)-M.change_ring(QQ)
    # Constraints on s elements : s ∈ {-1, 1} <==> s in roots of s^2-1 in QQ
    Cs = [u^2-1 for u in s]
    # Constraints on a elements : a ∈ {O, 1} <==> a in roots of a*(a-1) in QQ
    Ca = [u*(u-1) for u in a]
    # Solutions ideal described by :
    J1 = R1.ideal(flatten(Eq.list()) + Cs + Ca)
    # The hard, long part :
    from time import time as stime
    t0 = stime()
    D = J1.dimension()
    t1 = stime()
    if D==0: # Expected case on the basis of finite set of possible values
        t2 = stime()
        S = J1.variety()
        t3 = stime()
        return {"Time for dimension":t1-t0,
                "Time for solution":t3-t2,
                "Dimension":D,
                "Solution(s)":S}
    if D<0: # No solution.
        return {"Time for dimension":t1-t0,
                "Time for solution":t3-t2,
                "Dimension":D,
                "Solution(s)":[]}
    if D>0:
        # Infinite solutions. This shouldn't happen,
        # but I learned the virtues of defensive programming...
        return {"Time for dimension":t1-t0,
                "Time for solution":t3-t2,
                "Dimension":D,
                "Parent ring":R1,
                "Ideal to reprocess":J1}

But I have some serious doubts regarding the usefulness of this implementation :

• Decomposing your $M$ as the sum of two matrices involves 34 indeterminates and 50 polynomials ; determining the variety of its ideal takes 4.17 seconds (0,3 seconds to check its dimension, 4 s for finding the two elements of its variety).

• Decomposing it in the sum of 3 matrices involves 51 indeterminates and 67 polynomials, and hasn't finished in about 11 hours.

Combinatorial explosion is fast...

Unless my reasoning or my programming is|are faulty (which is of course possible), I do not realistically contemplate the use of this solution to decompose M in a sum of 6 matrices...

Hoping this non-help helps finding a better solution...