Endelig-genererte abelske grupper og Hermite-normalformen til heltallsmatriser

Dette blir en litt vanskeligere post enn ellers. Jeg vil kort skrive om et lite Python-skript jeg skrev for noen uker siden, og forklare hva det gjør.

Noen av leserne av denne bloggen (jeg vet ikke om dette er en ikke-tom mengde) har kanskje hatt et kurs i abstrakt algebra på universitetet, og lært om abelske grupper. Dette er matematiske objekter hvor du kan legge sammen og trekke fra, og slik at rekkefølgen ikke har noe å si (med andre ord er a+b = b+a alltid). Husker man spesielt godt, husker man også “fundamentalteoremet om endelig-genererte abelske grupper”. Dette teoremet sier at om du ikke trenger uendelig mange elementer for å beskrive gruppen din, kan den alltid beskrives ved hjelp av en n \times mm-matrise med heltallselementer.

Vi tar et par eksempler:

  • La G være heltallene modulo $11$. Dette betyr at 11 \equiv 0 i denne gruppen. Denne gruppen kan beskrives ved 1 \times 1-matrisen [11].
  • La G være gruppen av par (a,b), hvor a regnes modulo 11, og b er et vanlig heltall. Matrisen som presenterer denne gruppen er gitt ved \begin{pmatrix}11 & 0 \\ 0 & 0\end{pmatrix}.

Generelt vil gruppen \mathbb Z_{r_1} \oplus \mathbb Z_{r_2} \cdots \oplus \mathbb Z_{r_n} \oplus \mathbb Z^k beskrives ved en diagonalmatrise D med r_1,\ldots,r_n,0,0,\ldots,0 langs diagonalen. Dette er vel og fint, men presentasjonen er ikke unik! La meg forklare litt.

Denne presentasjon stammer fra å tenke på en heltallsmatrise som en avbildning \mathbb Z^n \xrightarrow{A} \mathbb Z^m. Den abelske gruppen \mathbb Z^m/im(A) (kokjernen til A) er da gruppen presentert av A. Ethvert basisskifte i \mathbb Z^n og \mathbb Z^m vil gi en isomorf kokjerne, så derfor er ikke A unik. Basisskifter svarer til å multiplisere A på venstre- eller høyresiden med inverterbare heltallsmatriser med determinant \pm 1. Dette svarer igjen til å utføre rad- eller søyleoperasjoner på A.

Så det vi ønsker, er å transformere A på en slik måte at det er lett å lese hvilken endelig-generert abelsk gruppe kokjernen er isomorf med. Helt ideelt ville det vært å få A på diagonalform, men det viser seg at vi ikke trenger å være så kravstore. Det holder at A er på triangulær-form, nemlig at alle elementer under diagonalen er null. Har vi A på denne formen, kan vi lese av diagonalelementene for å finne kokjernen.

Det finnes en algoritme for å gjøre dette, og formen matrisen ender opp på kalles for Hermite-normalform (strengt tatt er Hermite-normalformen en diagonalmatrise, men det trenger vi ikke). Det er en morsom øvelse å implementere denne i Python.

Det første jeg måtte gjøre, var å implementere min egen matrise-klasse i Python. Pakken Numpy har egne matriseobjekter, men disse konverterer alle elementene til flyttall, noe jeg ikke ville gjøre. Dette er ikke spesielt vanskelig, og også en morsom øvelse. Under gjengir jeg den relevante koden (det er mange flere ting man kan gjøre med matriser som jeg ennå ikke har implementert):

import ntheory
from operator import mul

class Matrix:
	"""
	Constructs a matrix object from a double list.
	"""
	def __init__(self, L):  
		self.L = L
		self.m = len(L[0]) # number of columns
		self.n = len(L) # number of rows

	def __add__(self,M):
		'''
		Returns the sum of self and M.
		'''
		newL = []
		for r in range(self.m):
			row = [self.L[r][i] + M.L[r][i] for i in range(self.n)]
			newL += [row]
		return Matrix(newL)

	def __sub__(self,N):
		return (self + (-N))

	def __neg__(self):
		newL = [[-r for r in R] for R in self.L]
		return Matrix(newL)

	def __mul__(self,N):
		'''
		Returns the product of self and N.
		'''
		NT = N.transpose()
		newL = []
		for i in range(self.n):
			newR = []
			for j in range(N.m):
				newR += [sum([self.L[i][k]*NT.L[j][k] for k in range(self.m)])]
			newL += [newR]
		return Matrix(newL)

	def transpose(self):
		'''
		Returns the transpose of self.
		'''
		newL = [[] for i in range(self.m)]
		for i in range(len(self.L)):
			for j in range(self.m):
				newL[j] += [self.L[i][j]]
		return Matrix(newL)

	def trace(self):
		'''
		If self is square, return trace.
		'''
		if self.n != self.m:
			return "NOT SQUARE"
		return sum([self.L[i][i] for i in range(self.n)])

	def switchRows(self,i,j):
		'''
		Returns the matrix obtained by switching rows i,j in self.
		Counting starts at 0.
		'''
		newL = list(self.L)
		newL[i], newL[j] = newL[j], newL[i]
		return Matrix(newL)

	def _prodDiagonal(self):
		return reduce(mul,[self.L[i][i] for i in range(len(self.L))],1)

	def det(self):
		return triangular(self)._prodDiagonal()


	def __str__(self):
		s = "{0}x{1}-matrix: ".format(self.m,self.n) + str(self.L[0]) + "\n"
		for row in self.L[1:]:
			 s+= 12*" " + str(row) + "\n"
		return s[:-1] + "."

def identity(n):
	'''
	Returns an n x n identity matrix.
	'''
	L = []
	for i in range(n):
		L += [[int(j == i) for j in range(n)]]
	return Matrix(L)

def concat(M,N):
	'''
	Input: M nxm matrix.
	       N n-1 x m-1 matrix.
	Output: A new matrix Q with N the
	 submatrix obtained by removing first col and row.
	 '''
	L = [M.L[0]]
	for i in range(1,len(M.L)):
		R = [M.L[i][0]]
		for j in range(len(N.L[0])):
			R += [N.L[i-1][j]]
		L += [R]
	return Matrix(L)

def submatrix(M,c=0,r=0):
	'''
	The submatrix obtained by removing col c and row r.
	'''
	L = []
	for i in range(len(M.L)):
		if i != r:
			R = []
			for j in range(len(M.L[0])):
				if j != c:
					R += [M.L[i][j]]
			L += [R]
	return Matrix(L)

Modulen “ntheory”, er en samling med tallteoretiske funksjoner jeg har skrevet tidligere. Den funksjonen jeg importerer regner ut største felles divisor mellom to tall a og b og returnerer en såkalt “Bezout-relasjon” mellom dem. Dette er et uttrykk på formen ax+by=d, hvor d er største felles divisor mellom a og $latex $b$. Koden til denne funksjonen ser slik ut:

def bezout(a,b):
    '''
    Returns a Bezout identity (as a tuple of two elts) for
    two numbers a,b. I.e. two numbers x,y such that
      ax+by=gcd(a,b)
    '''
    if b == 0:
        return (1,0)
    else:
        r = a % b
        q = (a-r)/b
        (s,t) = bezout(b,r)
        return (t,s-q*t)

Dette ble kanskje litt mye kode. Det viktigste er at dette er nok kode for å gjøre noen enkle men fundamentale operasjoner med heltallsmatriser. Nå er vi klare til å forklare litt om algoritmen bak funksjonen som transformerer en matrise til sin Hermite-normalform. Jeg skal ikke forklare det i detalj her, men ideen bak er dette: ved å multiplisere på venstresiden med elementærmatriser, kan vi på en systematisk måte fjerne elementene under diagonalen. Dette gjøres ved å finne det største elementet i første søyle, og bruke en Bezout-relasjon for å eliminere alle andre elementer.

Without further ado, her er koden:

def triangular(M):
	'''
	Input: an integer matrix M.
	Output: an upper triangulization U of M over the integers.
	Algorithm:
	 1. Find smallest nonzero element X in first column. Move this element to the top row.
	    If there are no nonzero elements go to step 3.
	 2. Use Bezout/Koszul-row operations to remove the nonzero elements below X. 
	 3. Repeat with the submatrix obtained by removing the first row and column.
	 4. Concatenate the result recursively.
	'''
	if len(M.L) == 1:
		if M.L[0][0] < 0:
			return -M
		return M
	(index, smallest) = (0,M.L[0][0])
	for i in range(len(M.L)):
		if abs(M.L[i][0])  < smallest:
			(index, smallest) = (i,M.L[i][0])
	if smallest == 0:
		return concat(M,triangular(submatrix(M)))
	if index != 0:
		N = M.switchRows(0,index)   ### row operation
	else:
		N = Matrix(M.L)
	for i in range(1,len(M.L)):
		d = ntheory.gcd(smallest, N.L[i][0])
		bez = ntheory.bezout(smallest, N.L[i][0])
		I = identity(len(N.L))
		I.L[0][0] = bez[0]
		I.L[0][i] = bez[1]
		I.L[i][0] = N.L[i][0]/d
		I.L[i][i] = -N.L[0][0]/d
		N = I * N
		if N.L[0][0] < 0:
			N = N.mult(-1)
	return concat(N,triangular(submatrix(N)))

Jeg skal la den interesserte og kapable leseren analysere funksjonen selv, og avslutter med noe noen eksempler på hvordan den fungerer.

La A være matrisen

A = \begin{pmatrix}1 & 2 \\ 3 & 4\end{pmatrix}.

Vi kjører da kommandoene “print(A)” og “print triangular(A)” i Python. Resultet er under:
2x2-matrix:
[1, 2]
[3, 4].
2x2-matrix:
[1, 2]
[0, 2].

Dermed kan vi se kokjernen er lik \mathbb Z/2. Nok et eksempel: La A være matrisen
A = \begin{pmatrix}1 & 2 & 3 \\ 4 & 5 & 6\end{pmatrix}

Resultatet blir da:

3x2-matrix:
[1, 2, 3]
[4, 5, 6].
3x2-matrix:
[1, 2, 3]
[0, 3, 6].

Dermed er kokjernen lik \mathbb Z/3 (det krever litt øvelse for å se sant umiddelbart, altså).

Til slutt en liten bemerkning: det finnes allerede matematisk programvare som regner ut ting som dette. For eksempel kan man i Macaulay2 skrive “prune coker A”, og svaret kommer med en gang. Hovedfordelen med å skrive slike funksjoner selv, er at en blir mer kjent med matematikken bak, og forståelsen blir større.

(dette ble en litt rotere bloggpost, men det får være. Hovedformålet var å subtilt skryte av at jeg bruke en lørdag til å implementere en interessant algoritme)