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^{\mathbb{N}}\mid x_i=0\text{ for almost all }i\}$. In particular any two sparse rows over the same base ring can be added.

Creation

# sparse_row — Method.

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$.

source

# sparse_row — Method.

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$.

source

# sparse_row — Method.

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$.

source

Basic operations

Rows support the usual operations:

• +, -, ==

• multiplication by scalars

• div, divexact

# getindex — Method.

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

Given a sparse row $(a_i{)}_i$ and an index $j$ return $a_j$.

source

# add_scaled_row — Method.

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

Returns the row $cA+B$.

source

# add_scaled_row — Method.

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

Returns the row $cA+B$.

source

# transform_row — Method.

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

Returns the tuple $(aA+bB,cA+dB)$.

source

# length — Method.

length(A::SRow)

Returns the number of nonzero entries of $A$.

source

Change of base ring

# change_base_ring — Method.

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

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

source

Maximum, minimum and 2-norm

# maximum — Method.

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

Returns the largest entry of $A$.

source

# maximum — Method.

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

Returns the largest entry of $A$.

source

# minimum — Method.

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

Returns the smallest entry of $A$.

source

  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$.

source

# minimum — Method.

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

Returns the smallest entry of $A$.

source

# norm2 — Method.

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

Returns $A\cdot A^t$.

source

Functionality for integral sparse rows

# lift — Method.

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

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

source

# mod! — Method.

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

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

source

# mod_sym! — Method.

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

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

source

# mod_sym! — Method.

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

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

source

# maximum — Method.

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

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

source

Conversion to/from julia and AbstractAlgebra types

# Vector — Method.

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

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

source

# sparse_row — Method.

sparse_row(A::MatElem)

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

source

# dense_row — Method.

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

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

source

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 ${\mathrm{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.

sparse_matrix(R::Ring) -> SMat

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

source

# sparse_matrix — Method.

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

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

source

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

# sparse_matrix — Method.

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$.

source

# sparse_matrix — Method.

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

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

source

# sparse_matrix — Method.

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

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

source

The normal way however, is to add rows:

# push! — Method.

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

Appends the sparse row B to A.

source

Sparse matrices can also be concatenated to form larger ones:

# vcat! — Method.

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

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

source

# vcat — Method.

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

Vertically joins $A$ and $B$.

source

# hcat! — Method.

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

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

source

# hcat — Method.

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

Horizontally concatenates $A$ and $B$.

source

(Normal julia $cat$ is also supported)

There are special constructors:

# identity_matrix — Method.

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

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

source

# zero_matrix — Method.

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

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

source

# zero_matrix — Method.

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

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

source

# block_diagonal_matrix — Method.

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.

source

Slices:

# sub — Method.

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$.

source

Transpose:

# transpose — Method.

transpose(A::SMat) -> SMat

Returns the transpose of $A$.

source

Elementary Properties

# sparsity — Method.

sparsity(A::SMat) -> Float64

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

source

# density — Method.

density(A::SMat) -> Float64

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

source

# nnz — Method.

nnz(A::SMat) -> Int

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

source

# number_of_rows — Method.

number_of_rows(A::SMat) -> Int

Return the number of rows of $A$.

source

# number_of_columns — Method.

number_of_columns(A::SMat) -> Int

Return the number of columns of $A$.

source

# isone — Method.

isone(A::SMat)

Tests if $A$ is an identity matrix.

source

# iszero — Method.

iszero(A::SMat)

Tests if $A$ is a zero matrix.

source

# is_upper_triangular — Method.

is_upper_triangular(A::SMat)

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

source

# maximum — Method.

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

Finds the largest entry of $A$.

source

# minimum — Method.

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

Finds the smallest entry of $A$.

source

# maximum — Method.

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

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

source

# elementary_divisors — Method.

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

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

source

# solve_dixon_sf — Method.

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. $xA=bd$ 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.

source

# hadamard_bound2 — Method.

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

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

source

# echelon_with_transform — Method.

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

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

source

# reduce_full — Method.

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.

source

# hnf! — Method.

hnf!(A::SMat{ZZRingElem})

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

source

# hnf — Method.

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

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

source

# snf — Method.

snf(A::SMat{ZZRingElem})

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

source

# hnf_extend! — Method.

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$.

source

# is_diagonal — Method.

is_diagonal(A::SMat) -> Bool

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

source

# det — Method.

det(A::SMat{ZZRingElem})

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

source

# det_mc — Method.

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$.

source

# valence_mc — Method.

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 $transpose(A)\ast 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.

source

# saturate — Method.

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

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

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

source

# hnf_kannan_bachem — Method.

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

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

source

# diagonal_form — Method.

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.

source

Manipulation/ Access

# getindex — Method.

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

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

source

# getindex — Method.

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

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

source

# setindex! — Method.

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$.

source

# swap_rows! — Method.

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

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

source

# swap_cols! — Method.

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

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

source

# scale_row! — Method.

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

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

source

# add_scaled_col! — Method.

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$.

source

# add_scaled_row! — Method.

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$.

source

# transform_row! — Method.

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$.

source

# diagonal — Method.

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

The diagonal elements of $A$ in an array.

source

# reverse_rows! — Method.

reverse_rows!(A::SMat)

Inplace inversion of the rows of $A$.

source

# mod_sym! — Method.

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

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

source

# find_row_starting_with — Method.

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

Tries to find the index $i$ such that $A_{i,p}\ne 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.

source

# reduce — Method.

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.

source

# reduce — Method.

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$.

source

# reduce — Method.

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$.

source

# rand_row — Method.

rand_row(A::SMat) -> SRow

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

source

Changing of the ring:

# map_entries — Method.

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$.

source

# change_base_ring — Method.

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

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

source

Arithmetic

Matrices support the usual operations as well

• +, -, ==

• div, divexact by scalars

• multiplication by scalars

Various products:

# * — Method.

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

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

source

# * — Method.

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

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

source

# * — Method.

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

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

source

# * — Method.

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

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

source

# dot — Method.

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

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

source

# dot — Method.

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

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

source

# dot — Method.

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

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

source

Other:

# sparse — Method.

sparse(A::SMat) -> SparseMatrixCSC

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

source

# ZZMatrix — Method.

ZZMatrix(A::SMat{ZZRingElem})

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

source

# ZZMatrix — Method.

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 $\mathbb{Z}Z$.

source

# Matrix — Method.

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

The same matrix, but as a julia matrix.

source

# Array — Method.

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

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

source