Sparse Matrices
While Julia provides a data structure for sparse matrix computations, the employed design decisions make it difficult to use this implementation for computations over finite fields. This is mainly due to the fact that in the Julia implementation, it is assumed that one can determine the zero and one elements from the data type alone. However, a finite field data type generally also depends on additional parameters, such as the characteristic of the field.
Since the algorithms underlying ConleyDynamics.jl only require basic row and column operations, a specialized sparse matrix implementation is provided in the package. It is briefly described in the following.
Sparse Matrix Format
Sparse matrices in this package have to be of the composite data type SparseMatrix
, which is structured as follows:
ConleyDynamics.SparseMatrix
— TypeSparseMatrix{T}
Composite data type for a sparse matrix with entries of type T
.
The struct has the following fields:
const nrow::Int
: Number of rowsconst ncol::Int
: Number of columnsconst char::Int
: Characteristic of typeT
const zero::T
: Number 0 of typeT
const one::T
: Number 1 of typeT
entries::Vector{Vector{T}}
: Matrix entries corresponding tocolumns
columns::Vector{Vector{Int}}
:column[k]
points to nonzero entries in column krows::Vector{Vector{Int}}
:rows[k]
points to nonzero entries in the k-th row
In this struct, the type T
has to be either Int
or Rational{Int}
, depending on whether the sparse matrix is interpreted as a matrix with entries in the finite field $GF(p)$ for some prime $p$, or over the field of rationals, respectively. The data type has the following fields:
nrow::Int
designates the number of rows.ncol::Int
gives the number of columns.char::Int
specifies the characteristic of the underlying field $F$. Ifchar=0
, then the field is the rationals $\mathbb{Q}$, and one has to haveT = Rational{Int}
. If, on the other hand, the finite field $F = GF(p)$ is used, thenchar=p
has to be a prime number. In this case, the data type of the matrix entries has to beT = Int
.zero::T
provides 0 in the data typeT
.one::T
provides 1 in the data typeT
.columns::Vector{Vector{Int}}
is a vector of integer vectors, which contains the row indices of nonzero matrix entries in each column. More precisely,columns[k]
contains an increasing list of row indices, which give the locations of all nonzero entries in columnk
. Note that the list for each colum has to be strictly increasing.rows::Vector{Vector{Int}}
is a vector of integer vectors, which contains the column indices of nonzero matrix entries in each row. It is the precise dual to the previous field. This time,rows[k]
contains an increasing list of column indices, which correspond to the nonzero entries of the matrix in thek
-th row.entries::Vector{Vector{T}}
is a vector of vectors which contains the actual matrix entries. It is organized in exactly the same way as the fieldcolumns
. In other words, for everyk = 1,..,ncol
the matrix entry in columnk
and rowcolumns[k][j]
is given byentries[k][j]
, wherej
indexes the nonzero column entries from top to bottom.
This data structure is clearly redundant, in the sense that the field rows
is not needed to uniquely determine the matrix. However, the type SparseMatrix
is fundamental for almost every aspect of ConleyDynamics.jl, as it is used to encode the incidence coefficient map $\kappa$, and therefore also the matrix representation of the boundary operator $\partial$. And for many operations on or queries of Lefschetz complexes, one needs fast access to both the cells in the boundary and the coboundary of a given cells. While the boundary can easily be accessed via the field columns
, the fast coboundary access is aided by the field rows
.
We would like to point out that in view of the different underlying fields, sparse matrices should only be manipulated using the specific commands provided by the package. These are described in detail below. If there is a need for additional functionality beyond these first methods, it can be added at a later point in time.
Creating Sparse Matrices
The package provides a number of methods for creating sparse matrices with the data type SparseMatrix
. These are geared towards their usage within ConleyDynamics.jl and are therefore by no means exhaustive:
sparse_from_full
is usually invoked in the formA = sparse_from_full(AF, p=PP)
. The first input argumentAF
has to be a regular Julia integer matrix. This matrix is then converted to sparse format and returned asA
. If the optional parameterp
is omitted, the resulting sparse matrix is over the rational numbers $\mathbb{Q}$, otherwise it is over the finite field with characteristicPP
.full_from_sparse
converts a given sparse matrix into a standard full matrix in Julia. The data type of the entries is eitherRational{Int}
orInt
, depending on whether the sparse input matrix is considered over the rationals $\mathbb{Q}$ or over a finite field, respectively. When invoking this command, be mindful of the size of the sparse matrix!sparse_from_lists
creates a sparse matrix solely based on its nonzero entries and their locations. It expects the following required input arguments, in the order they are listed:nr::Int
: Number of rowsnc::Int
: Number of columnstchar
: Field characteristic, which has to be 0 if $F = \mathbb{Q}$ and a positive prime otherwisetzero::T
: Number 0 of typeT
tone::T
: Number 1 of typeT
r::Vector{Int}
: Vector of row indicesc::Vector{Int}
: Vector of column indicesv::Vector{T}
: Vector of matrix entries
r
,c
, andv
have the same length and that the matrix has entryv[k]
at the location(r[k],c[k])
. Zero entries will be ignored, and multiple entries for the same matrix position raise an error. Furthermore, iftchar>0
, then the entries inv
are all replaced by their values modulotchar
. As mentioned before, iftchar=0
then the entry type has to beT = Rational{Int}
, otherwise we haveT = Int
.lists_from_sparse
takes a sparse matrix and disassembles it into the separate ingredients specified in the discussion of the previous function. In this sense, it is precisely the inverse method ofsparse_from_lists
.sparse_identity
creates a sparse identity matrix. It is invoked asA = sparse_identity(n, p=PP)
, and returns a sparse identity matrixA
withn
rows andn
columns. If the optional characteristic parameter specified and positive, then the matrix is considered over the finite field with characteristicPP
, otherwise it is over the rationals
\mathbb{Q}
`.
Of these methods, the function sparse_from_lists
provides the easiest and quickest way to create a sparse matrix.
Sparse Matrix Access
Access to the entries of sparse matrices is provided via the following commands:
sparse_get_entry
extracts the matrix entryval
of the matrixA
located in rowri
and columnci
, if it is invoked using the commandval = sparse_get_entry(A,ri,ci)
.sparse_set_entry!
sets the matrix entry of the matrixA
located in rowri
and columnci
to the value 'val', if it is invoked using the commandsparse_set_entry!(A,ri,ci,val)
. Internally, this commands makes sure that the above-defined format of the fields of a sparse matrix is preserved. Note that the data type ofval
has to match the type ofA.zero
. Moreover, if the matrix is considered over a finite field the valueval
has to be given as an integer between0
andA.char-1
.sparse_get_column
is invoked asAcol = sparse_get_column(A,ci)
, and it returns the fullci
-th column of the matrixA
as aVector{T}
of lengthA.nrow
.sparse_get_nz_column
returns the row indices for the nonzero entries in theci
-th column of the sparse matrixA
, if invoked asrivec = sparse_get_nz_column(A,ci)
.sparse_minor
creates a minor from a given sparse matrixA
. For this, one needs to specify the row and column indices of the minor in the integer vectorsrvec
andcvec
, respectively, and then invoke the function using the commandAM = sparse_minor(A,rvec,cvec)
. Note that the entries inrvec
andcvec
do not have to be in increasing order, but they are not allowed to contain repeated indices.
One can also read and set sparse matrix values using the overloaded methods y = A[i,j]
and A[i,j] = val
. In the latter case, it is up to the user to make sure that val
respects the underlying sparse matrix field.
Elementary Matrix Operations
The following commands perform the basic sparse matrix operations that are needed for the functionality of the package:
sparse_add_column!
is invoked using the formsparse_add_column!(A,ci1,ci2,cn,cd)
, and it replaces theci1
-th columncolumn[ci1]
ofA
bycolumn[ci1] + (cn/cd) * column[ci2]
. This operation automatically performs the computations over the field $F$ underlying the sparse matrixA
. In other words, if this field is finite, then it determines the inverse of the argumentcd
as part of the computation.sparse_add_row!
is invoked using the formsparse_add_row!(A,ri1,ri2,cn,cd)
, and it replaces theri1
-th rowrow[ri1]
ofA
byrow[ri1] + (cn/cd) * row[ri2]
. As before, this operation automatically performs the computations over the field $F$ underlying the sparse matrixA
.sparse_permute
creates a new sparse matrix by permuting the row and column indices. It is invoked using the commandAP = sparse_permute(A,pr,pc)
, and the integer vectorspr
andpc
have to describe the row and column permutations, respectively.sparse_remove!
is invoked assparse_remove!(A,ri,ci)
and removes the sparse matrix entry in theri
-th row andci
-th colum, i.e., it effectively sets the entry equal to zero.sparse_multiply
computes the matrix product of two sparse matrices. Exceptions are raised if the matrix product is not defined, or if the involved sparse matrices are defined over different fields. One can also use the operator formA*B
to compute the product of sparse matrices.
As mentioned earlier, additional operations can easily be implemented if they become necessary.
Sparse Matrix Information
Finally, ConleyDynamics.jl provides the following functions for quickly extracting certain information from sparse matrices:
sparse_size
is invoked assize = sparse_size(A,dim)
, and it returns the number of rows ifdim=1
, or the number of columns fordim=2
.sparse_low
returns the largest row indexri
of a nonzero entry in theci
-th column of the matrixA
, if used in the formri = sparse_low(A,ci)
. In other words, it returns the row index of the lowest nonzero matrix entry in the column.sparse_is_sut
checks whether a given sparse matrix is strictly upper triangular, and returns eithertrue
orfalse
.sparse_fullness
returns the fullness of a sparse matrix as a floating point number. Here fullness refers to the ratio of the number of nonzero matrix elements and the total number of matrix entries.sparse_sparsity
computes the sparseness of a sparse matrix, which is defined as $1$ minus its fullness, i.e., it is the ratio of the number of zero matrix elements and the total number of matrix entries.sparse_show
can be used to display a sparse matrix in traditional matrix form at the Julia REPL prompt.