In ArrayFire v3.4, we have added a number of features that will help our vast range of users. Our users have been asking us for a way to reduce the memory footprint of their data on the GPU. And now we have answered our users' requests. Introducing Sparse Matrices in ArrayFire!

[custom_headline type="left" level="h3" looks_like="h3"]Sparse Data Structure[/custom_headline]

A sparse data structure is one where all the non-zero elements are not stored. Sparse matrices are useful when the number of zero-values elements are &gt;&gt; number of non-zero elements (sparsity of the matrix is high).
The sparse data structure in general is stored as 3 arrays.
<ul>
 	<li>A data or values array containing all the non-zero elements</li>
 	<li>A vector for row indices (based on storage format).</li>
 	<li>A vector for column indices (based on storage format).</li>
</ul>
There are many ways to store sparse matrices, the most prominent of which are
<ul>
 	<li><a href="https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_row_.28CSR.2C_CRS_or_Yale_format.29">Compressed Sparse Row</a></li>
 	<li><a href="https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_column_.28CSC_or_CCS.29">Compressed Storage Column</a></li>
 	<li><a href="https://en.wikipedia.org/wiki/Sparse_matrix#Coordinate_list_.28COO.29">Coordinate Format</a>.</li>
</ul>
ArrayFire supports CSR and COO format, with CSC support coming in future releases.

A more detailed description of Sparse Matrix can be found on <a href="https://en.wikipedia.org/wiki/Sparse_matrix">Wikipedia</a>.

[custom_headline type="left" level="h3" looks_like="h3"]Sparse Arrays in ArrayFire[/custom_headline]

In ArrayFire, Sparse arrays use the same user facing object - af::array. This makes working with sparse arrays all the more seamless. We have added a function <a href="http://arrayfire.org/docs/group__method__mat.htm#ga5945f95549bae25f75694a69d98ca584">af::array::issparse()</a> that allows users to check if an array is sparse or not. Although externally, there is no difference between a sparse and dense array, internally, the sparse array is stored as 3 independent arrays along with other meta data. These arrays, as well as other properties can be fetched using the <a href="http://arrayfire.org/docs/group__sparse__func.htm">sparse functions</a>.

As of v3.4, the following features are supported for sparse arrays:
<ul>
 	<li><a href="http://arrayfire.org/docs/group__sparse__func__create.htm">Creation for sparse arrays.</a>
<ul>
 	<li><a href="http://arrayfire.org/docs/group__sparse__func__create.htm#gae4621e58bfaf6d80668309824fde86c0">From 3 (data, row, column) sparse host or device memory vectors.</a></li>
 	<li><a href="http://arrayfire.org/docs/group__sparse__func__create.htm#gae9a843c7e93b05e6b090ca32b271ebe9">From 3 (data, row, column) sparse af::arrays.</a></li>
 	<li><a href="http://arrayfire.org/docs/group__sparse__func__create.htm#gac0d367d8c9daf27b36a6981bbfcdb970">From dense af::array.</a></li>
 	<li>To create a sparse array from host or device dense array, first create a dense af::array, and then pass it to create sparse af::array from dense af::array.</li>
</ul>
</li>
 	<li><a href="http://arrayfire.org/docs/group__sparse__func__convert__to.htm">Storage format conversions.</a>
<ul>
 	<li><a href="http://arrayfire.org/docs/group__sparse__func__create.htm#gac0d367d8c9daf27b36a6981bbfcdb970">Dense to CSR/COO.</a></li>
 	<li><a href="http://arrayfire.org/docs/group__sparse__func__dense.htm">CSR/COO to dense.</a></li>
</ul>
</li>
 	<li><a href="http://arrayfire.org/docs/group__blas__func__matmul.htm#ga63306b6ed967bd1055086db862fe885b">Sparse-Dense Matrix Multiplication (SPMM/SPMV)</a>
<ul>
 	<li>Sparse Matrix - Dense Matrix</li>
 	<li>Sparse Matrix - Dense Vector</li>
 	<li>Options for Sparse Matrix to be transposed.</li>
</ul>
</li>
</ul>
The Sparse-Dense Matrix Multiplcation is a implemented as a part of <a href="http://arrayfire.org/docs/group__blas__func__matmul.htm#ga63306b6ed967bd1055086db862fe885b">af::matmul()</a>. This is part of making the usage of sparse arrays seamless with regular arrays.

As with all ArrayFire features, Sparse support is available across all backends and all architectures.

[custom_headline type="left" level="h3" looks_like="h3"]Example[/custom_headline]

The example we shall use is a conjugate gradient solver. The inputs will be the same except that for the sparse conjugate gradient we will convert the A matrix to sparse.
<pre class="">#include 
#include 

using namespace af;

static size_t dimension = 16 * 1024;
static const int maxIter = 10;
static const int sparsityFactor = 7;

static array A;
static array spA;  // Sparse A
static array x0;
static array b;

void setupInputs()
{
    // Generate a random input: A
    array T = randu(dimension, dimension, f32);
    // Create 0s in input.
    // Anything that is no divisible by sparsityFactor will become 0.
    A = floor(T * 1000);
    A = A * ((A % sparsityFactor) == 0) / 1000;
    // Make it positive definite
    A = transpose(A) + A + A.dims(0)*identity(A.dims(0), A.dims(0), f32);

    // Make A sparse as spA
    spA = sparse(A);

    // Generate x0: Random guess
    x0 = randu(A.dims(0), f32);

    //Generate b
    b = matmul(A, x0);

    std::cout &lt;&lt; "Sparsity of A = "
              &lt;&lt;  100.f * (float)sparseGetNNZ(spA) / (float)spA.elements()
              &lt;&lt; std::endl;
    std::cout &lt;&lt; "Memory Usage of A = " &lt;&lt;  A.bytes() &lt;&lt; " bytes" &lt;&lt; std::endl;
    std::cout &lt;&lt; "Memory Usage of spA = "
              &lt;&lt; sparseGetValues(spA).bytes()
               + sparseGetRowIdx(spA).bytes()
               + sparseGetColIdx(spA).bytes()
              &lt;&lt; " bytes" &lt;&lt; std::endl;
}

void conjugateGradient(const af::array in)
{
    array x = constant(0, b.dims(), f32);
    array r = b - matmul(in, x);
    array p = r;

    for (int i = 0; i &lt; maxIter; ++i) {
        array Ap = matmul(in, p);
        array alpha_num = dot(r, r);
        array alpha_den = dot(p, Ap);
        array alpha = alpha_num/alpha_den;
        r = r - tile(alpha, Ap.dims())*Ap;
        x = x + tile(alpha, Ap.dims())*p;
        array beta_num = dot(r, r);
        array beta = beta_num/alpha_num;
        p = r + tile(beta, p.dims())*p;
    }
    array res = x0 - x;

    std::cout&lt;&lt;"Final difference in solutions:\n";
    af_print(dot(res, res));
}

int main(int argc, char *argv[])
{
    af::info();
    setupInputs();

    std::cout &lt;&lt; "Verifying Dense Conjugate Gradient:" &lt;&lt; std::endl;
    conjugateGradient(A);

    std::cout &lt;&lt; "Verifying Sparse Conjugate Gradient:" &lt;&lt; std::endl;
    conjugateGradient(spA);

    return 0;
}
</pre>
[custom_headline type="left" level="h3" looks_like="h3"]Performance and Memory Footprint[/custom_headline]

Using the CUDA backend on a Tesla K40c, the conjugate gradient function took a time of 64.94ms for dense matrix and 39.31ms for sparse Matrix. That is almost a 2x speedup! Add in the fact that not all computation in the conjugate gradient function is sparse, the speedup for sparse vs dense matrix multiplication is extremely good.

With a <code>sparsityFactor</code> of 7, the sparsity of A is about 26.38%. This means a tremendous amount of space and performance can be achieved by using sparse arrays. For a floating point matrix of size 16M x 16M, a dense matrix uses 1 GB of memory while the sparse array, stored in CSR format, only used 527 MB - halving the memory consumption!
<table>
<tbody>
<tr>
<th>Matrix Elements</th>
<th>Memory Size (Dense)</th>
<th>Memory Size (Sparse)</th>
<th>CG Time (Dense)</th>
<th>CG Time (Sparse)</th>
</tr>
<tr>
<th>16M x 16M</th>
<th>1024 MB</th>
<th>527 MB</th>
<th>64.94ms</th>
<th>39.31ms</th>
</tr>
<tr>
<th>4M x 4M</th>
<th>64 MB</th>
<th>33 MB</th>
<th>4.752ms</th>
<th>3.935ms</th>
</tr>
</tbody>
</table>
Sparse arrays will now let users store more data on the GPU while allowing faster computation, thus allowing to solve for bigger, more complex problem.

[custom_headline type="left" level="h3" looks_like="h3"]Download[/custom_headline]

ArrayFire v3.4 can be downloaded from these locations:
<ul>
 	<li><a href="http://sandbox.arrayfire.com/download/">Binaries</a></li>
 	<li><a href="https://github.com/arrayfire">Source code</a></li>
</ul>
[custom_headline type="left" level="h3" looks_like="h3"]Community[/custom_headline]

ArrayFire is continually improving through the addition of new sparse matrix enhancements. We welcome your feedback:
<ul>
 	<li>General discussion forums on the <a href="https://groups.google.com/forum/#!forum/arrayfire-users">ArrayFire Google Group</a></li>
 	<li>Live discussion chat on the <a href="https://gitter.im/arrayfire/arrayfire?utm_source=badge&amp;utm_medium=badge&amp;utm_campaign=pr-badge&amp;utm_content=badge">ArrayFire Gitter</a></li>
 	<li>Issue reports on the <a href="https://github.com/arrayfire/arrayfire/issues">ArrayFire GitHub</a></li>
</ul>
Finally, as you find success with ArrayFire, we invite you to contribute a post to this blog to share with the broader community. Email <a href="mailto:scott@arrayfire.com">scott@arrayfire.com</a> to contribute to this blog.