22.1.2 Creating Sparse Matrices

There are several means to create sparse matrix.

Returned from a function

There are many functions that directly return sparse matrices. These include speye, sprand, spdiag, etc.

Constructed from matrices or vectors

The function sparse allows a sparse matrix to be constructed from three vectors representing the row, column and data. Alternatively, the function spconvert uses a three column matrix format to allow easy importation of data from elsewhere.

Created and then filled

The function sparse or spalloc can be used to create an empty matrix that is then filled by the user

From a user binary program

The user can directly create the sparse matrix within an oct-file.

There are several basic functions to return specific sparse matrices. For example the sparse identity matrix, is a matrix that is often needed. It therefore has its own function to create it as speye (n) or speye (r, c), which creates an n-by-n or r-by-c sparse identity matrix.

Another typical sparse matrix that is often needed is a random distribution of random elements. The functions sprand and sprandn perform this for uniform and normal random distributions of elements. They have exactly the same calling convention, where sprand (r, c, d), creates an r-by-c sparse matrix with a density of filled elements of d.

Other functions of interest that directly creates a sparse matrices, are spdiag or its generalization spdiags, that can take the definition of the diagonals of the matrix and create the sparse matrix that corresponds to this. For example

     s = diag (sparse(randn(1,n)), -1);

creates a sparse (n+1)-by-(n+1) sparse matrix with a single diagonal defined.

The recommended way for the user to create a sparse matrix, is to create two vectors contain the row and column index of the data and a third vector of the same size containing the data to be stored. For example

      function x = foo (r, j)
       idx = randperm (r);
       x = ([zeros(r-2,1); rand(2,1)])(idx);
      endfunction
     
      ri = [];
      ci = [];
      d = [];
     
      for j=1:c
         dtmp = foo (r, j);  # Returns vector of length r with (:,j) values
         idx = find (dtmp != 0.);
         ri = [ri; idx];
         ci = [ci; j*ones(length(idx),1)];
         d = [d; dtmp(idx)];
      endfor
      s = sparse (ri, ci, d, r, c);

creates an r-by-c sparse matrix with a random distribution of 2 elements per row. The elements of the vectors do not need to be sorted in any particular order as Octave will sort them prior to storing the data. However, pre-sorting the data will make the creation of the sparse matrix faster.

The function spconvert takes a three or four column real matrix. The first two columns represent the row and column index respectively and the third and four columns, the real and imaginary parts of the sparse matrix. The matrix can contain zero elements and the elements can be sorted in any order. Adding zero elements is a convenient way to define the size of the sparse matrix. For example

     s = spconvert ([1 2 3 4; 1 3 4 4; 1 2 3 0]')
     => Compressed Column Sparse (rows=4, cols=4, nnz=3)
           (1 , 1) -> 1
           (2 , 3) -> 2
           (3 , 4) -> 3

An example of creating and filling a matrix might be

     k = 5;
     nz = r * k;
     s = spalloc (r, c, nz)
     for j = 1:c
       idx = randperm (r);
       s (:, j) = [zeros(r - k, 1); rand(k, 1)] (idx);
     endfor

It should be noted, that due to the way that the Octave assignment functions are written that the assignment will reallocate the memory used by the sparse matrix at each iteration of the above loop. Therefore the spalloc function ignores the nz argument and does not preassign the memory for the matrix. Therefore, it is vitally important that code using to above structure should be as vectorized much as possible to minimize the number of assignments and reduce the number of memory allocations.

The above problem can be avoided in oct-files. However, the construction of a sparse matrix from an oct-file is more complex than can be discussed in this brief introduction, and you are referred to section Oct-Files, to have a full description of the techniques involved.