Sparse Linear Algebra

Sparse Linear Algebra

ADCME augments TensorFlow APIs by adding sparse linear algebra support. In ADCME, sparse matrices are represented by SparseTensor. This data structure stores indices, rows and cols of the sparse matrices and keep track of relevant information such as whether it is diagonal for performance consideration. The default is row major (due to TensorFlow backend).

When evaluating SparseTensor, the output will be SparseMatrixCSC, the native Julia representation of sparse matrices

A = run(sess, s) # A has type SparseMatrixCSC{Float64,Int64}

Sparse Matrix Construction

ii = [1;2;3;4]
jj = [1;2;3;4]
vv = [1.0;1.0;1.0;1.0]
s = SparseTensor(ii, jj, vv, 4, 4)
using SparseArrays
s = SparseTensor(sprand(10,10,0.3))
D = Array(sprand(10,10,0.3)) # a dense array
d = constant(D)
s = dense_to_sparse(d)

There are also special constructors.

Diagonal matrix with diagonal vspdiag(v)
Empty matrix with size m, nspzero(m, n)
Identity matrix with size mspdiag(m)

Matrix Traits

  1. Size of the matrices

    size(s) # (10,20)
    size(s,1) # 10
  2. Return row, col, val arrays (also known as COO arrays)

    ii,jj,vv = find(s)

Arithmetic Operations

  1. Add Subtract

    s = s1 + s2
    s = s1 - s2
  2. Scalar Product

    s = 2.0 * s1
    s = s1 / 2.0
  3. Sparse Product

    s = s1 * s2
  4. Transposition

    s = s1'

Sparse Solvers

  1. Solve a linear system (s is a square matrix)

    sol = s\rhs
  2. Solve a least square system (s is a tall matrix)

    sol = s\rhs

The least square solvers are implemented using Eigen sparse linear packages, and the gradients are also implemented. Thus, the following codes will work as expected (the gradients functions will correctly compute the gradients):

ii = [1;2;3;4]
jj = [1;2;3;4]
vv = constant([1.0;1.0;1.0;1.0])
rhs = constant(rand(4))
s = SparseTensor(ii, jj, vv, 4, 4)
sol = s\rhs
run(sess, sol)
run(sess, gradients(sum(sol), rhs))
run(sess, gradients(sum(sol), vv))

Assembling Sparse Matrix

In many applications, we want to accumulate row, col and val to assemble a sparse matrix in iterations. For this purpose, we provide the SparseAssembler utilities.

SparseAssembler(handle::Union{PyObject, <:Integer}, n::Union{PyObject, <:Integer}, tol::Union{PyObject, <:Real}=0.0)

Creates a SparseAssembler for accumulating row, col, val for sparse matrices.

  • handle: an integer handle for creating a sparse matrix. If the handle already exists, SparseAssembler return the existing sparse matrix handle. If you are creating different sparse matrices, the handles should be different.
  • n: Number of rows of the sparse matrix.
  • tol (optional): Tolerance. SparseAssembler will treats any values less than tol as zero.

Example 1

handle = SparseAssembler(100, 5, 1e-8)
op1 = accumulate(handle, 1, [1;2;3], [1.0;2.0;3.0])
op2 = accumulate(handle, 2, [1;2;3], [1.0;2.0;3.0])
J = assemble(5, 5, [op1;op2])

J will be a SparseTensor object.

Example 2

handle = SparseAssembler(0, 5)
op1 = accumulate(handle, 1, [1;2;3], ones(3))
op2 = accumulate(handle, 1, [3], [1.])
op3 = accumulate(handle, 2, [1;3], ones(2))
J = assemble(5, 5, [op1;op2;op3]) # op1, op2, op3 are parallel
Array(run(sess, J))≈[1.0  1.0  2.0  0.0  0.0
                1.0  0.0  1.0  0.0  0.0
                0.0  0.0  0.0  0.0  0.0
                0.0  0.0  0.0  0.0  0.0
                0.0  0.0  0.0  0.0  0.0]
accumulate(handle::PyObject, row::Union{PyObject, <:Integer}, cols::Union{PyObject, Array{<:Integer}}, vals::Union{PyObject, Array{<:Real}})

Accumulates row-th row. It adds the value to the sparse matrix

for k = 1:length(cols)
    A[row, cols[k]] += vals[k]

handle is the handle created by SparseAssembler.

See SparseAssembler for an example.


The function accumulate returns a op::PyObject. Only when op is executed, the nonzero values are populated into the sparse matrix.

assemble(m::Union{PyObject, <:Integer}, n::Union{PyObject, <:Integer}, ops::PyObject)

Assembles the sparse matrix from the ops created by accumulate. ops is either a single output from accumulate, or concated from several ops

op1 = accumulate(handle, 1, [1;2;3], [1.0;2.0;3.0])
op2 = accumulate(handle, 2, [1;2;3], [1.0;2.0;3.0])
op = [op1;op2] # equivalent to `vcat([op1, op2]...)`

m and n are rows and columns of the sparse matrix.

See SparseAssembler for an example.