Writing a matrix as a sum of matrices

edit

To start to follow the suggestion at https://ask.sagemath.org/question/605..., consider this function:

def sum_list(length, total):
    """
    INPUT:

    - length -- how many terms
    - total -- desired total

    Return list of lists, each of which has ``length`` entries chosen
    from {0, 1, -1} and adds up to ``total``.
    """
    C = IntegerListsLex(length=length, n=length+total, min_part=0, max_part=2)
    return [[a-1 for a in L] for L in C]

(It's only two lines long; everything between the """ lines is a comment.)

Sample use to get all lists of length 4 that add up to 1:

sage: sum_list(4, 1)
[[1, 1, 0, -1],
 [1, 1, -1, 0],
 [1, 0, 1, -1],
 [1, 0, 0, 0],
 [1, 0, -1, 1],
 [1, -1, 1, 0],
 [1, -1, 0, 1],
 [0, 1, 1, -1],
 [0, 1, 0, 0],
 [0, 1, -1, 1],
 [0, 0, 1, 0],
 [0, 0, 0, 1],
 [0, -1, 1, 1],
 [-1, 1, 1, 0],
 [-1, 1, 0, 1],
 [-1, 0, 1, 1]]

Now try to build up matrices using sum_list(LENGTH, 1) and sum_list(LENGTH, -1) based on the entries of your given matrix.

The following is a draft, and it could probably be sped up, but it seems to work. It treats a solution (A,B) as different from (B,A): the order matters.

def sum_list(length, total):
    """
    INPUT:

    - length -- how many terms
    - total -- desired total

    Return list of lists, each of which has ``length`` entries chosen
    from {0, 1, -1} and adds up to ``total``.
    """
    C = IntegerListsLex(length=length, n=length+total, min_part=0, max_part=2)
    return [[a-1 for a in L] for L in C]

def mat_list(mat, length):
    """
    INPUT:

    - mat -- matrix, assumed to have entries in {1, -1}
    - length -- how many terms

    Return list of lists, each of which has ``length`` entries and
    adds up to ``mat``, and each entry is either a `(0, 1)` or a `(0,
    -1)` matrix.
    """
    sum_plus = sum_list(length, 1)
    sum_minus = sum_list(length, -1)
    old_attempts = []
    for x in mat.list():
        if x == 1:
            sums = sum_plus
        elif x == -1:
            sums = sum_minus
        else:
            raise ValueError('each entry of matrix should be 1 or -1')
        if not old_attempts: # initial step
            new_attempts = [[[s] for s in S] for S in sums]
        else:
            new_attempts = []
            for mats in old_attempts:
                NEW = []
                for S in sums:
                    new_M = [M + [s] for (M, s) in zip(mats, S)]
                    # keep only entries that correspond to (0,1) or (0,-1) matrices
                    if all(max(entries) - min(entries) < 2 for entries in new_M):
                        NEW.append([M + [s] for (M, s) in zip(mats, S)])
                new_attempts.extend(NEW)
        old_attempts = new_attempts

    matrices = [[matrix(mat.base_ring(), mat.nrows(), mat.ncols(), entries) for entries in L] for L in old_attempts]
    return matrices

With this code:

sage: mat = matrix(ZZ, [[1, 1, 1, 1], [1, -1, 1, -1], [1, 1, -1, -1], [1, -1, -1, 1]])
sage: %time L = mat_list(mat, 3)
CPU times: user 8.04 s, sys: 141 ms, total: 8.18 s
Wall time: 8.18 s
sage: len(L)
179328

EDIT: here is a modification, in case you don't want to distinguish the same lists with a different ordering:

def mat_list_NEW(mat, length):
    """
    INPUT:

    - mat -- matrix, assumed to have entries in {1, -1}
    - length -- how many terms

    Return list of lists, each of which has ``length`` entries and
    adds up to ``mat``, and each entry is either a `(0, 1)` or a `(0,
    -1)` matrix.
    """
    sum_plus = sum_list(length, 1)
    sum_minus = sum_list(length, -1)
    old_attempts = []
    for x in mat.list():
        if x == 1:
            sums = sum_plus
        elif x == -1:
            sums = sum_minus
        else:
            raise ValueError('each entry of matrix should be 1 or -1')
        if not old_attempts: # initial step
            new_attempts = set(tuple([(s,) for s in S]) for S in sums)
        else:
            new_attempts = set()
            for mats in old_attempts:
                NEW = set()
                for S in sums:
                    new_M = tuple(sorted(M + (s,) for (M, s) in zip(mats, S)))
                    # keep only entries that correspond to (0,1) or (0,-1) matrices
                    if all(max(entries) - min(entries) < 2 for entries in new_M):
                        NEW.add(new_M)
                new_attempts.update(NEW)
        old_attempts = new_attempts

    matrices = [[matrix(mat.base_ring(), mat.nrows(), mat.ncols(), entries) for entries in L] for L in old_attempts]
    return matrices

The differences: use the set structure so that duplicates are automatically removed, and sort the lists before adding them to the set, so that it's easy to recognize duplicates. Elements of a set have to be hashable so use tuples instead of lists. With this new version:

sage: mat = matrix(ZZ, [[1, 1, 1, 1], [1, -1, 1, -1], [1, 1, -1, -1], [1, -1, -1, 1]])
sage: %time L = mat_list_NEW(mat, 3)
CPU times: user 1.47 s, sys: 41 ms, total: 1.51 s
Wall time: 1.53 s
sage:  len(L)
29889

But the lengths of the lists are still going to get large as the number of summands increases...

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...

This draft code may also be seen(though there are repetitions)

def add(l1,l2):
       if len(l1)!=len(l2):
                 return('Error in lengths..!')
       new = []
       for ee in range(len(l1)):
                 #print(vector(l1[ee]))
                 new.append(list(vector(l1[ee])+vector(l2[ee])))
       return new

def mat_set(l1):
       temp = []
       for a in (l1):
                 temp.append(matrix(m , m, a))
       return temp

def add_n(a):
       test = []
       mat_A = [[0]*(len(M_list))]*n
       mat_A = matrix(mat_A)
       d = 0
       while d<len(M_list):
                mat_A += matrix(main[d][a[d]])
                d=d+1
       return list(mat_A)

def test_1_m1(com):
          for a in com:
                  cc=0
                  cd=0
                  for ai in a:
                          if ai==1:
                                cc=1
                          if ai==-1:
                                cd=1
                          if (cc==1 & cd ==1):
                                return 0
          return 1

m = int(input("Enter the order m of your matrix (m x m):"))
print("Enter the elements of your matrix (size: "+str(m)+"x"+str(m)+") row wise")
temp = []
for i in range(m^2):
        temp.append(int(input()))
M = matrix(m,m,temp)
show("M =",M)
temp = []
for i in range(M.nrows()-1):
        temp +=list(vector(M[i:]))
M_list=temp


n = int(input("In how many matrices you want to break M:"))
T=Tuples((0,-1,1),n)

sum_1  = []
sum_m1 = []
for a in T:
       if sum(a)==1:
               sum_1.append(a)
       if sum(a)==-1:
               sum_m1.append(a)

temp_comb = []
for i in range(n):
       temp_comb.append([0]*len(M_list))

main = []
for x in range(len(M_list)):
         M_elem = x
         main_temp = []
         if M_list[M_elem]==1:
                   ss = sum_1
         if M_list[M_elem]==-1:
                   ss = sum_m1
         for w in range(len(ss)):
                   temp_comb = []
                   for i in range(n):
                             temp_comb.append([0]*len(M_list))
                   for j in range(n):
                             temp_comb[j][M_elem] = ss[w][j]
                   main_temp.append(temp_comb)
         main.append(main_temp)

V = (Integers(len(sum_1)))^(len(M_list))
main_com = []
for a in V:    
        com = add_n(a)
        if vector([0]*(len(M_list))) not in com:
                 if test_1_m1(com)==1:
                           main_com.append(com)
                           show(mat_set(com))