Sparsity! spits out a full matrix while the jacobian is sparse

[FuZhiyu]

I met an issue while experimenting with sparsity!. My application commonly involves constructing a sparse matrix (A) based on the input vectors (x) and then updating x by A (updating an HJB equation, for those in econ). A minimal example looks like this:

function testsparse!(out, x)
    A = Tridiagonal(x[2:end], x, x[1:end-1])
    out .= A * x
    return nothing
end
x = [1:100;]
out = similar(x)
sparsity!(testsparse!, out, x)

It can be validated by ForwardDiff that the jacobian is a tridiagonal matrix. However, as reported in issue #7, sparsity returns a matrix with zero entries.

A new bug arises when I compute A * x by hand:

function testsparse2!(out, x)
    A = Tridiagonal(x[2:end], x, x[1:end-1])

    for i in 1:100
        out[i] = 0
        for j in 1:100
            out[i]+= A[i, j] * x[j]
        end
    end
    return nothing
end
x = [1:100;]
out2 = similar(x)

ForwardDiff.jacobian(testsparse2!, out2, x)
sparsity!(testsparse2!, out, x)

ForwardDiff returns the same jacobian as that for testsparse!, while sparsity! returns a 100 x 100 matrix with all (10,000) entries being true.

My julia version is 1.1.1, if it matters. Thanks in advance!

[ChrisRackauckas]

I think this is the same issue and I think it's ranges. @shashi has been off for a bit but when he's back we can see if he can fix it.

[shashi]

Reduced this to:

julia> using LinearAlgebra: MulAddMul, _modify!

julia> function foo(y,x)
           _modify!(MulAddMul(), x[1], y, 1)
       end
foo (generic function with 1 method)

julia> sparsity!(foo, x,y)
Explored path: SparsityDetection.Path(Bool[], 1)
1×1 SparseArrays.SparseMatrixCSC{Bool,Int64} with 0 stored entries

[shashi]

The second case is expected since you're indexing into all elements of the array, not only the non-zero ones.

[jlchan]

Just wanted to check - is this the same issue? I pick up a zero sparsity pattern with a sparse matvec function:

Q = droptol!(sparse(round.(randn(10,10))),1e-10) # random sparse matrix
in = rand(10)
out = similar(in)
function spmatvec(dx,x)
  dx .= Q*x
end
sparsity_pattern = sparsity!(spmatvec,out,in)

The output is:

Explored path: SparsityDetection.Path(Bool[], 1)
10×10 SparseMatrixCSC{Bool,Int64} with 0 stored entries

If I use a dense matrix instead, it throws an error:

ERROR: LoadError: error compiling recurse: unsupported or misplaced expression "replaceglobalref" in function recurse

Reduced this to:

julia> using LinearAlgebra: MulAddMul, _modify!

julia> function foo(y,x)
           _modify!(MulAddMul(), x[1], y, 1)
       end
foo (generic function with 1 method)

julia> sparsity!(foo, x,y)
Explored path: SparsityDetection.Path(Bool[], 1)
1×1 SparseArrays.SparseMatrixCSC{Bool,Int64} with 0 stored entries