mutable-lattice 0.5.0

integer linear algebra using mutable sublattices of Z^n

mutable_lattice is a Python package that implements fast membership testing and mutation for sublattices of an integral lattice Z^n, with miscellaneous other integer linear algebra features.

Installation

Install this package with python -m pip install mutable_lattice.

Run the tests with python -m mutable_lattice.test.

Demo

The following constructs a sublattice of Z^3, adds some vectors to it, then tests which vectors are in the integer span of the previously added vectors.

>>> from mutable_lattice import Vector, Lattice
>>> L = Lattice(3)
>>> L.add_vector(Vector([2, 2, 2]))
>>> L.add_vector(Vector([2, 3, 3]))
>>> L.get_basis()
[Vector([2, 0, 0]), Vector([0, 1, 1])]
>>> Vector([-4, 7, 7]) in L
True
>>> Vector([3, 0, 0]) in L
False
>>> L.add_vector(Vector([3, 3, 3]))
>>> L.get_basis()
[Vector([1, 0, 0]), Vector([0, 1, 1])]
>>> Vector([1, 0, 0]) in L
True
>>> L.add_vector(Vector([1, 0, 1]))
>>> L.get_basis()
[Vector([1, 0, 0]), Vector([0, 1, 0]), Vector([0, 0, 1])]
>>> Vector([314, -159265, 3589793238462643383279]) in L
True

This behavior was originally implemented to help construct small projective resolutions for computing integral homology of finite monoids, but is useful for integer linear algebra in general.

There are more general integer linear algebra libraries included in SageMath (including PARI, GAP, and IML), and some of these may be preferred, especially when matrix entries include large integers, but this mutable_lattice package specializes in fast add_vector and __contains__ operations, and is especially fast with machine-word-sized integers. To this end, the mutable_lattice library presently does not use any algorithms based on modular arithmetic (e.g. Pernet and Stein 2009), and we instead rely only on (generalized) row operations.

Features

mutable_lattice exposes two types: Vector and Lattice.

Efficient Vectors of integers

The Vector class is implemented in C, and stores a sequence of integers, packed efficiently using tagged pointers to distinguish between machine-sized integers and pointers to Python int objects. This optimizes performance for integers that fit in one machine word (up to around 2^30 or 2^62), while also seamlessly intermixing integers of arbitrary size.

Vector arithmetic (addition, subtraction, negation, and integer scaling) behaves conventionally:

>>> v = Vector([10, 20, 30])
>>> w = Vector([7, 7, 7])
>>> v + w
Vector([17, 27, 37])
>>> v - w
Vector([3, 13, 23])
>>> 10*v
Vector([100, 200, 300])
>>> (-2) * v
Vector([-20, -40, -60])
>>> -w
Vector([-7, -7, -7])
>>> v * 10**20
Vector([1000000000000000000000, 2000000000000000000000, 3000000000000000000000])

One can access entries of a Vector with the __getitem__, tolist, and __iter__ methods:

>>> v = Vector([10, 20, 30])
>>> v[0]
10
>>> v[-1]
30
>>> v.tolist()
[10, 20, 30]
>>> [-x for x in v]
[-10, -20, -30]
>>> len(v)
3

Vectors are mutable and have __iadd__, __isub__, __imul__, and __setitem__ methods:

>>> v = Vector([1, 1, 1, 1, 1])
>>> v[1] = 4
>>> v
Vector([1, 4, 1, 1, 1])
>>> v += Vector([1, 1, 1, 1, 1])
>>> v
Vector([2, 5, 2, 2, 2])
>>> v *= 2
>>> v
Vector([4, 10, 4, 4, 4])
>>> v -= Vector([4, 4, 4, 4, 4])
>>> v
Vector([0, 6, 0, 0, 0])

Vector.zero(n) is equivalent to Vector([0]*n).

Efficient Lattices of integer vectors: Overview

The Lattice class by default stores a basis of vectors in Hermite normal form (HNF), i.e., the integer version of reduced row echelon form.

The Lattice.get_basis() method returns a list of basis Vectors for the Lattice. The Lattice.__str__() method prints this basis as the rows of a matrix. Each time a new vector is added via Lattice.add_vector(v), HNF is restored.

>>> L = Lattice(5)
>>> print(L)
<zero Lattice in Z^5>

>>> L.add_vector(Vector([1, 1, 1, 1, 1]))
>>> L.get_basis()
[Vector([1, 1, 1, 1, 1])]
>>> print(L)
[1 1 1 1 1]

>>> L.add_vector(Vector([10, 0, 10, 0, 10]))
>>> L.get_basis()
[Vector([1, 1, 1, 1, 1]), Vector([0, 10, 0, 10, 0])]
>>> print(L)
[ 1  1  1  1  1]
[ 0 10  0 10  0]

>>> L.add_vector(Vector([1, 0, 0, 0, 0]))
>>> L.get_basis()
[Vector([1, 0, 0, 0, 0]), Vector([0, 1, 1, 1, 1]), Vector([0, 0, 10, 0, 10])]
>>> print(L)
[ 1  0  0  0  0]
[ 0  1  1  1  1]
[ 0  0 10  0 10]

The Lattice.rank attribute gives the number of vectors currently stored in a Lattice, and the Lattice.ambient_dimension attribute gives the length of each vector. So L.__str__() is displaying a matrix of height L.rank and width L.ambient_dimension:

>>> print(L) # Continuing from above
[ 1  0  0  0  0]
[ 0  1  1  1  1]
[ 0  0 10  0 10]
>>> L.rank
3
>>> L.ambient_dimension
5

The Lattice.__contains__(v) method identifies whether the argument Vector is in the Lattice, i.e., whether v is in the integer span of the stored basis:

>>> print(L) # Continuing from above
[ 1  0  0  0  0]
[ 0  1  1  1  1]
[ 0  0 10  0 10]
>>> Vector([777, 1, 11,  1, 11]) in L
True
>>> Vector([0, 0, 10, 10, 10]) in L
False

To reconstruct the linear combination of the Lattice basis that produces the given vector, use the Lattice.coefficients_of(v) method. For the inverse operation, use the Lattice.linear_combination(w) method to evaluate a linear combination of the Lattice basis given a Vector of coefficients w:

>>> print(L) # Continuing from above
[ 1  0  0  0  0]
[ 0  1  1  1  1]
[ 0  0 10  0 10]
>>> L.coefficients_of(Vector([0, 0, 10, 10, 10]))
Traceback (most recent call last):
    ...
ValueError: Vector not present in Lattice
>>> L.coefficients_of(Vector([777, 1, 11, 1, 11]))
Vector([777, 1, 1])
>>> L.linear_combination(Vector([777, 1, 1]))
Vector([777, 1, 11, 1, 11])

Lattice.__add__ and Lattice.__iadd__ methods allow adding two lattices of the same ambient dimension, so L1 + L2 consists of all vectors v1 + v2 where v1 in L1 and v2 in L2.

>>> L1 = Lattice(4)
>>> L1.add_vector(Vector([2, 0, 0, 0]))
>>> L1.add_vector(Vector([0, 2, 0, 0]))
>>> L1.add_vector(Vector([0, 0, 2, 0]))
>>> L2 = Lattice(4)
>>> L2.add_vector(Vector([0, 4, 0, 0]))
>>> L2.add_vector(Vector([0, 0, 4, 0]))
>>> L2.add_vector(Vector([0, 0, 0, 4]))
>>> print(L1 + L2)
[2 0 0 0]
[0 2 0 0]
[0 0 2 0]
[0 0 0 4]

Comparison methods such as L1 < L2, L1 <= L2, L1 == L2 are also supported for detecting when one lattice is a subset of another. The expression L1 & L2 computes the intersection of the two Lattices.

Lattice(n, data=..., /, *, maxrank=-1, HNF_policy=1)

The Lattice(n, data) constructor creates a sublattice of Z^n and accepts an optional second positional argument data which must be a list of any length, of Vectors or lists of length n, which are added to the lattice as if by calling add_vector(v) or add_vector(Vector(v)) for each v in data.

For memory locality and to avoid reallocations, Lattice objects are allocated in one contiguous memory block with enough room for all n*n integer entries. If n is large, this may be needlessly memory intensive. If the optional keyword-only integer argument maxrank is provided then only enough memory is allocated for n*maxrank integer entries, potentially saving memory. If maxrank > n or if maxrank == -1 then maxrank is replaced by n.

For more uniformly fast performance in more situations, by default all Lattices are stored as a basis in HNF form. If a Lattice is constructed instead with HNF_policy=0 then the basis of the Lattice is stored in some row echelon form, but the pivots are not necessarily positive, and entries of the matrix above pivots are not normalized to HNF. This lazier policy 0 may be faster if your basis is particularly sparse or small (say n <= 10 but always be sure to measure for your particular application), but the required row operations can otherwise quickly cause integers to explode in size, which is partially mitigated by using HNF_policy=1. If you have a Lattice not stored in HNF then calling L.HNFify() will perform row operations to convert to HNF.

The full lattice that contains every vector in Z^n can be constructed using the classmethod Lattice.full(n). L.is_full() returns whether L == Lattice.full(L.ambient_dimension).

SNF invariants

Given a Lattice L in Z^n, some (integer-)invertible n-by-n integer matrix takes L to a Lattice with a basis [Vector([d0, 0, 0, ...]), Vector([0, d1, 0, ...]), Vector([0, 0, d2, ...]), ...] in which the ith Vector is positive in its ith entry and zero elsewhere. The matrix with these rows (and an additional divisibility constraint) is the Smith normal form (SNF) of the matrix with the original basis of L for rows. The Lattice.nonzero_invariants() method returns a list [d0, d1, d2, ...] of these diagonal entries of the Smith normal form. The length of L.nonzero_invariants() is always L.rank. A related method L.invariants() returns the same list with zeros appended to the end until the length of the result is L.ambient_dimension.

>>> L = Lattice(4)
>>> L.add_vector(Vector([10, 10, 10, 10]))
>>> L.add_vector(Vector([0, 20, 20, 20]))
>>> L.add_vector(Vector([0, 0, 30, 30]))
>>> L.nonzero_invariants()
[10, 10, 60]
>>> L.invariants()
[10, 10, 60, 0]

Note that the quotient group Z^n/L is isomorphic to the direct sum of cyclic groups Z/dZ for d in L.invariants(), where Z/0Z = Z.

The (co)homology of a chain complex of integer matrices between free abelian groups can be computed using this method:

>>> # Chain complex from a Klein bottle CW complex
>>> cell_counts = {-1: 0, 0: 1, 1: 2, 2: 1, 3: 0}
>>> boundary = {
...     0: [Vector([])],
...     1: [Vector([0]), Vector([0])],
...     2: [Vector([2, 0])],
...     3: [],
... }
>>> invariant_lists = {
...    i: Lattice(cell_counts[i-1], M).nonzero_invariants()
...    for i, M in boundary.items()
... }
>>> for i in (0, 1, 2):
...     free_rank = cell_counts[i] - len(invariant_lists[i]) - len(invariant_lists[i+1])
...     torsion = [d for d in invariant_lists[i+1] if d != 1]
...     print(f"H_{i} =", torsion + [0] * free_rank)
H_0 = [0]
H_1 = [2, 0]
H_2 = []
>>> for i in (0, 1, 2):
...     free_rank = cell_counts[i] - len(invariant_lists[i]) - len(invariant_lists[i+1])
...     torsion = [d for d in invariant_lists[i] if d != 1]
...     print(f"H^{i} =", torsion + [0] * free_rank)
H^0 = [0]
H^1 = [0]
H^2 = [2]

We do not currently provide functionality to produce the unimodular matrices S, T that put a matrix A into its Smith normal form SAT.

Kernel computations

Given a list of k vectors, the relations_among([v1, ..., vk]) function constructs the Lattice of linear dependencies among these vectors.

>>> from mutable_lattice import relations_among
>>> L = relations_among([Vector([1, 2]), Vector([-1, 0]), Vector([0, 1]), Vector([0, 0])])
>>> print(L)
[ 1  1 -2  0]
[ 0  0  0  1]

Thinking of the provided list [v1, ..., vk] as a matrix M of row vectors, the result of relations_among() is the Lattice of vectors w such that wM=0, i.e, the left-kernel of M.

To instead compute the classic right-kernel that solves Mw=0, first transpose the matrix with the transpose(n, [v1, ..., vk]) function (which requires each Vector in the list to have length n):

>>> from mutable_lattice import transpose
>>> M = [Vector([1, -1, 0, 0]), Vector([2, 0, 1, 0])]
>>> print(relations_among(M))
<zero Lattice in Z^2>
>>> transpose(4, M)
[Vector([1, 2]), Vector([-1, 0]), Vector([0, 1]), Vector([0, 0])]
>>> print(relations_among(transpose(4, M)))
[ 1  1 -2  0]
[ 0  0  0  1]

v.shuffled_by_action(a, result_size=..., /)

It is occasionally useful to shuffle around the entries of a Vector by a permutation:

>>> perm = Vector([0, 2, 4, 1, 3])
>>> Vector([0, 10, 20, 30, 40]).shuffled_by_action(perm)
Vector([0, 30, 10, 40, 20])

Note that each entry perm[i] above specifies where the ith entry moves to after the shuffle.

More generally, we can shuffle the entries of the Vector Using some function from range(len(v)) to range(result_size), represented as a "action" Vector. If result_size is not provided, it defaults to len(v). If the action a has duplicate entries a[j1] == a[j2] == ..., then v.shuffled_by_action(a)[a[j1]] will be their sum v[j1] + v[j2] + ..., so each action Vector specifies a linear operation.

>>> Vector([0, 10, 20, 30, 40]).shuffled_by_action(Vector([0, 1, 0, 1, 0]))
Vector([60, 40, 0, 0, 0])
>>> Vector([0, 10, 20, 30, 40]).shuffled_by_action(Vector([0, 1, 0, 1, 0]), 2)
Vector([60, 40])

This package does not provide a built-in way of applying more general linear operations on Vectors via matrices, nor a way to multiply matrices, but these can be emulated using appropriate Vector addition and scaling.

L.decompose(keep_together=[], /)

It is sometimes useful to decompose a lattice into a direct sum of smaller lattices, shuffled together into disjoint rows and columns:

>>> L = Lattice(6, [[10,0,2,0,0,0], [0,20,0,5,1,0], [0,0,0,30,2,0], [0,0,0,0,0,3]])
>>> print(L)
[10  0  2  0  0  0]
[ 0 20  0  5  1  0]
[ 0  0  0 30  2  0]
[ 0  0  0  0  0  3]
>>> indexes, summands = L.decompose()
>>> indexes
[[0, 2], [1, 3, 4], [5]]
>>> print(summands[0])
[10  2]
>>> print(summands[1])
[20  5  1]
[ 0 30  2]
>>> print(summands[2])
[3]

The result of L.decompose() is a pair (indexes, summands) of two lists of the same length. The indexes list contains a partition of range(L.ambient_dimension), stored as a list of lists of integers. Each list of integers is in ascending order, but the ordering between the lists in the indexes list is undefined. The list entry summands[i] is the result of restricting L only to the indexes (columns) defined by indexes[i].

By default, L.decompose() decomposes L into as many of these direct summands as possible. However, it is also possible to require some indices to be placed into the same direct summand:

>>> L = Lattice(5, [[1, 0, 0, 0, 0], [0, 2, 0, 0, 0], [0, 0, 3, 0, 0], [0, 0, 0, 4, 4]])
>>> indexes, summands = L.decompose([[0, 1], [2, 3]])
>>> indexes
[[0, 1], [2, 3, 4]]
>>> summands[0]
Lattice(2, [[1, 0], [0, 2]])
>>> summands[1]
Lattice(3, [[3, 0, 0], [0, 4, 4]], maxrank=2)

The single optional argument of L.decompose(keep_together) should be an iterable of "component" iterables, and the integers in each component iterable are then ensured to be placed in the same component. The summands are allocated "tightly", with maxrank equal to their actual rank.