Last active
December 21, 2015 21:09
-
-
Save patrickmclaren/6366274 to your computer and use it in GitHub Desktop.
Determine the conditions under which a polynomial matrix has rank less than or equal to n.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
############################################################################## | |
# 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)) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment