confusing Jacobian options...

[JianghuiDu]

I read the document carefully but still can't understand the relationships among the jacobian options and what the default options are. Is it possible to show a chart of how these options are related?

For example, how does autodiff=true inside the algorithm combine with jac_prototype and linsolve?

1. Does autodiff=true do matrix coloring and jacobian-vector product automatically?
2. Does autodiff=true come with the default linsolve=DefaultLinSolve?
3. In my own problem the jacobian is a banded matrix. But when I specify the jac_prototype using Bandedmatrices and autodiff=true , the result is much slower than doing autodiff=true alone. Does this mean specialized banded matrix factorization was not used?
4. Similarly, combination of autodiff=true with jac_prototype=JacVecOperator() is much slower than autodiff=true alone, why?
5. There is also an autodiff option inside JacVecOperator(). Is it redundant to enable both autodiff?

[ChrisRackauckas]

Does autodiff=true do matrix coloring and jacobian-vector product automatically?

Yes it does matrix coloring, no it doesn't do Jv automatically because that's normally much slower.

Does autodiff=true come with the default linsolve=DefaultLinSolve?

Yes

In my own problem the jacobian is a banded matrix. But when I specify the jac_prototype using Bandedmatrices and autodiff=true , the result is much slower than doing autodiff=true alone. Does this mean specialized banded matrix factorization was not used?

It means the OpenBLAS banded matrix solver is slower than RecursiveFactorization.jl

Similarly, combination of autodiff=true with jac_prototype=JacVecOperator() is much slower than autodiff=true alone, why?

Newton-Krylov methods are slow without really good problem-specific preconditioners.

There is also an autodiff option inside JacVecOperator(). Is it redundant to enable both autodiff?

We need to get rid of JacVecOperator entirely.

[JianghuiDu]

Thanks for the explanation!
So when autodiff=true, the FULL Jacobian is computed using matrix coloring, right?
But is it computed using finite difference or automatic differentiation when specifying jac_prototype at the same time?
Is it possible to compute a structured (Banded, known sparse pattern etc.) Jacobian directly using automatic differentiation?
Or is the FULL Jacobian converted to the structured type?

[ChrisRackauckas]

autodiff=true is orthogonal to the other options. If autodiff = true then autodiff is used. If autodiff is false then autodiff is not used. If jac_prototype is set and autodiff = true, then sparse automatic differentiation with matrix coloring is used. If jac_prototype is set and autodiff = false, then sparse numerical differentiation with matrix coloring is used. It's hard to describe it all because there's a lot of codegen going on to handle the full set of combinations, but the arguments are all orthogonal and what you get are the combinations of the ideas.

Is it possible to compute a structured (Banded, known sparse pattern etc.) Jacobian directly using automatic differentiation?

Yes. And that's what happens by default. jac_prototype is the matrix that is actually used as the Jacobian internally. sparsity_pattern is what's used for the AD coloring scheme. By default, sparsity_pattern = jac_prototype, and so it'll both use the matrix you give as the Jacobian and it'll specialize the differentiation to its form. So for example, on a Banded matrix with size m bands, the analytical solution to the coloring will require m differentiation passes to form the entire banded matrix (instead of n = number of columns, so much less!), and then it'll use the BandedMatrix type you give internally so it'll use the banded matrix LU factorization.

So why is your case slower? It's probably because the BandedMatrix solver is slower than our specialized LU factorization. By default it will be using OpenBLAS. However, the differential equation default linear solver on dense matrices is using a special library called RecursiveFactorization.jl. What's the difference? OpenBLAS:

using BenchmarkTools, LinearAlgebra
A = rand(40,40)
@btime lu!(A) # 124.600 μs (2 allocations: 432 bytes)
@btime lu(A) # 128.400 μs (3 allocations: 13.05 KiB)

vs RecursiveFactorization.jl:

using BenchmarkTools, LinearAlgebra, RecursiveFactorization
A = rand(40,40)
@btime RecursiveFactorization.lu!(A) # 3.112 μs (3 allocations: 464 bytes)
@btime RecursiveFactorization.lu(A) # 3.525 μs (3 allocations: 13.05 KiB)

Yes, we have written a pure Julia-based factorization method which outperforms the classic "BLAS is too fast don't even try" by about 100x. And that's what's being called by default. Now, the BandedMatrix stuff is using the OpenBLAS banded solver, so it's not as slow as the one I show there, but the fact that it's still using OpenBLAS at all will eat away at the perform gain. It should still be asymptotically better (it should be an O(m) algorithm compared to the O(n^3) of a dense LU), but for small enough matrices we have optimized the crap out of the pure Julia linear algebra tools that it probably still comes out on top.

See https://www.youtube.com/watch?v=KQ8nvlURX4M for more details on that. The best thing to do would be to create a pure Julia banded matrix solver, but that will take time.

So what would be the truly optimal thing to do? Set the sparsity pattern or the color vector directly to specialize the automatic differentiation but then use the dense matrix solver until we improve the banded matrix solvers.

[JianghuiDu]

Sounds great to have a native Julia BLAS!
So forwarddiff_color_jacobian from SparseDiffTools is the function to compute structured Jacobian using AD? But why does it takes in both jac_prototype and sparsity_pattern? When would these two be different?
Also, what's the difference between JacVecOperator() from DiffEqOperator and JacVec from SparseDiffTools?

[ChrisRackauckas]

So forwarddiff_color_jacobian from SparseDiffTools is the function to compute structured Jacobian using AD? But why does it takes in both jac_prototype and sparsity_pattern? When would these two be different?

When you want to use color differentiation to speed up the AD part, but want to use a dense matrix for the linear solves. Your case would actually be more optimized with that split.

Also, what's the difference between JacVecOperator() from DiffEqOperator and JacVec from SparseDiffTools?

They are pretty much the same 😓, except JacVecOperator is an AbstractDiffEqOperator. But I want to make that automatic in DiffEq.