With the release of <strong>ArrayFire v3.4</strong>, 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!

<h4>Sparse Data Structure</h4>
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>.

<h4>Sparse in ArrayFire</h4>
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>.

<h4>Sparse Functionality</h4>
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.

<h4>Example</h4>
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>
#include <arrayfire.h>
#include <iostream>

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 << "Sparsity of A = "
              <<  100.f * (float)sparseGetNNZ(spA) / (float)spA.elements()
              << std::endl;
    std::cout << "Memory Usage of A = " <<  A.bytes() << " bytes" << std::endl;
    std::cout << "Memory Usage of spA = "
              << sparseGetValues(spA).bytes()
               + sparseGetRowIdx(spA).bytes()
               + sparseGetColIdx(spA).bytes()
              << " bytes" << 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 < 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<<"Final difference in solutions:\n";
    af_print(dot(res, res));
}

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

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

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

    return 0;
}
</pre>

<h4>Performance and Memory Footprint</h4>
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!

<h4>Conclusion></h4>
Nothing excites us more than delivering on user requests for features and seeing them solve complex problems in diverse fields. 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.

The <a href="https://github.com/arrayfire/arrayfire">ArrayFire GitHub</a> page is the place for all things ArrayFire. We look forward to your comments, feature requests and anything else that you may want to talk about.