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 (RefLink[SparseArray, paclet:ref/SparseArray]).

This tutorial discusses how to create and work with (RefLink[SparseArray, paclet:ref/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 (RefLink[SparseArray, paclet:ref/SparseArray]).

A (RefLink[SparseArray, paclet:ref/SparseArray]) object can be created by giving a list of those elements that are nonzero:

sp=SparseArray[ {{1,1}->4.1,{1,3}->5.2,{3,2}->3.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 (RefLink[MatrixForm, paclet:ref/MatrixForm]):

sp//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:

1.5 sp

%//MatrixForm

Listable operations also work on sparse arrays to thread over all elements:

Sin[sp]

%//MatrixForm

All combinations of matrix multiplication using the function (RefLink[Dot, paclet:ref/Dot]) are supported. This demonstrates the dot product of two sparse arrays:

sp.sp

%//MatrixForm

This demonstrates the dot product of a sparse array with a dense vector:

sp.{1,2,3}

The dot product of a sparse array with a dense matrix is supported:

sp.DiagonalMatrix[ {1,2,3}]

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:

sp

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:

sp+1.5

If the (RefLink[N, paclet:ref/N]) command is applied to a sparse matrix, it works on all the elements. This builds an example 3×3 sparse matrix:

sp=SparseArray[ {{1,1}->4,{1,3}->5,{3,2}->3}];
sp//MatrixForm

When (RefLink[N, paclet:ref/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:

N[sp]//MatrixForm

Here (RefLink[N, paclet:ref/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 (RefLink[N, paclet:ref/N])[0,20] is still 0:

N[sp,20]//MatrixForm

SparseArray

The main function for generating a sparse array is (RefLink[SparseArray, paclet:ref/SparseArray]). This can operate on a matrix to generate its sparse representation:

mat={{1,0,3},{0,3,0},{0,0,1}};
SparseArray[ mat]

%//MatrixForm

(RefLink[SparseArray, paclet:ref/SparseArray]) can also take a list of rules showing the values for certain parts:

sp=SparseArray[ {{1,1}->3.,{2,3}->2.,{2,1}->4.}]

sp//MatrixForm

An equivalent syntax for (RefLink[SparseArray, paclet:ref/SparseArray]) groups the indices and element values into their own lists:

sp1=SparseArray[ {{1,1},{2,3},{2,1}}->{3.,2.,4.}]

The results of the two forms of (RefLink[SparseArray, paclet:ref/SparseArray]) are identical:

sp===sp1

A fuller discussion of the relative advantages of the two forms is given in the section "Rule Inputs for SparseArray".

(RefLink[SparseArray, paclet:ref/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}:

SparseArray[  {{1,1}->3.,{2,3}->2.,{2,1}->4.}, {5,5}]

%//MatrixForm

In this example the default value is set to 1; typically the default value is 0:

SparseArray[  {{1,1}->3.,{2,3}->2.,{2,1}->4.}, {5,5},1.]

%//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:

Sin[%]

%//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 (RefLink[Band, paclet:ref/Band]) if this can be done; this is discussed in the section "Banded Sparse Arrays":

SparseArray[ {{i_,i_}->i,
	{i_,j_}/;i==j+1->j},{10,10}]

%//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 (RefLink[Random, paclet:ref/Random]) to generate sparse arrays with entries that are pseudorandom numbers. You can also use (RefLink[Random, paclet:ref/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:

sp=SparseArray[
	Table[ {Random[Integer,{1,10}],Random[Integer,{1,10}]}->1,{50}],{10,10}];
sp//MatrixForm

If you want to use (RefLink[Random, paclet:ref/Random]) on the right-hand side of the rules in a sparse array, or in any other expression that will evaluate, you should use (RefLink[RuleDelayed, paclet:ref/RuleDelayed]) (entered with :>) to form the rules. If you do not do this, the same number will be used throughout the sparse matrix:

sp=SparseArray[
	 {i_,j_}/;Abs[i-j]<2:>Random[],{4,4}];
sp//MatrixForm

The general principle is that the rules you use with (RefLink[SparseArray, paclet:ref/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 (RefLink[SparseArray, paclet:ref/SparseArray]). These are demonstrated here:

sp1=SparseArray[ {{1,1}->3.,{2,3}->2.,{2,1}->4.}];
sp2=SparseArray[ {{1,1},{2,3},{2,1}}->{3.,2.,4.}];

Both generate identical sparse arrays:

sp1===sp2

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 (RefLink[SparseArray, paclet:ref/SparseArray]) with many rules to mix explicit values with patterns:

sp1=SparseArray[ {{1,2}->3.,{2,3}->2.,{3,1}->4., {i_,i_}->-1}];
sp1//MatrixForm

This demonstrates the form of (RefLink[SparseArray, paclet:ref/SparseArray]) that uses only one rule. The timing is measured to demonstrate performance:

ind=Table[{i,i},{i,50000}];
elems=Range[50000];
SparseArray[ ind->elems]//Timing

It is possible to convert the single rule to multiple rules using (RefLink[Thread, paclet:ref/Thread]). This takes more time than generating the sparse array:

tind=Thread[ind-> elems];//Timing

This uses the multiple rule form of (RefLink[SparseArray, paclet:ref/SparseArray]). It is slower than the single rule form:

SparseArray[tind]//Timing

One reason why the single rule form of (RefLink[SparseArray, paclet:ref/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 (RefLink[SparseArray, paclet:ref/SparseArray]) is slower, as shown:

ind=Developer`FromPackedArray[ind];
elems=Developer`FromPackedArray[elems];
SparseArray[ ind->elems]//Timing

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 (RefLink[Band, paclet:ref/Band]).

This creates a diagonal matrix:

SparseArray[Band[{1,1}]->x,{5,5}]//MatrixForm

This has bands above and below the diagonal:

SparseArray[{Band[{1,2}]->x,Band[{2,1}]->y},{5,5}]//MatrixForm

This inserts a band along the anti-diagonal:

SparseArray[Band[{1,-1},{-1,1},{1,-1}]->x,{5,5}]//MatrixForm

In general, if you can use (RefLink[Band, paclet:ref/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:

sp=IdentityMatrix[4,SparseArray];
sp//MatrixForm

This generates a 4×4 sparse diagonal matrix:

diag={a,b,c,d};
sp=DiagonalMatrix[SparseArray[diag]];
sp//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 (RefLink[IdentityMatrix, paclet:ref/IdentityMatrix]), you can use the (RefLink[WorkingPrecision, paclet:ref/WorkingPrecision]) option:

sp=IdentityMatrix[4,SparseArray,WorkingPrecision->MachinePrecision];
sp//MatrixForm

For a diagonal matrix, you can make sure that the diagonals already are machine-precision numbers:

diag=Range[1.,5];
sp=DiagonalMatrix[SparseArray[diag]];
sp//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 (RefLink[Normal, paclet:ref/Normal]):

sp = SparseArray[ {{1,2}->3,{3,4}->2}]

Normal[sp]

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:

sp=SparseArray[ {{1,2}->3,{3,4}->2},{50000,50000}]

If this is converted to a dense matrix, an exception is thrown because it is not possible to represent this on typical computers:

Normal[sp]

ArrayRules

(RefLink[SparseArray, paclet:ref/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. (RefLink[ArrayRules, paclet:ref/ArrayRules]) generates the rules for a sparse array:

sp = SparseArray[ {{2,1}->5,{1,2}->11},{3,3}]

(RefLink[ArrayRules, paclet:ref/ArrayRules]) returns the rules that represent the sparse array. It returns a rule of blank patterns {_,_} to show the default value:

ArrayRules[ sp]

When a scalar is added to the sparse array, the default value is changed. For example, here the default value is 5:

ArrayRules[sp+5]

(RefLink[ArrayRules, paclet:ref/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:

ArrayRules[ sp,5]

(RefLink[ArrayRules, paclet:ref/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:

Drop[Map[ First, ArrayRules[ sp]],-1]

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 (RefLink[Part, paclet:ref/Part]). Typically (RefLink[Part, paclet:ref/Part]) is entered with [[ ]] notation.

Ways to get pieces of matrices.

Here is a sample sparse matrix.

This gets the third element in the first row:

sp[[1,3]]

This gets the third row; the row is returned as a sparse vector:

sp[[3]]

It can also obtain a column by using (RefLink[All, paclet:ref/All]) to specify all rows; the column is returned as a sparse vector:

sp[[All,3]]

Negative indices are used to refer to the end of the matrix. The following gets the last element of the last row:

sp[[-1,-1]]

You can get ranges of a matrix using ;;. This gets the second through the fourth rows:

sp[[2;;4]]

This gets the second through the fourth columns:

sp[[All,2;;4]]

You can also give a step. This gets every other column:

sp[[All,1;;-1;;2]]

The function (RefLink[Tr, paclet:ref/Tr]) works on the diagonal elements of a matrix. The one-argument form adds them up:

Tr[sp]

(RefLink[Tr, paclet:ref/Tr]) can also take a function as its second argument, which will apply to the diagonal elements. If (RefLink[List, paclet:ref/List]) is used, this returns the diagonal elements:

Tr[sp,List]

Getting Multiple Pieces

It is possible to extract multiple elements by using indices in lists. This is demonstrated with the following sample matrix.

The following gets the first and third elements of the third row:

sp[[3,{1,3}]]

The following gets the first and third rows:

sp[[{1,3}]]

The following gets the first and third elements of the first and second rows:

sp[[{1,2},{1,3}]]

Setting Pieces of Matrices

Setting elements, rows, and columns so that a sparse matrix is updated is quite straightforward by using the Mathematica function (RefLink[Part, paclet:ref/Part]) on the left-hand side of an assignment.

Resetting parts of matrices.

Here is a sample sparse matrix:

sp=SparseArray[ {{1,1}->1,{1,3}->2,{3,2}->3}];
sp//MatrixForm

To update parts of a matrix, you can use (RefLink[Part, paclet:ref/Part]) on the left-hand side of an assignment. This sets the third element of the third row:

sp[[3,3]]=x;
sp//MatrixForm

This sets the second row of the matrix:

sp[[2]]={1,2,3};
sp//MatrixForm

Here, the second column is set:

sp[[All,2]]={-1,-2,-3};
sp//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:

sp[[All,1;;-1;;2]]=z;
sp//MatrixForm

Setting Multiple Pieces

It is possible to set multiple elements by using indices in lists. This is demonstrated with the following sample matrix:

sp=SparseArray[ {{1,1}->1,{1,3}->2,{3,2}->3}];
sp//MatrixForm

The following sets the first and third elements of the second row:

sp[[2,{1,3}]]=x;
sp//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:

sp[[2,{1,3}]]={x,y};
sp//MatrixForm

The following sets the second and third rows:

sp[[{2,3}]]={1,2,3};
sp//MatrixForm

The following gives two different values for the second and third rows:

sp[[{2,3}]]={{1,2,3},{-1,-2,-3}};
sp//MatrixForm

The following sets the first and third elements of the second and third rows:

sp=SparseArray[ {{1,1}->1,{1,3}->2,{3,2}->3}];
sp[[{2,3},{1,3}]]=x;
sp//MatrixForm

The following sets the first and third elements of the second and third rows with different values:

sp[[{2,3},{1,3}]]={{a,b},{c,d}};
sp//MatrixForm

Extracting Submatrices

The range syntax is useful to extract a submatrix.

Extracting submatrices.

sp=SparseArray[ {{1,1}->1,{2,1}->2,{1,3}->2,{3,2}->3,{4,5}->4,{5,4}->2}];
sp//MatrixForm

This extracts the submatrix from (Subscript[m, 2,1]) to (Subscript[m, 3,2]):

sp[[2;;3,1;;2]]//MatrixForm

This extracts the submatrix of rows 1 to 3:

sp[[1;;3]]//MatrixForm

This extracts the submatrix of columns 2 to 4:

sp[[All,2;;4]]//MatrixForm

You can use negative indices to count from the end. This returns the matrix with the first and last columns dropped:

sp[[All,2;;-2]]//MatrixForm

Deleting Rows and Columns

If you want to delete rows or columns, you can use (RefLink[Drop, paclet:ref/Drop]).

Deleting rows and columns.

sp=SparseArray[ {{1,1}->1,{2,1}->2,{1,3}->2,{3,2}->3,{4,5}->4,{5,4}->2}];
sp//MatrixForm

This drops rows 2 through 4:

Drop[sp,{2,4}]//MatrixForm

This drops columns 2 through 4:

Drop[sp,{},{2,4}]//MatrixForm

In this example, rows 2 and 3, and columns 1, 2, and 3 are all dropped:

Drop[sp,{2,3},{1,3}]//MatrixForm

Inserting Rows and Columns

If you want to insert a row, you can use (RefLink[Insert, paclet:ref/Insert]).

Inserting a row.

sp=SparseArray[ {{1,1}->1,{2,1}->2,{1,3}->2,{3,2}->3,{4,5}->4,{5,4}->2}];
sp//MatrixForm

This inserts a row before row 3:

Insert[sp,{x,x,x,x,x},3]//MatrixForm

If you want to insert a column, you must transpose the matrix, insert the column as a row, and then transpose the matrix back. This inserts a column before column 5:

Transpose[Insert[Transpose[sp],{x,x,x,x,x},5]]//MatrixForm

Extending Matrices

You can increase the size of a matrix by padding it with (RefLink[PadLeft, paclet:ref/PadLeft]) and (RefLink[PadRight, paclet:ref/PadRight]).

(RefLink[PadLeft, paclet:ref/PadLeft]) adds the elements to the beginning of the matrix. This example returns a 4×4 sparse matrix:

PadLeft[sp,{4,4}]

You can see the result with (RefLink[MatrixForm, paclet:ref/MatrixForm]):

%//MatrixForm

(RefLink[PadRight, paclet:ref/PadRight]) adds the elements to the end of the matrix; the result is a sparse array:

PadRight[sp,{4,4}]//MatrixForm

You can use (RefLink[PadLeft, paclet:ref/PadLeft]) to pad on the right if you specify negative indices:

PadLeft[sp,{-4,-4}]//MatrixForm

One important use of the padding functions is to replicate and tile a matrix. In this example, the input matrix is extended to have two versions in every row and three in every column.

(RefLink[PadLeft, paclet:ref/PadLeft]) and (RefLink[PadRight, paclet:ref/PadRight]) have a number of extra features. These are described in the Mathematica documentation.

Transpose

Transposition of elements is a general matrix operation:

sp=SparseArray[ {{1,1}->1,{1,3}->2,{2,2}->3}];
sp//MatrixForm

(RefLink[Transpose, paclet:ref/Transpose]) swaps elements at specific indices; the result is a sparse array:

Transpose[sp]

%//MatrixForm

If you wish to compute the conjugate transpose of a sparse matrix, you can use (RefLink[ConjugateTranspose, paclet:ref/ConjugateTranspose]):

sp=SparseArray[{{1,1}->1.,{2,3}->-I,{3,2}->I,{3,3}->10.5}];
spCT=ConjugateTranspose[sp];
spCT//MatrixForm

If a matrix is equal to its conjugate transpose, it is said to be Hermitian:

sp===spCT

You can also test with the function (RefLink[HermitianMatrixQ, paclet:ref/HermitianMatrixQ]):

HermitianMatrixQ[sp]

Rotating Elements

Another structural operation is to rotate elements within an index. This can be done with the functions (RefLink[RotateLeft, paclet:ref/RotateLeft]) and (RefLink[RotateRight, paclet:ref/RotateRight]):

sp=SparseArray[ {{1,1}->1,{1,3}->2,{3,2}->3}];
sp//MatrixForm

This rotates left by one step in the first level, thereby operating on the rows; the result is a sparse array:

RotateLeft[sp,1]

%//MatrixForm

This rotates left by one step in the second level that operates on the columns; the result is a sparse array:

RotateLeft[sp,{0,1}]//MatrixForm

This rotates the columns in the opposite direction:

RotateRight[sp,{0,1}]//MatrixForm

Testing Matrices

Mathematica provides a number of functions for testing sparse matrices and extracting size information.

The predicate (RefLink[MatrixQ, paclet:ref/MatrixQ]) can be used to test sparse matrices. It can also be used to test dense matrices as shown previously.

This tests that sp is a sparse matrix:

MatrixQ[sp]

(RefLink[MatrixQ, paclet:ref/MatrixQ]) takes an optional second argument that specifies a test to apply to every element. In this example, every element is tested to see if it is an integer:

MatrixQ[ sp,IntegerQ]

In this example, every element must be an integer greater than 1. Because some of the elements are 0, the result is (RefLink[False, paclet:ref/False]):

MatrixQ[ sp,IntegerQ[#] && #>1 &]

In addition to (RefLink[MatrixQ, paclet:ref/MatrixQ]), there are predicates (RefLink[VectorQ, paclet:ref/VectorQ]) and (RefLink[ArrayQ, paclet:ref/ArrayQ]), which are useful for testing vectors and tensors:

VectorQ[SparseArray[ {{1}->2,{5}->2,{9}->3}]]

ArrayQ[ SparseArray[ {{1,1,5}->2,{2,1,3}->2,{3,4,2}->3}]]

(RefLink[ArrayQ, paclet:ref/ArrayQ]) can also take a rank argument to test the depth of the array. In this example, the sparse argument is not a rank-4 tensor, and the result is (RefLink[False, paclet:ref/False]):

ArrayQ[ SparseArray[ {{1,1,5}->2,{2,1,3}->2,{3,4,2}->3}],4]

(RefLink[ArrayQ, paclet:ref/ArrayQ]) also takes a third argument that tests each element. In this example, the result is (RefLink[True, paclet:ref/True]) because all the elements are (RefLink[NumberQ, paclet:ref/NumberQ]):

ArrayQ[ SparseArray[ {{1,1,5}->2,{2,1,3}->2,{3,4,2}->3}], 3,NumberQ]

The command (RefLink[Dimensions, paclet:ref/Dimensions]) is useful for extracting size information:

Dimensions[ sp]

(RefLink[Dimensions, paclet:ref/Dimensions]) returns a list of length 2 when the input is a matrix, stating that two indices are used to reference any element in the matrix. Another way to test the number of indices required to reference elements is with (RefLink[ArrayDepth, paclet:ref/ArrayDepth]); this is equivalent to (RefLink[Length, paclet:ref/Length])[(RefLink[Dimensions, paclet:ref/Dimensions])[sp]]:

ArrayDepth[ sp]

To compare if the elements of two matrices are equal, it is possible to use (RefLink[Equal, paclet:ref/Equal]), typically entered using a == shorthand notation. For example, comparing a matrix with itself returns (RefLink[True, paclet:ref/True]):

sp=SparseArray[ {{1}->2,{5}->2,{9}->3}];
sp==sp

(RefLink[Equal, paclet:ref/Equal]) uses the value of numbers, so it can be used to compare integer and real values:

sp==N[sp]

It should be noted that (RefLink[Equal, paclet:ref/Equal]) works on any Mathematica expression. If you want to compare two matrices for equality using properties of the matrix as a whole, it may be better to compare matrix norms. These are discussed under "Matrix Norms".

Further Structural Operations

This section discusses some further structural operations that are useful for working with matrices.

This generates a sample matrix for demonstration:

sp=SparseArray[ {{1,1}->a,{2,3}->b,{3,3}->c}];
sp//MatrixForm

This flattens the matrix into a vector:

vec=Flatten[sp];
vec//Normal

The vector can be partitioned back into a matrix with rows of length 3 using (RefLink[Partition, paclet:ref/Partition]):

Partition[vec,3]

Matrices can be joined together with the function (RefLink[Join, paclet:ref/Join]):

Join[ sp,sp]//MatrixForm

Alternatively, you can join the new matrix as new columns:

Join[ sp,sp,2]//MatrixForm

A new row can be inserted at the end of a matrix with (RefLink[Append, paclet:ref/Append]):

Append[sp,{x,x,x}]//MatrixForm

It should be noted that this can also be done with (RefLink[Insert, paclet:ref/Insert]). This is shown in the section "Inserting Rows and Columns".

Element-wise Operations

If you want to work on the elements of a sparse matrix, you can do this easily with Mathematica. First, a matrix of floating-point numbers is created:

sp=SparseArray[ {{1,1}->1.2,{1,3}->2.3,{3,2}->3.1}];
sp//MatrixForm

In general, arithmetic operations applied to a matrix work on each element. This adds 5 to the matrix:

sp+5//MatrixForm

Here, every element of the matrix is squared:

(sp^2)//MatrixForm

If one matrix is divided by another, the division is done element by element. If the dimensions of the two matrices do not agree, then an error ensues:

sp/(1+sp)//MatrixForm

To apply the (RefLink[Sin, paclet:ref/Sin]) function to every element, you apply (RefLink[Sin, paclet:ref/Sin]) to the entire matrix:

Sin[sp]//MatrixForm

If both of the arguments of an operation are matrices, the operation is carried out on corresponding elements. This result is the same as multiplying the matrix by two:

sp +sp//MatrixForm

If one argument is a matrix and the other is a vector, the operation is carried out between rows of the matrix and elements of the vector. This is an efficient way of generating a diagonal sparse matrix:

sp=SparseArray[ {i_,i_}->1,{5,5}];
vec=SparseArray[Range[5]];
sp//MatrixForm

sp vec

%//MatrixForm

These operations are all fast because they only need to operate on the elements that are actually stored (including the default element).

Note that multiplying two matrices using the operator (RefLink[Times, paclet:ref/Times]) is equivalent to multiplying the corresponding elements. On the other hand, matrix multiplication can be done with the function (RefLink[Dot, paclet:ref/Dot]), which is described in the section "Matrix Multiplication". Element-wise multiplication is shown in the following example:

sp = SparseArray[{{1,2}->3, {2,3}->4}];
sp sp//MatrixForm

Listability

If you want to apply your own function to each element in a matrix, you may give your function the attribute (RefLink[Listable, paclet:ref/Listable]). This function squares each element and divides the result by 3:

fun[ sp]//MatrixForm

When a listable operation is applied to a sparse array, it only operates on the elements that are actually stored (including the default element), so the value of fun[0] is only computed once.

Map

Instead of using listability, you can use (RefLink[Map, paclet:ref/Map]) to apply a function to every element in a matrix.

(RefLink[Map, paclet:ref/Map]) will apply the function fun to every element in the matrix:

Map[f,sp,{-1}]

You can see the matrix that the sparse array represents. Note that because the function f did not have any definition, it is just kept wrapped around every element:

%//MatrixForm

Here, a function that squares its argument and divides the result by 5 is applied to every element in the matrix:

Map[#^2/5 &,sp,{-1}]

%//MatrixForm

When (RefLink[Map, paclet:ref/Map]) applies a function to a sparse array, it only operates on the elements that are actually stored (including the default element).

Visualization of Sparse Matrices

This section reviews the functions that are available for formatting and plotting sparse matrices. Because sparse matrices are well integrated into the system, most of the examples in this section are very similar to the way that dense matrices work. Visualization techniques for dense matrices are described under "Visualization of Matrices".

Formatting Sparse Matrices

Sparse matrices can be formatted with the function (RefLink[MatrixForm, paclet:ref/MatrixForm]):

sp = SparseArray[ {{1,2}->4,{4,2}->2,{2,3}->1}];
MatrixForm[sp]

(RefLink[MatrixForm, paclet:ref/MatrixForm]) also works for vectors and higher-rank sparse arrays; the braces can help in understanding the grouping:

MatrixForm[ SparseArray[ {{1,2,2,1}->4,{4,1,2,2}->2,{2,3,3,1}->1}]]

The view that you get from (RefLink[MatrixForm, paclet:ref/MatrixForm]) is dense, which can help you see the sparsity pattern. However, it can quickly lead to very large output, especially as the rank of the array increases.

Plotting Sparse Matrices

A convenient way to plot sparse matrices is with the function (RefLink[MatrixPlot, paclet:ref/MatrixPlot]), which is described in greater detail under "Plotting Matrices".

This reads in a sparse matrix, using the (RefLink[Import, paclet:ref/Import]) function demonstrated in the previous section:

sp=Import["LinearAlgebraExamples/Data/dwg961b.cua","HarwellBoeing"]

The matrix can be plotted. This is a useful way to see the structure of the matrix:

MatrixPlot[sp]

If the matrix is multiplied with itself, the result is less sparse, but it still has a related structure:

MatrixPlot[sp.sp]

Import and Export of Sparse Matrices

Mathematica provides a number of different tools for I/O. You can save your data in a file so that you or a colleague can continue to work with it in Mathematica later. For this you might want to use some of the expression I/O functions; these are discussed under "Expression Input and Output".

If you want to work with matrices from a source external to Mathematica using specific data formats, the functions (RefLink[Import, paclet:ref/Import]) and (RefLink[Export, paclet:ref/Export]) are useful. The (RefLink[Import, paclet:ref/Import]) function supports a variety of different formats, some of which are relevant to sparse matrices.

There are two formats that are specific to sparse matrices and are supported in Mathematica: "Harwell–Boeing" and "Matrix Market". This example imports the matrix in the file dwg961b.cua using the Harwell–Boeing format:

sp=Import["LinearAlgebraExamples/Data/dwg961b.cua","HarwellBoeing"]

The result is a 961×961 matrix of numbers with 10591 nonzero elements:

MatrixQ[ sp,NumberQ]

You can confirm the dimensions of the matrix:

Dimensions[sp]

If the matrix is multiplied with itself, the result may not be as sparse as before:

sp.sp

Many examples of matrices in Harwell–Boeing and Matrix Market formats can be found at http://math.nist.gov/MatrixMarket/index.html.

Matrix Multiplication

Matrix multiplication is carried out in Mathematica with the function (RefLink[Dot, paclet:ref/Dot]), which is typically entered with a dot shorthand syntax. The operation is fully supported for sparse arrays:

sp = SparseArray[ {{1,2}->4,{4,2}->2,{2,3}->1,{2,4}->1}];
MatrixForm[sp]

The sparse matrix can be multiplied by itself; the result is also a sparse matrix:

sp.sp

%//MatrixForm

This multiplies a sparse matrix by a dense vector. In this case the result is dense:

sp.{4,5,2,1}

If the vector is made into a sparse array, the result of the multiplication will be sparse:

sp.SparseArray[ {4,5,1,1}]

%//MatrixForm

In general, the result of matrix multiplication will be sparse if both arguments are sparse. This is not true if the default values are different:

(sp+4).sp

What is particularly important about sparse matrix vector multiplication is that it is fast. In this example, a 100000×100000 tri-diagonal sparse matrix is multiplied by a dense vector:

sp=SparseArray[{Band[{1,1}],Band[{1,2}],Band[{2,1}]}->1.,{100000,100000},0.];
vec=Table[Random[],{100000}];

sp.vec;//Timing

Another issue with multiplication of sparse and dense matrices is that there are differences in the speed of computation depending on whether the dense matrix is on the left or the right:

mat=Table[Random[],{1000},{1000}];
sp=SparseArray[{Band[{1,1}]->1.},{1000,1000}];

The first example is slower than the second, even though the result is the same in both cases:

Timing[Head[mat.sp]]

Timing[Head[sp.mat]]

One possibility to improve the speed is to convert the dense matrix to a sparse representation. Now the result is a sparse matrix:

sp1=SparseArray[mat];
Timing[Head[sp1.sp]]

Another solution would be to reverse the order by working with the transpose of each matrix:

Timing[Head[Transpose[ Transpose[sp].Transpose[mat]]]]

Outer Product

The outer product is a way to build a higher-rank tensor from ones of lower rank. Mathematica provides this functionality with the function (RefLink[Outer, paclet:ref/Outer]). One use of this is to combine two vectors to form a matrix as an outer product. This works for sparse vector input to generate a sparse array:

vec1=SparseArray[{a,0,c}];
vec2=SparseArray[{0,y,z}];
Outer[Times,vec1,vec2]

%//MatrixForm

This generation of sparse outer products can be very advantageous if you want to use Mathematica sparse arrays as general sparse data structures.

Matrix Permutations

Many matrix techniques rely on ordering a matrix in particular ways. For example, some techniques try to order the matrix to put elements on the diagonal, while others try to group certain elements into dense blocks. The Mathematica function (RefLink[Part, paclet:ref/Part]) is very well suited to applying permutations to the rows and columns of a matrix.

Applying permutations to matrices.

This generates a random matrix:

sp=SparseArray[Table[Random[Integer,{0,3}],{4},{4}]];
sp//MatrixForm

Now the matrix will be reordered so that the rows with smallest 2-norm are first. ("Norms" are discussed under "Matrix Computations".) First, the norm of each row is computed:

Map[ Norm, sp]

This computes the permutation that is necessary to reorder the matrix:

perm=Ordering[%,All,Less]

This applies the ordering to the rows of the matrix; the result has rows with smallest 2-norms first:

spperm=sp[[perm]];
spperm//MatrixForm

Now the inverse permutation is applied using part assignment. Note that this modifies the matrix held by the symbol pmat. Part assignment is described previously:

spperm[[perm]]=spperm;
spperm//MatrixForm

Converting Equations to Sparse Arrays

Matrices are so important in many areas of science and technology because they are an efficient way to represent linear systems of equations. Mathematica is unique among technical computing systems in that it combines very efficient ways to work with matrices and also with the equations that the matrices represent. It is easy to go from a matrix to a system of equations and back. Here, a sparse matrix is multiplied by a vector of the unknowns and a system of equations is formed.

sp=SparseArray[ { {1,1}->3.1,{1,3}->3.1,{2,3}->1.1,{3,1}->3.2,{3,2}->3.2,{2,1}->1.9}];
eqns=Thread[sp.{x,y,z}=={2,6,8}]

sp//MatrixForm

These equations can be solved with the algebraic equation solver, (RefLink[Solve, paclet:ref/Solve]):

Solve[eqns,{x,y,z}]

However, it is a bit more involved to go from the equation to the matrix representation. Mathematica provides the function (RefLink[CoefficientArrays, paclet:ref/CoefficientArrays]) to make this transformation easier. This takes the equations that were generated; the result is a list of a sparse vector and a sparse matrix:

CoefficientArrays[eqns,{x,y,z}]

You can see the details of the sparse arrays with (RefLink[MatrixForm, paclet:ref/MatrixForm]):

Map[MatrixForm,%]

(RefLink[CoefficientArrays, paclet:ref/CoefficientArrays]) is general and will work for nonlinear as well as linear polynomial equations:

{s0,s1,s2,s3} = CoefficientArrays[ {3+x^2+5y^3,1+x^3+2y^2}, {x,y}]

When the input is a nonlinear polynomial, the result will include higher-rank tensors:

Map[MatrixForm,%]

You can still regain the original expressions with dot products:

s0+s1.{x,y}+s2.{x,y}.{x,y}+s3.{x,y}.{x,y}.{x,y}

The ability to work directly with systems of linear equations can be very advantageous for certain applications. For example, generating finite difference solutions. This is demonstrated in the example "Finite Difference Solutions".

SparseArray Data Format

There are several different formats that can be used to hold sparse arrays. Each has its own advantages and disadvantages. Mathematica uses the Compressed Sparse Row (CSR) format as an internal storage format. This can be explained for a sample matrix shown here.

sp=SparseArray[{{1,1}->a11,{1,3}->a13,{2,1}->a21,
{2,4}->a24,{3,2}->a32,{3,3}->a33},{3,4}];
sp//MatrixForm

You can see the internal form in the (RefLink[InputForm, paclet:ref/InputForm]) of the matrix:

sp//InputForm

This format is general enough to describe arbitrary rank tensors. Other advantages of the format include the fact that data elements in the same row are stored next to each other, which leads to better cache performance. A disadvantage is that it is not optimized for inserting new elements into the matrix.