# Sparse matrices

Sparse matrix is a matrix populated primarily with zeros. Sparse matrices are stored using special data structures, which allow to store only nonzero elements and save a lot of memory and CPU time when working with such matrices.

Specialized data structures allow to solve computational problems which are impossible to solve when matrices are stored using "conventional", dense storage.

Usually it is linear systems with large number of variables, but almost empty matrix. Such system arise in many areas - from partial differential equations to large scale RBF interpolation problems. All these problems are impossible to solve without sparse matrices. Hence, sparse subpackage of the ALGLIB package contains data structures, which allow to store sparse matrices, and algorithms, which can work with them.

Below we will study different sparse storage formats employed by sparsematrix structure, and functions which can be used on this structure. We also briefly review algorithms which can be applied to sparse matrices.

# Sparse storage formats

Structure used for sparse matrix storage should provide following functionality:

• maintain list of nonzero elements and their values
• retrieve element by its indexes i/j
• enumerate all nonzero elements in the matrix
• insert new element in the matrix
• multiply sparse matrix by vector (or dense matrix)
• other service operations

There are several popular data structures which are used to store sparse matrices, each with its own benefits and drawbacks. ALGLIB package implements two such structures: hash table storage (also known as Dictionary of keys) and CRS (Compressed row storage).

## Hash table storage

Sparse matrix can be stored as hash table which maps keys (row/column indexes) to values of corresponding elements. Such storage format has following features:

• matrix creation - requires to know only matrix size. It is not necessary to know amount of nonzero elements, although you can benefit from such knowledge. It will allow you to allocate required amount of space during matrix creation and to avoid subsequent reallocations of the hash table during initialization of the matrix elements.
• time to insert element - O(1). Elements can be inserted in arbitrary order. However, when hash table becomes nearly full we have to spend O(TableSize) time creating new, larger table and copying elements from the old one.
• memory requirements - O(γ·(2·sizeof(int)+sizeof(double))) per one nonzero element in nearly full table. Here γ is a small constant from [1.2,1.5] which is a consequence of the fact that hash table require some amount of "spare memory".
• time to read element - O(1). Hash tables are quite fast and allow arbitrary access to elements.
• time to enumerate all nonzero elements - O(TableSize), where TableSize is a total size of hash table (including both busy and free elements). Hash table allows to enumerate its entries, but it does not provide any ordering guarantee. Elements of the matrix will be retrieved in arbitrary order.
• linear operations with matrix (matrix-vector products, matrix-matrix products) - not supported. Hash table stores matrix elements in random order, which prevents efficient implementation of the linear algebra operations. Thus, in order to use matrix in linear operations you have to convert it to CRS format.

You may see that hash table storage allows easy creation and modification of sparse matrices. We may access (read and write) elements in arbitrary order, we do not have to know table size in advance, we have fast random access to the table. However, this storage format is not good for and does not allow to use sparse matrix in linear algebra operations like matrix-vector product. You have to convert table to CRS format in order to do so.

Thus, hash table storage can be classified as intermediate representation which is used when you want to initialize sparse matrix.

## CRS representation

CRS (Compressed Row Storage) is another representation which can be used to store sparse matrix. Matrix elements are ordered by rows (from top to bottom) and columns (within row - from left to right). Such storage format has following features:

• matrix creation - requires prior knowledge of both matrix size and nonzero elements count. It is require to tell in advance how many nonzero elements is stored in each row of the matrix.
• time to insert element - O(1). Elements can not be inserted in arbitrary locations. You can insert them only in strict order - within row from left to right, rows are processed from top to bottom.
• memory requirements - O(sizeof(int)+sizeof(double)) per one nonzero element.
• time to read element - O(log(NonzeroElementsInRow)). Retrieval of the arbitrary element requires binary search in the corresponding row.
• time to enumerate all nonzero elements - O(NonzeroCount), where NonzeroCount is a number of nonzero elements. CRS representation guarantees that elements will be retrieved in specific order - rows are enumerated from top to bottom, within row elements are enumerated from left to right.
• linear algebra operations (matrix-vector products, matrix-matrix products) - are supported.

You may see that CRS representation allows to use matrices in the linear algebra operations due to efficient storage format. However, hard initialization is weakness of such representation: a) you have to know in advance number of nonzero elements in the matrix, and b) you have to fill matrix in strict order (from left to right, from top to bottom).

Thus, CRS representation is good for numerical work, but is not good for easy initialization. Luckily, there exists hash table storage format, which is good for initialization, and can be converted to CRS.

## Use cases

Depending on problem being solved we can point out two ways of working with sparse matrix:

• matrix is created using hash table storage format, filled by elements in any convenient order, converted to CRS format. After conversion it is used in linear algebra operations (for example, by iterative solver). Such scenario is good for medium and large-scale problems which fit into RAM. For example, matrix with 10.000.000 nonzero elements will require 500MB of RAM when stored in hash-based format, and about 160MB after conversion to CRS format. We will need more than 660MB in total in order to create hash-based matrix and convert it to CRS, although after conversion only 160MB will be needed.
• in cases when matrix size is close to RAM size, it may be necessary to avoid hash-based storage, i.e. to use CRS format for matrix initialization. Performance will be order of magnitude higher, and memory requirements will be several times lower. However, this way of matrix initialization is hard to implement because each problem has its own natural way to generate matrix, which is often different from initialization order required by CRS format (from left to right, from top to bottom). However, such scenario sometimes allow to work with matrices whose size is close to RAM limit.

Note #1
Although hash table format allows to create matrix without knowing exact number of nonzero elements, it is desirable to pass to constructor at least approximate, order-of-magnitude correct estimate. It will help to avoid costly reallocation of hash table which occurs when table becomes nearly full.

# Operations with sparse matrices

## Creation - hash table storage format

In order to create matrix in a hash table storage format you should use sparsecreate function. This function creates empty (zero) M·N matrix. If you know at least approximate estimate of the nonzero elements count, you can pass it to the constructor function. Such estimate, if given, allows to reserve required amount of memory and avoid reallocation of internal hash table.

After matrix creation you can fill it by elements. sparseset function allows to modify matrix - call for nonzero element rewrites its contents, call for zero element inserts new entry in the hash table. sparseadd function allows to add value to the element - call for nonzero element increments its value, call for zero element inserts new entry in the table. sparse_d_1 example shows how to work with these functions.

After the initialization sparse matrix can be passed to one of the algorithms (say, sparse solver). Most algorithms require sparse matrix to be in the CRS format. In this case you should manually convert it with sparseconverttocrs function before passing elsewhere.

## Creation - CRS format

In order to create matrix in CRS format you should use sparsecreatecrs function. This function accepts following parameters: rows count M, columns count N, array NER which contains number of nonzero elements in each row. Thus, in order to create CRS matrix you must know exact number of nonzero elements in each row of the matrix. Such information often requires a lot of careful programming work, but as result memory requirements become approximately three times lower.

After creation of the CRS matrix you can start filling its elements. Unlike hash table storage, CRS representation requires you to fill elements in strictly specified order: row by row from top to bottom, within row - from left to right. The only function which can be used for initialization is sparseset. You can read sparse_d_2 example which shows how to create matrix in CRS format.

You can read matrix elements by means of two functions. First, you may use sparseget function, which allows to get element value (nonzero or zero) by its index. This function works for matrices in CRS and hash table formats. For hash table storage it has O(1) time complexity, for CRS format it has O(log(NonzeroElementsInRow)) running time (binary search within a row is used). The most important drawback of this function is that you can get value of any specific element (nonzero or zero), but you can't get a list of nonzero elements without scanning all M·N elements, zero or not.

Second option is to use sparseenumerate function. It allows to enumerate all nonzero elements, subsequently returning their indexes and values. It works for both hash table and CRS matrices, and in both cases only O(1) time is needed to get element. However, enumeration order is different for different formats. Hash-based matrices are enumerated in seemingly random order (one which is used by hash table), while CRS format guarantees that elements will be scanned from top row to the bottom one, and within row - from left to right.

## Passing sparse matrix to other algorithms

ALGLIB package contains algorithms which use user-generated sparse matrices - linear solvers, for example. Different algorithms may have different requirements to sparse storage format. Thus, you have to study algorithm description and pass sparse matrix in correct format - one which is required by algorithm.

## Linear algebra operations with sparse matrices

sparse subpackage contains functions for linear algebra operations with sparse matrix - different kinds of matrix-vector or matrix-matrix products are supported. These functions are require sparse matrix to be stored in CRS format because hash-based representation does not allow to perform efficient numerical operations. This subpackage includes following functions:

• sparsemv - calculates A·x for sparse A and vector x.
• sparsemtv - calculates A T·x for sparse A and vector x.
• sparsemv2 - simultaneously calculates A·x and A T·x for square sparse A and vector x.
• sparsesmv - calculates A·x for sparse symmetric A (given by upper or lower triangle) and vector x.
• sparsemm - calculates A·X for sparse A and dense matrix X.
• sparsemtm - calculates A T·X for sparse A and dense matrix X.
• sparsemm2 - simultaneously calculates A·X and A T·X for square sparse A and dense matrix X.
• sparsesmm - calculates A·x for sparse symmetric A (given by upper or lower triangle) and dense matrix X.

## Other functions

In addition to functions listed above, sparse subpackage contains:

• sparsecopy - copy constructor.
• sparseresizematrix - function which allows to reallocate hash-based matrix and free unneeded memory. It can be called for hash-based matrices with low fill factor (ratio of used entries vs. total entries in hash table). Such matrices can be result of long operations involving subsequent modification of matrix entries, with intermediate results being nonzero, but final value of most entries reduced down to zero.
• sparserewriteexisting - this function can be used for thread-safe modification of matrix by changing value of the already existing element. "Already existing" means that element already has nonzero value (thus, there is not need to change internal data structures and no need to perform synchronization). As long as different threads rewrite different elements, such modification can be performed safely and without locks.

ALGLIB Project offers you two editions of ALGLIB:

ALGLIB Free Edition:
offers full set of numerical functionality
extensive algorithmic optimizations
no low level optimizations

ALGLIB Commercial Edition:
flexible pricing
offers full set of numerical functionality
extensive algorithmic optimizations
high performance (SMP, SIMD)

## ALGLIB 3.13.0 for C++

C++ library.
Delivered with sources.
Monolithic design.
Extreme portability.
Editions:   FREE   COMMERCIAL

## ALGLIB 3.13.0 for C#

C# library with native kernels.
Delivered with sources.
VB.NET and IronPython wrappers.
Extreme portability.
Editions:   FREE   COMMERCIAL

## ALGLIB 3.13.0 for Delphi

Delphi wrapper around C core.
Delivered as precompiled binary.
Compatible with FreePascal.
Editions:   FREE   COMMERCIAL

## ALGLIB 3.13.0 for CPython

CPython wrapper around C core.
Delivered as precompiled binary.
Editions:   FREE   COMMERCIAL