Skip to content

Automated sparsity detection #10

@ChrisRackauckas

Description

@ChrisRackauckas

@shashi got the first round going. From talks with @vchuravy and all, it seems like there are a few ways to go about this:

  • Trace with ModelingToolkit. Shown here. This is nice because it gets the analytical solution for the Jacobian as well, but can be expensive to form the computational graphs.
  • Shashi's PR sparsity of the jacobian of a program #9 . This improves upon the ModelingToolkit form because it doesn't form a whole computational graph and instead just stores an interaction vector, making it scale much better memory-wise. However, it runs with a given input X, and so it just takes the branches that input sees. If X is a good representative input, then it'll get the true sparsity pattern.
  • Branch elimination. Essentially it's what Shashi's PR was, except at every branch, you just eliminate the branch and inline both sides of it. This could give an overestimate to the sparsity pattern, but being slightly larger is okay (being slightly smaller isn't). However, because this doesn't have a representative value, it won't work well with while loops with state-dependent flow. I think one way to do this is to just inline the while loop (making it run like a do-while false), but that will always be a problem with this method. However, for a large portion of DiffEq codes, this will work and likely be all that's needed.
  • Concolic fuzzing. This will give the full true estimate, but would require a SAT solver. This is likely what you'd want if you were running it only once, but not likely something you'd default to running on the fly.

So I think the API on sparsity should have a dispatch between the different methods since they all have trade-offs. sparsity!(f!, Y, X, method=TraceGraph(),S=Sparsity(length(Y), length(X))). A last method could be a DefaultDectection() that works like:

  • Checks for any control flow. If none, just do TraceGraph().
  • Otherwise, if there's no state-dependent while loops, do BranchEliminationFlow().
  • Otherwise, do CFuzz()

As for tests, the test equations would use (du,u)->f(du,u,p,t) from definitions of differential equations. It would also be nice if there's just a dispatch on DEProblem that remakes the problem with sparse matrix support, but let's leave that for later. The Lorenz equation is good for unit tests: http://docs.juliadiffeq.org/latest/tutorials/ode_example.html#Example-2:-Solving-Systems-of-Equations-1 . For a bigger example, https://github.com/JuliaDiffEq/DiffEqProblemLibrary.jl has a few. The Bruss problem is particularly interesting

https://github.com/JuliaDiffEq/DiffEqProblemLibrary.jl/blob/master/src/ode/brusselator_prob.jl

since it matches a lot of things we'd typically see in GPU-based PDE code. The different forms of the PDE from http://juliadiffeq.org/DiffEqTutorials.jl/html/introduction/optimizing_diffeq_code.html also are interesting and we should make sure we support the optimized and unoptimized forms well.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions