Solving linear systems using Hermite Normal Forms

Python file for an implementation of Hermite normal form computations

This module offers an implementation of Hermite Normal Form computation. This algorithm can later be used to solve linear system within Euclidean domains.

This is an adaptation of method described in this paper for computing the Hermite Normal Form in Dedekind domains. Since Euclidean Domains are, in particular, Dedekind domains, we use the corresponding and adapted method.

The main differences with the algorithm described in that paper is that we do a row echelon form and that we do not force the diagonal to have 1. Moreover, we do not need to have maximal rank, obtaining zero rows in the end of the matrix.

Given a matrix \(M\) with \(n\) rows and \(m\) columns, a Hermite normal form (or HNF) is a matrix \(H\) equivalent to \(M\) (i.e., there is a unimodular matrix \(U\) such that \(UM = H\), also called transformation matrix) such that every element below the main diagonal is zero. This is similar to computing the echelon form as in a Gauss-Jordan elimination, but all operations stays in the same ring as the elements of \(M\) belong (given that it is an Euclidean domain).

AUTHORS:

• Antonio Jimenez-Pastor (2021-02-08): initial version

class ajpastor.misc.hermite.HermiteSolver(parent, matrix, inhomogeneous, euclidean=<function HermiteSolver.<lambda> at 0x7f192fa648b0>, xgcd=<function xgcd at 0x7f1940fed430>)

Bases: ajpastor.misc.linear_solver.LinearSystemSolver

This class represents the solving of a linear system using Hermite Normal Forms.

Solving a linear system using Hermite normal forms is possible by solving each of the equations and extending the result to the rest of the system. This can be done in every Euclidean Domain.

This algorithm is based on the Euclidean division algorithm and work anywhere where the elements of the parent structure have the methods __div__, __rem__ and xgcd.

INPUT:

• parent: the ring where the solutions will be searched. This parent can be a localization ring that can be provided as a triplet (R, g, d) where R is a euclidean domain, g is an empty list and d is a list of elements on R that will be localized.

• matrix: matrix for the system.

• inhomogeneous: the inhomogeneous vector for the system.

• euclidean: method for computing the euclidean division with remainder.

• xgcd: method for computing the Extended Euclidean GCD (by default, it takes the value of xgcd())

EXAMPLES:

sage: from ajpastor.misc.hermite import *

First we show that this solver gives the same Hermite normal forms as the usual method hermite_form:

sage: A = MatrixSpace(ZZ,2)([1,2,3,4])
sage: hs = HermiteSolver(ZZ, A, vector([0,0]))
sage: hs.echelon_form()
[ 1  0]
[ 0 -2]
sage: B = MatrixSpace(ZZ,5)(range(25))
sage: hs = HermiteSolver(ZZ, B, vector([0,0,0,0,0]))
sage: hs.echelon_form()
[  5   0  -5 -10 -15]
[  0   1   2   3   4]
[  0   0   0   0   0]
[  0   0   0   0   0]
[  0   0   0   0   0]
sage: C = matrix(ZZ,5,3,[1..15])
sage: hs = HermiteSolver(ZZ, C, vector([0,0,0,0,0]))
sage: hs.echelon_form()
[1 2 3]
[0 3 6]
[0 0 0]
[0 0 0]
[0 0 0]
sage: hs.transformation_matrix()
[ 1  0  0  0  0]
[ 4 -1  0  0  0]
[-1  2 -1  0  0]
[ 2 -3  0  1  0]
[ 3 -4  0  0  1]
sage: hs.U*hs.A == hs.H
True

Some special cases when we have 0 or 1 row:

sage: a = matrix(ZZ, 1,2,[0,-1])
sage: hs = HermiteSolver(ZZ, a, vector([0]))
sage: hs.echelon_form()
[ 0 -1]
sage: b = matrix(ZZ, 1, 3)
sage: hs = HermiteSolver(ZZ, b, vector([0]))
sage: hs.echelon_form()
[0 0 0]

Since this class works for any Euclidean domain, we can also run the tests from the hermite_form() method for matrices with univariate polynomials:

sage: M.<x> = GF(7)[]
sage: A = matrix(M, 2, 3, [x, 1, 2*x, x, 1+x, 2])
sage: hs = HermiteSolver(M, A, vector([0,0]))
sage: hs.echelon_form()
[      x       1     2*x]
[      0       x 5*x + 2]
sage: hs.transformation_matrix()
[1 0]
[6 1]
sage: hs.U*hs.A == hs.H
True
sage: B = matrix(M, 2, 3, [x, 1, 2*x, 2*x, 2, 4*x])
sage: hs = HermiteSolver(M, B, vector([0,0]))
sage: hs.echelon_form()
[  x   1 2*x]
[  0   0   0]
sage: hs.transformation_matrix()
[0 4]
[5 1]
sage: hs.U*hs.A == hs.H
True

But the main goal of this class is to solve linear systems:

sage: A = matrix(M, 2, 3, [x, 1, 2*x, x, 1+x, 2])
sage: v = vector([x,1+x])
sage: hs = HermiteSolver(M, A, v)
sage: hs.solution()
(6*x + 6, x, 4*x + 4)
sage: hs.syzygy()
[5*x^2 + 5*x + 2]
[    2*x^2 + 5*x]
[            x^2]
sage: hs.A*hs.syzygy()
[0]
[0]
sage: hs.A*hs.solution() == hs.b
True

This class also work with localized rings over Euclidean domains. The definition of the localization is given in the input of this documentation:

sage: R.<x> = QQ[]
sage: A = matrix([[2, 3*x, 1/x],[1/(x+1), 2*x + 1, 1/(x^2+x)]])
sage: v = vector([3,x])
sage: hs = HermiteSolver((R, [], [x, 1+x]), A, v)
sage: hs.solution()
(-x^2 - x + 3, 0, 2*x^3 + 2*x^2 - 3*x)
sage: hs.syzygy()
[         (-2*x^2 - 1)/(x + 1)]
[                   -1/(x + 1)]
[(4*x^3 + 3*x^2 + 2*x)/(x + 1)]
sage: hs.A*hs.syzygy()
[0]
[0]
sage: hs.A*hs.solution() == hs.b
True

Sometimes, the system has zero solutions. Then a NoSolutionError is raised:

sage: A = matrix([[2, 3*x, -3],[1/(x+1), 2*x + 1, 2*x-1/2]])
sage: v = vector([1,x])
sage: hs = HermiteSolver((R, [], [x, 1+x]), A, v)
sage: hs.solution()
Traceback (most recent call last):
...
NoSolutionError: There is no solution to equation ...

But we can still extract the Hermite Normal Form and the transformation matrix:

sage: hs.echelon_form()
[      1/(x + 1)         2*x + 1       2*x - 1/2]
[              0 4*x^2 + 3*x + 2 4*x^2 + 3*x + 2]
sage: hs.transformation_matrix()
[      0       1]
[     -1 2*x + 2]