Determine the conditions under which a polynomial matrix has rank less than or equal to n.

new_polynomial_matrix.py

##############################################################################
# System Libraries
##############################################################################
from collections import defaultdict
from sage.all import *

import itertools

##############################################################################
# Debugging
##############################################################################

import pdb

##############################################################################
# Helper Functions
##############################################################################

def _join_dicts_of_lists(*ddicts):
    new_dict = defaultdict(list)

    for ddict in ddicts:
        for k,v in ddict.items():
            new_dict[k].extend(v)

    return new_dict

def _clean_dict_of_lists(ddict):
    for k,v in ddict.items():
        new_v = []
        for el in v:
            if not el in new_v:
                new_v.append(el)

        ddict[k] = new_v

    return ddict

##############################################################################
# new_polynomial_matrix.py
##############################################################################

class MMatrix():
    def __init__(self, mat):
        self.matrix = matrix(copy(mat))

    def __add__(self, other):
        return MMatrix(self.matrix + other.matrix)

    def __sub__(self, other):
        return MMatrix(self.matrix - other.matrix)

    def __mul__(self, other):
        return MMatrix(self.matrix * other.matrix)

    def __getitem__(self, i):
        return copy(self.matrix[i])

    def __setitem__(self, i, new_row):
        try:
            self.matrix[i] = new_row
        except Exception:
            rows = []
            for idx, row in enumerate(self.matrix):
                if idx == i:
                    rows.append(new_row)
                else:
                    rows.append(list(row))

            self.matrix = matrix(rows)

    def __deepcopy__(self):
        return MMatrix(self.matrix)

    def row_size(self):
        return self.matrix.dimensions()[0]

    def column_size(self):
        return self.matrix.dimensions()[1]

    def row_reduce(self):
        """Returns a tuple, (reduced matrix, non zero entries)"""
        new_mat = self.__deepcopy__()
        non_zero = []

        for i in range(new_mat.max_pivots()):
            if new_mat[i][i] == 0:
                found_row = False
                for j in range(i+1, new_mat.row_size()):
                    if new_mat[i][j] != 0:
                        new_mat.swap_rows(i, j)
                        found_row = True
                        break
                if not found_row:
                    return (new_mat, non_zero)

            pivot = new_mat[i][i]
            non_zero.append(pivot)
            new_mat.scalar_mult(1/pivot, row=i)

            for j in range(i+1, new_mat.row_size()):
                new_mat.subtract_row(mat[i][j], i, j)

            new_mat.scalar_mult(pivot, row=i)

        return (new_mat, non_zero)

    def scalar_mult(self, scalar, row=None):
        if row is None:
            for i in range(self.row_size()):
                self[i] = [x * scalar for x in self[i]]
        else:
            self[row] = [x * scalar for x in self[row]]

    def swap_rows(self, i, j):
        row_i = self[i]
        row_j = self[j]

        self[i] = row_j
        self[j] = row_i

    def subtract_row(self, mult, i, j):
        row_i = [x * mult for x in self[i]]
        row_j = self[j]

        self[j] = [row_j[k] - row_i[k] for k in range(len(row_j))]

    def get_pivots(self):
        if not self._is_reduced():
            raise Exception

        return [self[i][i] for i in range(self.max_pivots())]

    def _is_reduced(self):
        for i in range(self.max_pivots()):
            for j in range(i+1, self.row_size()):
                if self[i][j] != 0:
                    return False

        return True

    def max_pivots(self):
        return min(self.row_size(), self.column_size())

    def apply_fn(self, fn, *args, **kwargs):
        new_mat = self.__deepcopy__()

        for i in range(new_mat.row_size()):
            for j in range(new_mat.column_size()):
                new_mat.update_entry(i, j, fn(mat[i][j], *args, **kwargs))

        return new_mat

    def update_entry(self, i, j, val):
        tmp_row = self[i]
        tmp_row[j] = val
        self[i] = tmp_row

def poly_div_rem(f, g):
    """Divide f by g, and return the remainder."""
    return f.quo_rem(g)[1]

def find_conditions_rank_leq_n(mat, n):
    conditions = defaultdict(list)

    for i in reversed(range(n + 1)):
        reduced_mat, non_zero = mat.row_reduce()

        for pivot_subset in itertools.combinations([k for k in range(reduced_mat.max_pivots())], i):
            # only consider unique, non-zero pivots
            non_zero_pivots = list(set([reduced_mat[x][x] for x in pivot_subset if reduced_mat[x][x] != 0]))

            conditions[len(non_zero_pivots)].append(non_zero_pivots)

        for el in non_zero:
            next_mat = reduced_mat.apply_fn(poly_div_rem, el)
            conditions = _join_dicts_of_lists(conditions, find_conditions_rank_leq_n(next_mat, i))

    conditions = _clean_dict_of_lists(conditions)

    return conditions

if __name__ == "__main__":
    R, (x, y) = PolynomialRing(QQ, 2, 'xy').objgens()

    f = x + y
    g = 0
    h = 0
    k = x + y

    mat = MMatrix([[f, g], [h, k]])
    n = 2

    print(find_conditions_rank_leq_n(mat, n))

polynomial_matrix.py

import itertools
import copy

from sage.all import *

import pdb

def row_reduce(mat, non_zero):
    """Let the entries belong to the field of rational functions,
    then we perform Gaussian elimination on the matrix, taking care to
    note that certain rows may be zero."""

    # Row size
    size = mat.dimensions()[0]

    # Iterate over pivots from upper-left to lower-right
    for i in range(size):
        # Check for row-exchange
        if mat[i][i] == 0:
            for k in range(i + 1, size):
                mat.swap_rows(i, k)

                # Found a row
                if mat[i][i] != 0:
                    break

            # No more nonzero rows
            if mat[i][i] == 0:
                continue

        # Ensure the pivot is of value 1
        multiple = mat[i][i].lc()
        for j in range(size):
            new_row = copy_row(mat, i)
            new_row[j] = mat[i][j] / multiple
            mat[i] = new_row

        # It is at this point that we must ensure that
        # the pivot we're dividing by is not equal to zero.
        # We store it, so we can keep track of it.
        non_zero.append(mat[i][i])

        # Clear the column, we're now moving top to bottom
        for j in range(size):
            # Don't clear the row we're working on though
            if i == j:
                continue

            quotient = mat[j][i] / mat[i][i]

            new_row = copy_row(mat, j)

            for n in range(size):
                new_row[n] = mat[j][n] - quotient*mat[i][n]

            mat[j] = new_row

    return

def copy_row(mat, i):
    """Return a copy of row i from matrix."""
    return list(mat[i])

def reduce_entry(f, mat):
    row_length = mat.dimensions()[0]
    col_length = mat.dimensions()[1]

    for i in range(row_length):
        for j in range(col_length):
            new_row = copy_row(mat, i)
            new_row[j] = mat[i][j].quo_rem(f)[1]
            mat[i] = new_row

def get_pivots(mat):
    n = mat.dimensions()[0]
    pivots = []
    for i in range(n):
        if mat[i][i] != 0:
            pivots.append(i)

    return pivots

def max_rank(mat):
    return len(get_pivots(mat))

def substitute_entry(f, g, mat):
    """Substitutes all occurrences of f for g in matrix."""
    row_length = mat.dimensions()[0]
    col_length = mat.dimensions()[1]

    for i in range(row_length):
        for j in range(col_length):
            if mat[i][j].quo_rem(f)[1] == 0:
                new_row = copy_row(mat, i)
                new_row[j] = g
                mat[i] = new_row

def ideal_for_rank_leq_n(n, mat, poly_ring):
    """Find conditions under which the matrix has rank less 
    than or equal to n. Conditions are given as an ideal."""

    row_length = mat.dimensions()[0]
    mrank = max_rank(mat)

    offset = 0
    for i in range(n + 1):
        if n - i > mrank:
             offset += 1
             continue

        for family in itertools.combinations(get_pivots(mat), i - offset):
            pivots = [mat[x][x] for x in family]
            pivots = poly_ring.ideal(pivots)

def find_all_conditions_for_rank_leq_n(n, mat, poly_ring):
    # We need an unreduced copy of the matrix
    cmat = copy(mat)

    # We must assume a few polynomials are non-zero in order to
    # row reduce.
    non_zero = []
    row_reduce(cmat, non_zero)
    ideal_for_rank_leq_n(n, cmat, poly_ring)

    for entry in non_zero:
        next_mat = copy(mat)
        reduce_entry(entry, next_mat)
        #substitute_entry(entry, 0, next_mat)
        find_all_conditions_for_rank_leq_n(n, next_mat, poly_ring.quotient([entry]))

def main():
    """Main function, everything comes from here."""

    R, (x, y) = FractionField(PolynomialRing(QQ, 2, 'xy')).objgens()
    
    f = 2*x**2 - 2*x
    g = 0
    h = 4*x**2 - 4*x
    k = 4*y**2 - 4*y

    pdb.set_trace()
    
    mat = matrix([[f, g], [h, k]])

    # R, (a, b, c) = FractionField(PolynomialRing(QQ, 3, 'abc')).objgens()
    # 
    # x_1_1 = a
    # x_1_2 = 0
    # x_1_3 = 0
    # x_2_1 = 0
    # x_2_2 = b
    # x_2_3 = 0
    # x_3_1 = 0
    # x_3_2 = 0
    # x_3_3 = c
    # 
    # mat = matrix([[x_1_1, x_1_2, x_1_3], [x_2_1, x_2_2, x_2_3], [x_3_1, x_3_2, x_3_3]])
    # 
    # pdb.set_trace()

    find_all_conditions_for_rank_leq_n(3, mat, R)
    
if __name__ == "__main__":
    main()