Sparse linear algebra

Introduction

This chapter deals with sparse linear algebra over commutative rings and fields.

Sparse linear algebra, that is, linear algebra with sparse matrices, plays an important role in various algorithms in algebraic number theory. For example, it is one of the key ingredients in the computation of class groups and discrete logarithms using index calculus methods.

Sparse rows

Building blocks for sparse matrices are sparse rows, which are modelled by objects of type SRow. More precisely, the type is of parametrized form SRow{T}, where T is the element type of the base ring $R$ . For example, SRow{ZZRingElem} is the type for sparse rows over the integers.

It is important to note that sparse rows do not have a fixed number of columns, that is, they represent elements of ${(x_{i})_{i} \in R^{N} ∣ x_{i} = 0 for almost all i}$ . In particular any two sparse rows over the same base ring can be added.

Creation

sparse_row Method

julia

sparse_row(R::Ring, J::Vector{Tuple{Int, T}}) -> SRow{T}

Constructs the sparse row $(a_{i})_{i}$ with $a_{i_{j}} = x_{j}$ , where $J = (i_{j}, x_{j})_{j}$ . The elements $x_{i}$ must belong to the ring $R$ .

sparse_row Method

julia

sparse_row(R::Ring, J::Vector{Tuple{Int, Int}}) -> SRow

Constructs the sparse row $(a_{i})_{i}$ over $R$ with $a_{i_{j}} = x_{j}$ , where $J = (i_{j}, x_{j})_{j}$ .

sparse_row Method

julia

sparse_row(R::NCRing, J::Vector{Int}, V::Vector{T}) -> SRow{T}

Constructs the sparse row $(a_{i})_{i}$ over $R$ with $a_{i_{j}} = x_{j}$ , where $J = (i_{j})_{j}$ and $V = (x_{j})_{j}$ .

Basic operations

Rows support the usual operations:

+, -, ==
multiplication by scalars
div, divexact

getindex Method

julia

getindex(A::SRow, j::Int) -> RingElem

Given a sparse row $(a_{i})_{i}$ and an index $j$ return $a_{j}$ .

add_scaled_row Method

julia

add_scaled_row(A::SRow{T}, B::SRow{T}, c::T) -> SRow{T}

Returns the row $c A + B$ .

add_scaled_row Method

julia

add_scaled_row(A::SRow{T}, B::SRow{T}, c::T) -> SRow{T}

Returns the row $c A + B$ .

transform_row Method

julia

transform_row(A::SRow{T}, B::SRow{T}, i::Int, j::Int, a::T, b::T, c::T, d::T)

Returns the tuple $(a A + b B, c A + d B)$ .

length Method

julia

length(A::SRow)

Returns the number of nonzero entries of $A$ .

Change of base ring

change_base_ring Method

julia

change_base_ring(R::Ring, A::SRow) -> SRow

Create a new sparse row by coercing all elements into the ring $R$ .

Maximum, minimum and 2-norm

maximum Method

julia

maximum(A::SRow{T}) -> T

Returns the largest entry of $A$ .

maximum Method

julia

maximum(A::SRow{T}) -> T

Returns the largest entry of $A$ .

minimum Method

julia

minimum(A::SRow{T}) -> T

Returns the smallest entry of $A$ .

julia

  minimum(A::RelNumFieldOrderIdeal) -> AbsNumFieldOrderIdeal{AbsSimpleNumField, AbsSimpleNumFieldElem}
  minimum(A::RelNumFieldOrderIdeal) -> RelNumFieldOrderIdeal

Returns the ideal $A \cap O$ where $O$ is the maximal order of the coefficient ideals of $A$ .

minimum Method

julia

minimum(A::SRow{T}) -> T

Returns the smallest entry of $A$ .

norm2 Method

julia

norm2(A::SRow{T} -> T

Returns $A \cdot A^{t}$ .

Functionality for integral sparse rows

lift Method

julia

lift(A::SRow{zzModRingElem}) -> SRow{ZZRingElem}

Return the sparse row obtained by lifting all entries in $A$ .

mod! Method

julia

mod!(A::SRow{ZZRingElem}, n::ZZRingElem) -> SRow{ZZRingElem}

Inplace reduction of all entries of $A$ modulo $n$ to the positive residue system.

mod_sym! Method

julia

mod_sym!(A::SRow{ZZRingElem}, n::ZZRingElem) -> SRow{ZZRingElem}

Inplace reduction of all entries of $A$ modulo $n$ to the symmetric residue system.

mod_sym! Method

julia

mod_sym!(A::SRow{ZZRingElem}, n::Integer) -> SRow{ZZRingElem}

Inplace reduction of all entries of $A$ modulo $n$ to the symmetric residue system.

maximum Method

julia

maximum(abs, A::SRow{ZZRingElem}) -> ZZRingElem

Returns the largest, in absolute value, entry of $A$ .

Conversion to/from julia and AbstractAlgebra types

Vector Method

julia

Vector(a::SMat{T}, n::Int) -> Vector{T}

The first n entries of a, as a julia vector.

sparse_row Method

julia

sparse_row(A::MatElem)

Convert A to a sparse row. nrows(A) == 1 must hold.

dense_row Method

julia

dense_row(r::SRow, n::Int)

Convert r[1:n] to a dense row, that is an AbstractAlgebra matrix.

Sparse matrices

Let $R$ be a commutative ring. Sparse matrices with base ring $R$ are modelled by objects of type SMat. More precisely, the type is of parametrized form SRow{T}, where T is the element type of the base ring. For example, SMat{ZZRingElem} is the type for sparse matrices over the integers.

In contrast to sparse rows, sparse matrices have a fixed number of rows and columns, that is, they represent elements of the matrices space ${Mat}_{n \times m} (R)$ . Internally, sparse matrices are implemented as an array of sparse rows. As a consequence, unlike their dense counterparts, sparse matrices have a mutable number of rows and it is very performant to add additional rows.

Construction

sparse_matrix Method

julia

sparse_matrix(R::Ring) -> SMat

Return an empty sparse matrix with base ring $R$ .

sparse_matrix Method

julia

sparse_matrix(R::Ring, n::Int, m::Int) -> SMat

Return a sparse $n$ times $m$ zero matrix over $R$ .

Sparse matrices can also be created from dense matrices as well as from julia arrays:

sparse_matrix Method

julia

sparse_matrix(A::MatElem; keepzrows::Bool = true)

Constructs the sparse matrix corresponding to the dense matrix $A$ . If keepzrows is false, then the constructor will drop any zero row of $A$ .

sparse_matrix Method

julia

sparse_matrix(R::Ring, A::Matrix{T}) -> SMat

Constructs the sparse matrix over $R$ corresponding to $A$ .

sparse_matrix Method

julia

sparse_matrix(R::Ring, A::Matrix{T}) -> SMat

Constructs the sparse matrix over $R$ corresponding to $A$ .

The normal way however, is to add rows:

push! Method

julia

push!(A::SMat{T}, B::SRow{T}) where T

Appends the sparse row B to A.

Sparse matrices can also be concatenated to form larger ones:

vcat! Method

julia

vcat!(A::SMat, B::SMat) -> SMat

Vertically joins $A$ and $B$ inplace, that is, the rows of $B$ are appended to $A$ .

vcat Method

julia

vcat(A::SMat, B::SMat) -> SMat

Vertically joins $A$ and $B$ .

hcat! Method

julia

hcat!(A::SMat, B::SMat) -> SMat

Horizontally concatenates $A$ and $B$ , inplace, changing $A$ .

hcat Method

julia

hcat(A::SMat, B::SMat) -> SMat

Horizontally concatenates $A$ and $B$ .

(Normal julia $c a t$ is also supported)

There are special constructors:

identity_matrix Method

julia

identity_matrix(::Type{SMat}, R::Ring, n::Int)

Return a sparse $n$ times $n$ identity matrix over $R$ .

zero_matrix Method

julia

zero_matrix(::Type{SMat}, R::Ring, n::Int)

Return a sparse $n$ times $n$ zero matrix over $R$ .

zero_matrix Method

julia

zero_matrix(::Type{SMat}, R::Ring, n::Int, m::Int)

Return a sparse $n$ times $m$ zero matrix over $R$ .

block_diagonal_matrix Method

julia

block_diagonal_matrix(xs::Vector{SMat})

Return the block diagonal matrix with the matrices in xs on the diagonal. Requires all blocks to have the same base ring.

Slices:

sub Method

julia

sub(A::SMat, r::AbstractUnitRange, c::AbstractUnitRange) -> SMat

Return the submatrix of $A$ , where the rows correspond to $r$ and the columns correspond to $c$ .

Transpose:

transpose Method

julia

transpose(A::SMat) -> SMat

Returns the transpose of $A$ .

Elementary Properties

sparsity Method

julia

sparsity(A::SMat) -> Float64

Return the sparsity of A, that is, the number of zero-valued elements divided by the number of all elements.

density Method

julia

density(A::SMat) -> Float64

Return the density of A, that is, the number of nonzero-valued elements divided by the number of all elements.

nnz Method

julia

nnz(A::SMat) -> Int

Return the number of non-zero entries of $A$ .

number_of_rows Method

julia

number_of_rows(A::SMat) -> Int

Return the number of rows of $A$ .

number_of_columns Method

julia

number_of_columns(A::SMat) -> Int

Return the number of columns of $A$ .

isone Method

julia

isone(A::SMat)

Tests if $A$ is an identity matrix.

iszero Method

julia

iszero(A::SMat)

Tests if $A$ is a zero matrix.

is_upper_triangular Method

julia

is_upper_triangular(A::SMat)

Returns true if and only if $A$ is upper (right) triangular.

maximum Method

julia

maximum(A::SMat{T}) -> T

Finds the largest entry of $A$ .

minimum Method

julia

minimum(A::SMat{T}) -> T

Finds the smallest entry of $A$ .

maximum Method

julia

maximum(abs, A::SMat{ZZRingElem}) -> ZZRingElem

Finds the largest, in absolute value, entry of $A$ .

elementary_divisors Method

julia

elementary_divisors(A::SMat{ZZRingElem}) -> Vector{ZZRingElem}

The elementary divisors of $A$ , i.e. the diagonal elements of the Smith normal form of $A$ .

solve_dixon_sf Method

julia

solve_dixon_sf(A::SMat{ZZRingElem}, b::SRow{ZZRingElem}, is_int::Bool = false) -> SRow{ZZRingElem}, ZZRingElem
solve_dixon_sf(A::SMat{ZZRingElem}, B::SMat{ZZRingElem}, is_int::Bool = false) -> SMat{ZZRingElem}, ZZRingElem

For a sparse square matrix $A$ of full rank and a sparse matrix (row), find a sparse matrix (row) $x$ and an integer $d$ s.th. $x A = b d$ holds. The algorithm is a Dixon-based linear p-adic lifting method. If \code{is_int} is given, then $d$ is assumed to be $1$ . In this case rational reconstruction is avoided.

hadamard_bound2 Method

julia

hadamard_bound2(A::SMat{T}) -> T

The square of the product of the norms of the rows of $A$ .

echelon_with_transform Method

julia

echelon_with_transform(A::SMat{zzModRingElem}) -> SMat, SMat

Find a unimodular matrix $T$ and an upper-triangular $E$ s.th. $T A = E$ holds.

reduce_full Method

julia

reduce_full(A::SMat{ZZRingElem}, g::SRow{ZZRingElem},
                      with_transform = Val(false)) -> SRow{ZZRingElem}, Vector{Int}

Reduces $g$ modulo $A$ and assumes that $A$ is upper triangular.

The second return value is the array of pivot elements of $A$ that changed.

If with_transform is set to Val(true), then additionally an array of transformations is returned.

hnf! Method

julia

hnf!(A::SMat{ZZRingElem})

Inplace transform of $A$ into upper right Hermite normal form.

hnf Method

julia

hnf(A::SMat{ZZRingElem}) -> SMat{ZZRingElem}

Return the upper right Hermite normal form of $A$ .

snf Method

julia

snf(A::SMat{ZZRingElem})

The Smith normal form (snf) of $A$ .

hnf_extend! Method

julia

hnf_extend!(A::SMat{ZZRingElem}, b::SMat{ZZRingElem}, offset::Int = 0) -> SMat{ZZRingElem}

Given a matrix $A$ in HNF, extend this to get the HNF of the concatenation with $b$ .

is_diagonal Method

julia

is_diagonal(A::SMat) -> Bool

True iff only the i-th entry in the i-th row is non-zero.

det Method

julia

det(A::SMat{ZZRingElem})

The determinant of $A$ using a modular algorithm. Uses the dense (zzModMatrix) determinant on $A$ for various primes $p$ .

det_mc Method

julia

det_mc(A::SMat{ZZRingElem})

Computes the determinant of $A$ using a LasVegas style algorithm, i.e. the result is not proven to be correct. Uses the dense (zzModMatrix) determinant on $A$ for various primes $p$ .

valence_mc Method

julia

valence_mc{T}(A::SMat{T}; extra_prime = 2, trans = Vector{SMatSLP_add_row{T}}()) -> T

Uses a Monte-Carlo algorithm to compute the valence of $A$ . The valence is the valence of the minimal polynomial $f$ of $t r a n s p o s e (A) * A$ , thus the last non-zero coefficient, typically $f (0)$ .

The valence is computed modulo various primes until the computation stabilises for extra_prime many.

trans, if given, is a SLP (straight-line-program) in GL(n, Z). Then the valence of trans * $A$ is computed instead.

saturate Method

julia

saturate(A::SMat{ZZRingElem}) -> SMat{ZZRingElem}

Computes the saturation of $A$ , that is, a basis for $Q \otimes M \cap Z^{n}$ , where $M$ is the row span of $A$ and $n$ the number of rows of $A$ .

Equivalently, return $T A$ for an invertible rational matrix $T$ , such that $T A$ is integral and the elementary divisors of $T A$ are all trivial.

hnf_kannan_bachem Method

julia

hnf_kannan_bachem(A::SMat{ZZRingElem}) -> SMat{ZZRingElem}

Compute the Hermite normal form of $A$ using the Kannan-Bachem algorithm.

diagonal_form Method

julia

diagonal_form(A::SMat{ZZRingElem}) -> SMat{ZZRingElem}

A matrix $D$ that is diagonal and obtained via unimodular row and column operations. Like a snf without the divisibility condition.

Manipulation/ Access

getindex Method

julia

getindex(A::SMat, i::Int, j::Int)

Given a sparse matrix $A = (a_{i j})_{i, j}$ , return the entry $a_{i j}$ .

getindex Method

julia

getindex(A::SMat, i::Int) -> SRow

Given a sparse matrix $A$ and an index $i$ , return the $i$ -th row of $A$ .

setindex! Method

julia

setindex!(A::SMat, b::SRow, i::Int)

Given a sparse matrix $A$ , a sparse row $b$ and an index $i$ , set the $i$ -th row of $A$ equal to $b$ .

swap_rows! Method

julia

swap_rows!(A::SMat{T}, i::Int, j::Int)

Swap the $i$ -th and $j$ -th row of $A$ inplace.

swap_cols! Method

julia

swap_cols!(A::SMat, i::Int, j::Int)

Swap the $i$ -th and $j$ -th column of $A$ inplace.

scale_row! Method

julia

scale_row!(A::SMat{T}, i::Int, c::T)

Multiply the $i$ -th row of $A$ by $c$ inplace.

add_scaled_col! Method

julia

add_scaled_col!(A::SMat{T}, i::Int, j::Int, c::T)

Add $c$ times the $i$ -th column to the $j$ -th column of $A$ inplace, that is, $A_{j} \to A_{j} + c \cdot A_{i}$ , where $(A_{i})_{i}$ denote the columns of $A$ .

add_scaled_row! Method

julia

add_scaled_row!(A::SMat{T}, i::Int, j::Int, c::T)

Add $c$ times the $i$ -th row to the $j$ -th row of $A$ inplace, that is, $A_{j} \to A_{j} + c \cdot A_{i}$ , where $(A_{i})_{i}$ denote the rows of $A$ .

transform_row! Method

julia

transform_row!(A::SMat{T}, i::Int, j::Int, a::T, b::T, c::T, d::T)

Applies the transformation $(A_{i}, A_{j}) \to (a A_{i} + b A_{j}, c A_{i} + d A_{j})$ to $A$ , where $(A_{i})_{i}$ are the rows of $A$ .

diagonal Method

julia

diagonal(A::SMat) -> ZZRingElem[]

The diagonal elements of $A$ in an array.

reverse_rows! Method

julia

reverse_rows!(A::SMat)

Inplace inversion of the rows of $A$ .

mod_sym! Method

julia

mod_sym!(A::SMat{ZZRingElem}, n::ZZRingElem)

Inplace reduction of all entries of $A$ modulo $n$ to the symmetric residue system.

find_row_starting_with Method

julia

find_row_starting_with(A::SMat, p::Int) -> Int

Tries to find the index $i$ such that $A_{i, p} \neq 0$ and $A_{i, p - j} = 0$ for all $j > 1$ . It is assumed that $A$ is upper triangular. If such an index does not exist, find the smallest index larger.

reduce Method

julia

reduce(A::SMat{ZZRingElem}, g::SRow{ZZRingElem}, m::ZZRingElem) -> SRow{ZZRingElem}

Given an upper triangular matrix $A$ over the integers, a sparse row $g$ and an integer $m$ , this function reduces $g$ modulo $A$ and returns $g$ modulo $m$ with respect to the symmetric residue system.

reduce Method

julia

reduce(A::SMat{ZZRingElem}, g::SRow{ZZRingElem}) -> SRow{ZZRingElem}

Given an upper triangular matrix $A$ over a field and a sparse row $g$ , this function reduces $g$ modulo $A$ .

reduce Method

julia

reduce(A::SMat{T}, g::SRow{T}) -> SRow{T}

Given an upper triangular matrix $A$ over a field and a sparse row $g$ , this function reduces $g$ modulo $A$ .

rand_row Method

julia

rand_row(A::SMat) -> SRow

Return a random row of the sparse matrix $A$ .

Changing of the ring:

map_entries Method

julia

map_entries(f, A::SMat) -> SMat

Given a sparse matrix $A$ and a callable object $f$ , this function will construct a new sparse matrix by applying $f$ to all elements of $A$ .

change_base_ring Method

julia

change_base_ring(R::Ring, A::SMat)

Create a new sparse matrix by coercing all elements into the ring $R$ .

Arithmetic

Matrices support the usual operations as well

+, -, ==
div, divexact by scalars
multiplication by scalars

Various products:

* Method

julia

*(A::SMat{T}, b::AbstractVector{T}) -> Vector{T}

Return the product $A \cdot b$ as a dense vector.

* Method

julia

*(A::SMat{T}, b::AbstractMatrix{T}) -> Matrix{T}

Return the product $A \cdot b$ as a dense array.

* Method

julia

*(A::SMat{T}, b::MatElem{T}) -> MatElem

Return the product $A \cdot b$ as a dense matrix.

* Method

julia

*(A::SRow, B::SMat) -> SRow

Return the product $A \cdot B$ as a sparse row.

dot Method

julia

dot(x::SRow{T}, A::SMat{T}, y::SRow{T}) where T -> T

Return the generalized dot product dot(x, A*y).

dot Method

julia

dot(x::MatrixElem{T}, A::SMat{T}, y::MatrixElem{T}) where T -> T

Return the generalized dot product dot(x, A*y).

dot Method

julia

dot(x::AbstractVector{T}, A::SMat{T}, y::AbstractVector{T}) where T -> T

Return the generalized dot product dot(x, A*y).

Other:

sparse Method

julia

sparse(A::SMat) -> SparseMatrixCSC

The same matrix, but as a sparse matrix of julia type SparseMatrixCSC.

ZZMatrix Method

julia

ZZMatrix(A::SMat{ZZRingElem})

The same matrix $A$ , but as an ZZMatrix.

ZZMatrix Method

julia

ZZMatrix(A::SMat{T}) where {T <: Integer}

The same matrix $A$ , but as an ZZMatrix. Requires a conversion from the base ring of $A$ to $Z Z$ .

Matrix Method

julia

Matrix(A::SMat{T}) -> Matrix{T}

The same matrix, but as a julia matrix.

Array Method

julia

Array(A::SMat{T}) -> Matrix{T}

The same matrix, but as a two-dimensional julia array.