Solving for the steady-state of a large nonlinear ODE in SciML

Yes, if the solvers do. And they already do: Number is not an AbstractArray.

@rveltz, no need to apologize. The more data we as a collective can give @briochemcm, the better off he will be.

1 Like

Yup, scalar equation solvers would get very upset if fed an array. @rveltz, do you have something in mind?

Things like ArrayFire, ApproxFun…

Thank you! After looking a bit more carefully at the code, I realize that there is no option for an out-of-place Jacobian function, so I need to spend a bit more time with the inner workings of AIBECS to return an in-place Jacobian function before I can try ptcsol (forming the Jacobian out-of-place was never the bottleneck for me so I just skipped trying to implement that until now).

Actually, I realize now that updating a sparse Jacobian inplace is not so straightforward… See, e.g. these many discussions around this topic:

@ctkelley, Are there any plans to allow for out-of-place (sparse) jacobian functions? I ask because I already have available a fast-but-out-of-place jacobian function :man_shrugging:


I think you forgot these alternatives:

  • use of BandedArrays (can really be much faster than sparse and ~ easy to update)
  • use of SparseDiffTools.jl which can be inplace
1 Like

What do you mean? You can’t cook an inplace one from your out-of place one?

First create a sparse matrix with “filled-locations” wherever the Jacobian can be non-zero. Then during your in-place update, fill those locations. If one of them is ==0, that’s ok, just set it to 0 (and make sure that the method you use for setting it, does not delete the “filled-location” spot).

But yes, consider SparseDiffTools.jl

1 Like

Once you’ve allocated the storage, ptcsol defaults to the usual sparse factorization from SparseSuite. I’ll have an example of this in a later chapter of the notebook.

For now I’m using special structrues. There’s an example for PTC and the buckling beam in

where the Jacobian is tridiagonal. There is also an Newton’s method example using BandedMatrices.jl. You need to know the sparsity pattern in advance so the solver does not allocate for the Jacobian. Let me think a bit about how hard it would be for me to use an allocating Jacobian evaluation function.

@rveltz was correct on both of his points. If you Jacobian is banded, using something like BandedMatrices.jl will make your life better and sparse differencing can figure out the sparsity pattern automatically.

I will tag the latest version of the notebook today and will register the package by the end of this week.

1 Like

Well I think I could, it’s just a bit more complicated than I anticipated because my out-of-place Jacobian function drops zeros (the non-structural entries that are sometimes zero). This is both good (because there are quite a few) and bad (because some branching inside the nonlinear function means some non-structural zeros appear and disappear and make it hard to use for an inplace version).

For the record and because it helps me think about, I’ll desribe my Jacobian function here and I might eventually copy this into the package docs anyway. Also maybe someone has good ideas to improve it and make it inplace, too, so here goes. This goes in the details so it is a little wordy — so apologies in advance.

TLDR: I have a function F that looks like F(x) = G(x) - T x where G is a perfect candidate for ForwardDiff and SparseDiffTools, while T is already available as the Jacobian of x \mapsto T x. How can I make an inplace Jacobian for J(x) = \nabla G(x) - T?

Long version:

When trying to solve for F(x) = 0, the variable x is actually a concatenation of smaller (but still large) vectors:

x = \begin{bmatrix} x_1\\ x_2\\ \vdots \\ x_N \\ \end{bmatrix}

where each x_i can be a ~200,000 sized column vector.

All the x_i have the same length and share the same indices, in the sense that these indices correspond to some linear indices of a 3D grid. Each different x_i thus actually represents a different 3D field that has been “vectorized”. The function F can be similarly split

F(x) = \begin{bmatrix} F_1(x_1, \ldots , x_N) \\ F_2(x_1, \ldots , x_N) \\ \vdots \\ F_N(x_1, \ldots , x_N) \\ \end{bmatrix}.

and each F_i is further split into a “nonloncal” linear part and a “local” nonlinear part, like so

F_i(x) = G_i(x_1, \ldots , x_N) - T_i x_i

where the T_i are sparse matrices but not diagonal, and their sparsity structure does not correspond to any of the sparse array packages I have seen (they have 1 to ~20 off-diagonals that can be quite far off the main diagonal and whose distance to the diagonal can vary with each row/column). The G_i, on the other hand, are nonlinear, but local, in that the j-th entry [G_i(x_1, \ldots , x_N)]_j only depends on the j-th entries of its x_i arguments. The AIBECS.jl user supplies the T_i and the G_i functions, which I then use to build F and its Jacobian as a function of a single larger vector x.

The jacobian of F can be built efficiently using autodifferentiation of the G_i only (ForwardDiff.jl in my case). If I denote by \partial_j G_i the “derivative” of G_i with respect to x_j to be the (square) jacobian of the x_j \mapsto G_i(x_1,\ldots,x_N) function, which is diagonal because the G_i are “local” as explained just above, then jacobian of F can be written as

J(x) = \begin{bmatrix} (\partial_1 G_1)(x)-T_1 & (\partial_2 G_1)(x) & \cdots & (\partial_N G_1)(x)\\ (\partial_1 G_2)(x) & (\partial_2 G_2)(x) -T_2 & \cdots & (\partial_N G_2)(x)\\ \vdots & & \ddots & \vdots \\ (\partial_1 G_N)(x) & (\partial_2 G_N)(x) & \cdots & (\partial_N G_N)(x) - T_N\\ \end{bmatrix}

In practice, I actually build it with

J(x) = \begin{bmatrix} (\partial_1 G_1)(x) & (\partial_2 G_1)(x) & \cdots & (\partial_N G_1)(x)\\ (\partial_1 G_2)(x) & (\partial_2 G_2)(x) & \cdots & (\partial_N G_2)(x)\\ \vdots & & \ddots & \vdots \\ (\partial_1 G_N)(x) & (\partial_2 G_N)(x) & \cdots & (\partial_N G_N)(x)\\ \end{bmatrix} - \begin{bmatrix} T_1 & & & \\ & T_2 & & \\ & & \ddots & \\ & & & T_N\\ \end{bmatrix}

where N ForwardDiff’d calls to the G_i's are sufficient to build the first term (which is essentially an array of diagonals) with some combination of sparse, Diagonal, hcat, and vcat operations, and where the second term uses the blockdiag function on the T_i. I believe the - operation drops all the non-structural zeros (but I need to double check), which is good and bad, as mentioned in the intro to this post.

I could use SparseDiffTools for the first term by supplying the sparsity pattern, but it is not efficient to use it on the whole F because the T_i fill in a lot of the sparsity pattern and I actually don’t need to use ForwardDiff on the T_i x_i terms since these are linear and the T_i are already provided.

So how should I go about making an inplace version of this Jacobian?

Using FD, you can update the non diagonal blocks inplace. So I would focus on updating the block (1,1) inplace. You save the indices of T_1 which are changed by \partial_1G_1 and write them inplace. This is what I do here

1 Like

@rveltz I just started thinking about all this again and noticed your link does not work anymore. I’m assuming the code was just moved somewhere else — could you update said link? :pray:

id say BifurcationKit.jl/Utils.jl at master · rveltz/BifurcationKit.jl · GitHub

1 Like

Then it is used here BifurcationKit.jl/PeriodicOrbitTrapeze.jl at master · rveltz/BifurcationKit.jl · GitHub

1 Like

I’ll have a matrix-free version of pseudo-transient continuation ready soon (mid-March at the latest). You’ll need a good preconditioner for it to work well.

Support for sparse matrix factorization PTC will be in the next release (withing days) and is already in the master branch. If you can store a sparse Jacobian, my stuff should do the job for you.

The current release should support general sparse Jacobians but the upcoming release has a working example and does the job for certain.

1 Like

Thinking about this again, is the jacobian a Fredholm operator? I am thinking that T is a differential operator and G is a compact operator. In this case, you can use a GMRES solver with a Matrix-free jacobian using T as a preconditioner. You can save the (incomplete) LU decomposition of T before hand.

I use the same strategy here and here (scroll to the end) or here!

1 Like

The Jacobian (Frechet derivitive) of a nonlinear differential operator, for example, is also a differential operator, which is not Fredholm. So if

F(u) = u'' + f(u) + boundary conditions, then the Frechet derivative acts on a function via
F'(u)v = v'' + f'(u) v. This is an unbounded operator and is not Fredholm

For integral opreators like F(u)(x) = u(x) = \int k(x,y,u) u(y) dy, you get
F'(u)v(x) = v(x) + \int \frac{\partial k}{\partial u} (x,y,u) v(y) dy, which is a compact perturbation of the identity and Fredholm.

1 Like

I haven’t followed the discussion but I think Romain means that (in this case) you precondition with the Laplacian, so (1-Delta)^-1 (u’’ + f’(u)) is Fredholm.

1 Like

Yes I was a bit quick, I mean you can precondition with the differential operator whose resolvent is often compact.

Thanks @antoine-levitt, we were typing at the same time it seems