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 RR. 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 {(xi)iRNxi=0 for almost all i}\{ (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_rowMethod
sparse_row(R::Ring, J::Vector{Tuple{Int, T}}) -> SRow{T}

Constructs the sparse row (ai)i(a_i)_i with aij=xja_{i_j} = x_j, where J=(ij,xj)jJ = (i_j, x_j)_j. The elements xix_i must belong to the ring RR.

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

Constructs the sparse row (ai)i(a_i)_i over RR with aij=xja_{i_j} = x_j, where J=(ij,xj)jJ = (i_j, x_j)_j.

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

Constructs the sparse row (ai)i(a_i)_i over RR with aij=xja_{i_j} = x_j, where J=(ij)jJ = (i_j)_j and V=(xj)jV = (x_j)_j.

source

Basic operations

Rows support the usual operations:

  • +, -, ==
  • multiplication by scalars
  • div, divexact
getindexMethod
getindex(A::SRow, j::Int) -> RingElem

Given a sparse row (ai)i(a_i)_{i} and an index jj return aja_j.

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

Returns the row cA+Bc A + B.

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

Returns the row cA+Bc A + B.

source
transform_rowMethod
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)(aA + bB, cA + dB).

source
lengthMethod
length(A::SRow)

Returns the number of nonzero entries of AA.

source

Change of base ring

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

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

source

Maximum, minimum and 2-norm

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

Returns the largest entry of AA.

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

Returns the largest entry of AA.

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

Returns the smallest entry of AA.

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

Returns the ideal AOA \cap O where OO is the maximal order of the coefficient ideals of AA.

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

Returns the smallest entry of AA.

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

Returns AAtA \cdot A^t.

source

Functionality for integral sparse rows

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

Return the sparse row obtained by lifting all entries in AA.

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

Inplace reduction of all entries of AA modulo nn to the positive residue system.

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

Inplace reduction of all entries of AA modulo nn to the symmetric residue system.

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

Inplace reduction of all entries of AA modulo nn to the symmetric residue system.

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

Returns the largest, in absolute value, entry of AA.

source

Conversion to/from julia and AbstractAlgebra types

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

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

source
sparse_rowMethod
sparse_row(A::MatElem)

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

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

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

source

Sparse matrices

Let RR be a commutative ring. Sparse matrices with base ring RR 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 Matn×m(R)\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_matrixMethod
sparse_matrix(R::Ring) -> SMat

Return an empty sparse matrix with base ring RR.

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

Return a sparse nn times mm zero matrix over RR.

source

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

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

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

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

Constructs the sparse matrix over RR corresponding to AA.

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

Constructs the sparse matrix over RR corresponding to AA.

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 AA and BB inplace, that is, the rows of BB are appended to AA.

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

Vertically joins AA and BB.

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

Horizontally concatenates AA and BB, inplace, changing AA.

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

Horizontally concatenates AA and BB.

source

(Normal julia catcat is also supported)

There are special constructors:

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

Return a sparse nn times nn identity matrix over RR.

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

Return a sparse nn times nn zero matrix over RR.

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

Return a sparse nn times mm zero matrix over RR.

source
block_diagonal_matrixMethod
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:

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

Return the submatrix of AA, where the rows correspond to rr and the columns correspond to cc.

source

Transpose:

transposeMethod
transpose(A::SMat) -> SMat

Returns the transpose of AA.

source

Elementary Properties

sparsityMethod
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
densityMethod
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
nnzMethod
nnz(A::SMat) -> Int

Return the number of non-zero entries of AA.

source
isoneMethod
isone(A::SMat)

Tests if AA is an identity matrix.

source
iszeroMethod
iszero(A::SMat)

Tests if AA is a zero matrix.

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

Finds the largest entry of AA.

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

Finds the smallest entry of AA.

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

Finds the largest, in absolute value, entry of AA.

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

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

source
solve_dixon_sfMethod
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 AA of full rank and a sparse matrix (row), find a sparse matrix (row) xx and an integer dd s.th. xA=bdx A = bd holds. The algorithm is a Dixon-based linear p-adic lifting method. If \code{is_int} is given, then dd is assumed to be 11. In this case rational reconstruction is avoided.

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

The square of the product of the norms of the rows of AA.

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

Find a unimodular matrix TT and an upper-triangular EE s.th. TA=ETA = E holds.

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

Reduces gg modulo AA and assumes that AA is upper triangular.

The second return value is the array of pivot elements of AA 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 AA into upper right Hermite normal form.

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

Return the upper right Hermite normal form of AA.

truncatetruncate

indicates if zero rows should be returned or disgarded, when fullhnffull_hnf is set to falsefalse only an upper triangular matrix is computed, ie. the upwards reduction is not performed.

autoauto

if set to true selects a different elemination strategy - here the upwards reduction is temporarily switched off.

limitlimit

marks the last column where size reduction happens, so calling hnfhnf on hcat(A,identitymatrix(SMat,ZZ,nrows(A)))hcat(A, identity_matrix(SMat, ZZ, nrows(A))) with limitlimit set to ncols(A)ncols(A) should produce a non-reduced transformation matrix in the right. If limitlimit is not set, the transformation matrix will also be partly triangular.

source
snfMethod
snf(A::SMat{ZZRingElem})

The Smith normal form (snf) of AA.

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

Given a matrix AA in HNF, extend this to get the HNF of the concatenation with bb.

source
is_diagonalMethod
is_diagonal(A::SMat) -> Bool

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

source
detMethod
det(A::SMat{ZZRingElem})

The determinant of AA using a modular algorithm. Uses the dense (zzModMatrix) determinant on AA for various primes pp.

source
det_mcMethod
det_mc(A::SMat{ZZRingElem})

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

source
valence_mcMethod
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 AA. The valence is the valence of the minimal polynomial ff of transpose(A)Atranspose(A)*A, thus the last non-zero coefficient, typically f(0)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 * AA is computed instead.

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

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

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

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

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

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

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

source

Manipulation/ Access

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

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

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

Given a sparse matrix AA and an index ii, return the ii-th row of AA.

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

Given a sparse matrix AA, a sparse row bb and an index ii, set the ii-th row of AA equal to bb.

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

Swap the ii-th and jj-th row of AA inplace.

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

Swap the ii-th and jj-th column of AA inplace.

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

Multiply the ii-th row of AA by cc inplace.

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

Add cc times the ii-th column to the jj-th column of AA inplace, that is, AjAj+cAiA_j \rightarrow A_j + c \cdot A_i, where (Ai)i(A_i)_i denote the columns of AA.

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

Add cc times the ii-th row to the jj-th row of AA inplace, that is, AjAj+cAiA_j \rightarrow A_j + c \cdot A_i, where (Ai)i(A_i)_i denote the rows of AA.

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

Applies the transformation (Ai,Aj)(aAi+bAj,cAi+dAj)(A_i, A_j) \rightarrow (aA_i + bA_j, cA_i + dA_j) to AA, where (Ai)i(A_i)_i are the rows of AA.

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

The diagonal elements of AA in an array.

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

Inplace reduction of all entries of AA modulo nn to the symmetric residue system.

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

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

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

Given an upper triangular matrix AA over the integers, a sparse row gg and an integer mm, this function reduces gg modulo AA and returns gg modulo mm with respect to the symmetric residue system.

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

Given an upper triangular matrix AA over a field and a sparse row gg, this function reduces gg modulo AA.

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

Given an upper triangular matrix AA over a field and a sparse row gg, this function reduces gg modulo AA.

source
rand_rowMethod
rand_row(A::SMat) -> SRow

Return a random row of the sparse matrix AA.

source

Changing of the ring:

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

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

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

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

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 AbA \cdot b as a dense vector.

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

Return the product AbA \cdot b as a dense matrix.

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

Return the product AbA \cdot b as a dense matrix.

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

Return the product ABA\cdot B as a sparse row.

source
mul_sparseMethod
mul_sparse(A::SMat{T}, B::SMat{T}) -> SMat{T}

Return the product ABA\cdot B as a sparse matrix.

The product of two sparse matrices is, in general, not sparse, so depending on the context, it might be more efficient to use mul_dense(::SMat{T}, ::SMat{T}) where {T} instead.

source
mul_denseMethod
mul_dense(A::SMat{T}, B::SMat{T}) -> MatElem{T}

Return the product ABA\cdot B as a dense matrix.

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

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

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

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

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

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

source

Other:

sparseMethod
sparse(A::SMat) -> SparseMatrixCSC

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

source
ZZMatrixMethod
ZZMatrix(A::SMat{ZZRingElem})

The same matrix AA, but as an ZZMatrix.

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

The same matrix AA, but as an ZZMatrix. Requires a conversion from the base ring of AA to ZZ\mathbb ZZ.

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

The same matrix, but as a julia matrix.

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

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

source