Working with Sparse Arrays
Sparse representations of matrices are useful because they do not store every element. If one particular value appears very frequently it can be very advantageous to use a sparse representation.
Mathematica offers a sparse representation for matrices, vectors, and tensors with
SparseArray.
This tutorial discusses how to create and work with
SparseArray objects in
Mathematica. If you are interested in carrying out linear algebra computations on sparse matrices, you should consult "
Matrix Computations".
Basic Operations
The basic object for representing a sparse matrix in
Mathematica is a
SparseArray.
| SparseArray[list] | a SparseArray version of an ordinary list |
| SparseArray[{{i1,j1}->v1,{i2,j2}->v2,...},{m,n}] |
| an m×n sparse array with element {ik, jk} having value vk |
| SparseArray[{{i1,j1},{i2,j2},...}->{v1,v2,...},{m,n}] |
| the same sparse array |
| SparseArray[data,{m,n},def] | an m×n sparse array with default element def |
| SparseArray[Band[b]->v,{m,n}] | an m×n banded sparse array |
| Normal[array] | the ordinary list corresponding to a SparseArray |
| ArrayRules[m] | positions of nonzero elements |
A SparseArray object can be created by giving a list of those elements that are nonzero.
| Out[1]= |  |
|
By default, a sparse matrix will print with a special output format. You can see the matrix that the sparse array represents by using MatrixForm.
Out[2]//MatrixForm= |
| |  |
|
Operations on sparse matrices are all equivalent to the operations on dense matrices. For example, arithmetic is supported and a sparse array is the result.
| Out[3]= |  |
Out[4]//MatrixForm= |
| |  |
|
Listable operations also work on sparse arrays to thread over all elements.
| Out[5]= |  |
Out[6]//MatrixForm= |
| |  |
|
All combinations of matrix multiplication using the function Dot are supported. This demonstrates the dot product of two sparse arrays.
| Out[7]= |  |
Out[8]//MatrixForm= |
| |  |
|
This demonstrates the dot product of a sparse array with a dense vector.
| Out[9]= |  |
|
The dot product of a sparse array with a dense matrix is supported.
| Out[10]= |  |
|
Sparse representations are useful because they do not store every element. If one particular value, typically this is zero, appears many times in the sparse array, it can be much more efficient if only elements that are different from this common value are stored. The default output format shows the number of nondefault elements and the dimensions.
| Out[11]= |  |
|
If a single number is added to the sparse array, it is added to all elements and also to the default element, which was zero. Now that the default element is no longer zero but 1.5, it is shown in the output.
| Out[12]= |  |
|
If the N command is applied to a sparse matrix, it works on all the elements. This builds an example 3×3 sparse matrix.
Out[14]//MatrixForm= |
| |  |
|
When N is applied to the sparse matrix, the result is a sparse matrix with elements (including the default element) that are all approximate machine numbers.
Out[15]//MatrixForm= |
| |  |
|
Here N with a precision argument is applied to the matrix. This generates a sparse matrix of approximate real numbers with 20 digits of precision. Note that N[0, 20] is still 0.
Out[16]//MatrixForm= |
| |  |
|
SparseArray
The main function for generating a sparse array is SparseArray. This can operate on a matrix to generate its sparse representation.
| Out[2]= |  |
Out[3]//MatrixForm= |
| |  |
|
SparseArray can also take a list of rules showing the values for certain parts.
| Out[4]= |  |
Out[5]//MatrixForm= |
| |  |
|
An equivalent syntax for SparseArray groups the indices and element values into their own lists.
| Out[6]= |  |
|
The results of the two forms of SparseArray are identical.
| Out[7]= |  |
|
A fuller discussion of the relative advantages of the two forms is given in the section "
Rule Inputs for SparseArray".
SparseArray can also accept the dimensions of the matrix to be created. Here a 5×5 matrix is created even though the maximum explicit index was {2,3}.
| Out[8]= |  |
Out[9]//MatrixForm= |
| |  |
|
In this example the default value is set to 1; typically the default value is 0.
| Out[10]= |  |
Out[11]//MatrixForm= |
| |  |
|
Allowing the default value to be changed makes many element-wise operations very fast. They just need to work on the elements that are actually present and the default value.
| Out[12]= |  |
Out[13]//MatrixForm= |
| |  |
|
You can also give patterns for the rules. This can often be a convenient way to build structured matrices. You can use the names of the patterns on the right-hand side of the rules. Typically, it is better to use Band if this can be done, this is discussed in the section " Banded Sparse Arrays".
| Out[14]= |  |
Out[15]//MatrixForm= |
| |  |
|
Mathematica uses symbolic algebraic techniques to simplify certain of the patterns for building sparse arrays. This allows it to construct the sparse array without testing every potential element to see if it is actually present in the array, thus providing a significant saving in computational time.
You can use the Mathematica function Random to generate sparse arrays with entries that are pseudorandom numbers. You can also use Random to generate the indices for the sparse array. In this example a 10×10 sparse matrix with at most 50 nondefault entries is generated.
Out[17]//MatrixForm= |
| |  |
|
If you want to use Random on the right-hand side of the rules in sparse array, or any other expression that will evaluate, you should use RuleDelayed (entered with ⧴) to form the rules. If you do not do this, the same number will be used throughout the sparse matrix.
Out[19]//MatrixForm= |
| |  |
|
The general principle is that the rules you use with
SparseArray work in the typical way for rules in
Mathematica.
Rule Inputs for SparseArray
There are two different ways that rules can be used as input syntax for SparseArray. These are demonstrated here. |
Both generate identical sparse arrays.
| Out[3]= |  |
|
The first form, which has many rules, is convenient if you want to mix explicit indices with patterns. It is also more readable for small examples. The second is more efficient, and is preferred if you only have explicit indices, for example, after reading data from a file.
This uses SparseArray with many rules to mix explicit values with patterns.
Out[5]//MatrixForm= |
| |  |
|
This demonstrates the form of SparseArray that uses only one rule. The timing is measured to demonstrate performance.
| Out[8]= |  |
|
It is possible to convert the single rule to multiple rules using Thread. This takes more time than generating the sparse array.
| Out[9]= |  |
|
This uses the multiple rule form of SparseArray. It is slower than the single rule form.
| Out[10]= |  |
|
One reason why the single rule form of SparseArray is faster is that the indices and elements can use packed arrays, an efficient storage technology. If they are converted from packed arrays then SparseArray is slower as shown.
| Out[13]= |  |
|
More information on packed arrays is found under "
Packed Arrays".
Banded Sparse Arrays
If you want to build sparse matrices that have some type of banded structure, this can be achieved with
Band.
| SparseArray[Band[b]->v,{m,n}] | an m×n banded sparse array |
| Band[{i,j}] | a diagonal band that starts with the position {i, j} |
| Band[{imin,jmin},{imax,imax}] | a diagonal band from {imin, jmin} to {imax, imax} |
| Band[{imin,jmin},{imax,imax}{di,dj}] | a diagonal band from {imin, jmin} moving with step {di, dj} |
This creates a diagonal matrix.
Out[1]//MatrixForm= |
| |  |
|
This has bands above and below the diagonal.
Out[2]//MatrixForm= |
| |  |
|
This inserts a band along the anti-diagonal.
Out[6]//MatrixForm= |
| |  |
|
In general, if you can use
Band to create a sparse array it will be more efficient.
Identity and Diagonal Sparse Matrices
There are a number of ways to create identity and diagonal sparse matrices in
Mathematica.
This generates a 4×4 sparse identity matrix.
Out[2]//MatrixForm= |
| |  |
|
This generates a 4×4 sparse diagonal matrix.
Out[5]//MatrixForm= |
| |  |
|
If you want to compute with floating-point numbers it can be advantageous to use matrices that contain floating-point entries; this is described in more detail under " Matrix Contents". For IdentityMatrix you can use the WorkingPrecision option.
Out[7]//MatrixForm= |
| |  |
|
For a diagonal matrix you can make sure that the diagonals already are machine precision numbers.
Out[10]//MatrixForm= |
| |  |
|
Notice how all the entries of the matrix are machine-precision numbers. This can help to improve the efficiency of your computations. The different types of matrices that
Mathematica can work with are described in more detail under "
Matrix Types".
Normal
To convert from a sparse array to the dense object that it represents, you can use Normal.
| Out[1]= |  |
| Out[2]= |  |
|
Of course, it is very straightforward to make a sparse array that cannot be represented on your system. For example, the following sparse matrix only has two nonzero elements, but it has 2499999998 zero elements.
| Out[3]= |  |
|
If this is converted to a dense matrix, an exception is thrown because it is not possible to represent this on typical computers.
| Out[4]= |  |
|
ArrayRules
SparseArray can accept a list of rules to form a sparse array. These rules hold the indices and values of nonzero elements. In the following example, the element at position {2,1} has the value 5. ArrayRules generates the rules for a sparse array.
| Out[1]= |  |
|
ArrayRules returns the rules that represent the sparse array. It returns a rule of blank patterns {_, _} to show the default value.
| Out[2]= |  |
|
When a scalar is added to the sparse array, the default value is changed. For example, here the default value is 5.
| Out[3]= |  |
|
ArrayRules can take a second argument, which is the default value shown on output. Here, a default rule with value 5 is used. Because there is only one element with this value, the list of rules is much longer.
| Out[4]= |  |
|
ArrayRules is a useful way to get information about a sparse array such as the nondefault elements or the default value. This example shows how to get the indices of the nondefault elements.
| Out[5]= |  |
|
Structural Operations
Structural operations on sparse arrays are all equivalent to the operations on dense matrices.
Getting Pieces of Matrices
Extracting elements, rows, and columns of a sparse matrix is quite straightforward with the
Mathematica function
Part. Typically
Part is entered with
[[ ]] notation.
| m[[i,j]] | the i, jth entry |
| m[[i]] | the ith row |
| m[[i;;i]] | rows i through j |
| m[[All,i]] | the ith column |
| m[[All,i;;j]] | columns i through j |
| m[[{i1,...,ir},{j1,...,js}]] | the r×s submatrix with elements having row indices ik and column indices jk |
| Tr[m,List] | list of the diagonal elements of m |
Ways to get pieces of matrices.
Here is a sample sparse matrix.
Out[2]//MatrixForm= |
| |  |
|
This gets the third element in the first row.
| Out[3]= |  |
|
This gets the third row; the row is returned as a sparse vector.
| Out[4]= |  |
|
It can also obtain a column by using All to specify all rows; the column is returned as a sparse vector.
| Out[5]= |  |
|
Negative indices are used to refer to the end of the matrix. The following gets the last element of the last row.
| Out[6]= |  |
|
You can get ranges of a matrix using ;;. This gets the second through the fourth rows.
| Out[7]= |  |
|
This gets the second through the fourth columns.
| Out[8]= |  |
|
You can also give a step. This gets every other column.
| Out[9]= |  |
|
The function Tr works on the diagonal elements of a matrix. The one-argument form adds them up.
| Out[7]= |  |
|
Tr can also take a function as its second argument, which will apply to the diagonal elements. If List is used, this returns the diagonal elements.
| Out[8]= |  |
|
Getting Multiple Pieces
It is possible to extract multiple elements by using indices in lists. This is demonstrated with the following sample matrix.
Out[2]//MatrixForm= |
| |  |
|
The following gets the first and third elements of the third row.
| Out[3]= |  |
|
The following gets the first and third rows.
| Out[4]= |  |
|
The following gets the first and third elements of the first and second rows.
| Out[5]= |  |
|
Setting Pieces of Matrices
Setting elements, rows, and columns so that a sparse matrix is updated is quite straightforward by using the
Mathematica function
Part on the left-hand side of an assignment.
| m={{a11,a12,...},{a21,a22,...},...} | assign m to be a matrix |
| m[[i,j]]=v | reset element {i, j} to be v |
| m[[i]]=v | reset all elements in row i to be v |
| m[[i]]={v1,v2,...} | reset elements in row i to be {v1, v2, ...} |
| m[[All,j]]=v | reset all elements in column j to be v |
| m[[All,j]]={v1,v2,...} | reset elements in column j to be {v1, v2, ...} |
Resetting parts of matrices.
Here is a sample sparse matrix.
Out[2]//MatrixForm= |
| |  |
|
To update parts of a matrix you can use Part on the left-hand side of an assignment. This sets the third element of the third row.
Out[4]//MatrixForm= |
| |  |
|
This sets the second row of the matrix.
Out[6]//MatrixForm= |
| |  |
|
Here, the second column is set.
Out[8]//MatrixForm= |
| |  |
|
You can also use the range syntax for setting pieces of the matrix. This sets every element in every other row to be z.
Out[10]//MatrixForm= |
| |  |
|
Setting Multiple Pieces
It is possible to set multiple elements by using indices in lists. This is demonstrated with the following sample matrix.
Out[2]//MatrixForm= |
| |  |
|
The following sets the first and third elements of the second row.
Out[4]//MatrixForm= |
| |  |
|
If the right-hand side of the assignment is a list that matches the number of elements being assigned, the assignment is done element by element. Thus, the following gives two different values for the first and third elements of the second row.
Out[6]//MatrixForm= |
| |  |
|
The following sets the second and third rows.
Out[8]//MatrixForm= |
| |  |
|
The following gives two different values for the second and third rows.
Out[10]//MatrixForm= |
| |  |
|
The following sets the first and third elements of the second and third rows.
Out[13]//MatrixForm= |
| |  |
|
The following sets the first and third elements of the second and third rows with different values.
Out[15]//MatrixForm= |
| |  |
|
Extracting Submatrices
The range syntax is useful to extract a submatrix.
| m[[i0;;i1,j0;;j1]] | extract the submatrix with rows i0 through i1 and columns j0 through j1 |
| m[[i0;;i1]] | extract the submatrix with rows i0 through i1 |
| m[[All,j0;;j1]] | extract the submatrix with columns j0 through j1 |
Extracting submatrices.
Out[2]//MatrixForm= |
| |  |
This extracts the submatrix from m2,1 to m3,2.
Out[3]//MatrixForm= |
| |  |
|
This extracts the submatrix of rows 1 to 3.
Out[4]//MatrixForm= |
| |  |
|
This extracts the submatrix of columns 2 to 4.
Out[5]//MatrixForm= |
| |  |
|
You can use negative indices to count from the end. This returns the matrix with the first and last columns dropped.
Out[6]//MatrixForm= |
| |  |
|
Deleting Rows and Columns
If you want to delete rows or columns, you can use
Drop.
| Drop[m,{i0,i1}] | delete rows i0 through i1 |
| Drop[m,{},{j0,j1}] | delete columns j0 through |