Lazy W operator support #443

MSeeker1340 · 2018-07-19T22:46:24Z

Implements the lazy W operators as suggested in #416 (comment).

I decided to write WOperator in OrdinaryDiffEq.jl so that DiffEqOperators is not required for OrdinaryDiffEq. The fact that the interface for AbstractDiffEqLinearOperator in DiffEqBase helps a lot.

(By the way, it occurs to me that a lot of methods in common_defaults.jl for DiffEqOperators should, in fact, be written in DiffEqBase instead. I should do a separate PR for this in the future).

Currently the lazy part of WOperator is already in place and what remains is to modify alg_cache and calc_W! to use it. (Also I need to migrate mass_matrix to *DEFunction in DiffEqBase). The suggestions here #416 (comment) is also worth looking into.

I also intend to add more tests for the methods in derivative_utils.jl to utility_tests.jl, starting with calc_J! and calc_W!.

coveralls · 2018-07-19T23:50:53Z

Coverage remained the same at 92.289% when pulling 284fa70 on MSeeker1340:lazy-W into 6eeabcb on JuliaDiffEq:master.

codecov · 2018-07-20T19:32:52Z

Codecov Report

Merging #443 into master will decrease coverage by 0.12%.
The diff coverage is 94.64%.

@@            Coverage Diff             @@
##           master     #443      +/-   ##
==========================================
- Coverage   95.38%   95.25%   -0.13%     
==========================================
  Files          78       78              
  Lines        9379     9483     +104     
==========================================
+ Hits         8946     9033      +87     
- Misses        433      450      +17

Impacted Files	Coverage Δ
src/perform_step/sdirk_perform_step.jl	`86.17% <ø> (ø)`	⬆️
src/caches/adams_bashforth_moulton_caches.jl	`100% <100%> (ø)`	⬆️
src/caches/sdirk_caches.jl	`100% <100%> (ø)`	⬆️
src/caches/rosenbrock_caches.jl	`99.03% <100%> (+0.08%)`	⬆️
src/caches/kencarp_kvaerno_caches.jl	`100% <100%> (ø)`	⬆️
src/caches/bdf_caches.jl	`100% <100%> (ø)`	⬆️
src/caches/euler_imex_caches.jl	`100% <100%> (ø)`	⬆️
src/derivative_utils.jl	`87.34% <82.69%> (-12.66%)`	⬇️
src/dense/low_order_rk_addsteps.jl	`90.51% <0%> (-2.73%)`	⬇️
... and 4 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 6eeabcb...284fa70. Read the comment docs.

MSeeker1340 · 2018-07-20T19:34:27Z

I will be using ImplicitEuler as my test solver as opposed to KenCarp3 as suggested in #416 (comment) because apparently KenCarp3 doesn' support problems with custom mass matrices at the moment.

(Speaking of which, I should also remind myself to give a similar warning for the ExpRK methods.)

The out-of-place version turns out to be fairly simple: all we need is add an additional clause in calc_W! that checks if a lazy W can be constructed (having a AbstractDiffEqLinearOperator jac_prototype), and construct one every time it is called. No change is required for perform_step! or diffeq_nlsolve!.

I'm still using the mass_matrix field in the problem for the moment, as I don't want to change too many things at once. A simple test script can be found here and the results are:

-------------------------
Linear univariate problem
Finite difference error = 0.000756116458477063
Exact Jacobian error = 0.000756116458477063
-------------------------
Linear bivariate problem
Finite difference error = 0.001981690609776791
Exact Jacobian error = 0.001981690609776791
-------------------------
Linear bivariate problem with mass matrix
Finite difference error = 0.2085322284227815
Exact Jacobian error = 0.2085322284227815

MSeeker1340 · 2018-07-20T19:40:55Z

@ChrisRackauckas I have a question regarding the W_transform parameter in calc_W!.

If I understand it correctly, a transformed W operator is

$$ W = \frac{1}{\gamma}MM - J $$

right, basically dividing the original W by $\gamma$? When is this form of W used?

I can probably add a transform field in WOperator to differentiate between the two.

ChrisRackauckas · 2018-07-20T19:42:16Z

I'm still using the mass_matrix field in the problem for the moment, as I don't want to change too many things at once. A simple test script can be found here and the results are:

This might be nice to keep in OrdinaryDiffEqExtendedTests.jl if not the main repo.

right, basically dividing the original W by $\gamma$? When is this form of W used?

Yes, some methods are written in one way and others in the other way. The transformed version is simply better, but I think Rosenbrock23 hasn't upgraded to it yet.

ChrisRackauckas · 2018-07-20T19:44:58Z

I decided to write WOperator in OrdinaryDiffEq.jl so that DiffEqOperators is not required for OrdinaryDiffEq. The fact that the interface for AbstractDiffEqLinearOperator in DiffEqBase helps a lot.

Note we will want to upgrade StochasticDiffEq.jl to get similar implicit handling soon after. We want to keep this part of them in sync.

ChrisRackauckas · 2018-07-20T21:39:20Z

I guess it doesn't need a cache anyone could just build a cache into the linsolve?

MSeeker1340 · 2018-07-22T21:29:01Z

I took an attempt at the second route in #416 (comment) (i.e. migrate the mass_matrix field from *DEProblem to *DEFunction). While this works, I need to change too many things at once both in DiffEqBase and OrdinaryDiffEq. Considering the amount of work left for v0.7 compatibility, I think this might not be a good idea for the time being.

The compromise I came up with is as follows: we keep mass_matrix in *DEProblem for the time being. An additional constructor for WOperator that doesn't take a mass_matrix is included, which is used by alg_cache to build a partially constructed W (the type of its mass matrix being Any). In calc_W!, the mass matrix of W is reset using the mass matrix in the problem. We don't need to worry about linsolve because the default solver, DEFAULT_LINSOLVE, boils down to calling lu! on W and ldiv! on the factorization result, both of which are already supported by FactorizedDiffEqArrayOperator.

Also, instead of dropping J in the cache and using update_coefficients! on the whole of W, calc_J! is used for the update purpose. This is also meant to introduce minimal change to calc_W!.

Simple test script for In-place ImplicitEuler using lazy W can be found here. Output:

Dense W error = 0.0010693101503117756
Lazy W error = 0.0010693101503117756
Test Summary:          | Pass  Total
In-place ImplicitEuler |    3      3

I think the plan going forward should be something like this:

We keep the current implementation of WOperator and calc_W!, so that minimal change to the codebase is required (thus making v0.7 upgrade easier), while still allowing the use of lazy W operator in linsolve.
We do a few more tests for ImplicitEuler to make sure everything is right, then update the other implicit methods, which can be done in a similar way as ImplicitEuler. The PR can be merged after this is done.
After v0.7 upgrade is done, we can go back and resume the original route, as well as taking a shot at Sparse Jacobians #416 (comment) and deciding on what the jacobian/W interface should be for SplitFunction and DynamicalODEFunction.

MSeeker1340 · 2018-07-22T21:31:41Z

I guess it doesn't need a cache anyone could just build a cache into the linsolve?

I think it's more natural to keep WOperator self-contained, e.g. for krylov linsolve we need the mul! method.

ChrisRackauckas · 2018-07-22T21:37:04Z

I took an attempt at the second route in #416 (comment) (i.e. migrate the mass_matrix field from *DEProblem to *DEFunction). While this works, I need to change too many things at once both in DiffEqBase and OrdinaryDiffEq. Considering the amount of work left for v0.7 compatibility, I think this might not be a good idea for the time being.

The upgrade is pretty much done. Be bold, make the change at least to ODEs. The change then needs to be propagated to SDEs and DDEs though.

MSeeker1340 · 2018-07-24T15:33:48Z

Newest version intended for SciML/DiffEqBase.jl#126.

ChrisRackauckas · 2018-07-24T18:08:39Z

src/derivative_utils.jl

@@ -78,7 +78,7 @@ internal cache (can be specified in the constructor; default to regular `Vector`
 It supports all of `AbstractDiffEqLinearOperator`'s interface.
 """
 mutable struct WOperator{T,
-  MType,
+  MType <: Union{UniformScaling,AbstractMatrix},


diffeqoperators are <: AbstractMatrix?

The previously removed concrete W calculation section in `calc_W!` is needed for problems whose function does not have a jacobian but instead calculates it via forward diff/finite difference using `calc_J!`

Building lazy W conflicts with how a *DEFunction with `invW` is handled in `calc_W!`. Thus the cache building criteria in `alg_cache` is changed accordingly.

ChrisRackauckas · 2018-07-29T12:34:19Z

src/derivative_utils.jl

+  if isa(f, Union{SplitFunction, DynamicalODEFunction})
+    error("WOperator does not support $(typeof(f)) yet")
+  end
+  # Convert mass matrix, if needed


That's fine for now. I am hoping to support time-dependent mass matrices but we can just add an issue that it's a feature to implement in the future.

ChrisRackauckas · 2018-07-29T12:54:23Z

src/derivative_utils.jl

+Base.size(W::WOperator, args...) = size(W.J, args...)
+function Base.getindex(W::WOperator, i::Int)
+  if W.transform
+    W.mass_matrix[i] / W.gamma - W.J[i]


Since this shows up everywhere and it's constant given by the algorithm, it might make sense to make this a type parameter so the code can simply dispatch off of it, or at least not have large performance cuts because of this conditioning in indexing. I'm not sure this indexing is really used all that much though.

I'm not sure this indexing is really used all that much though.

This is true both for factorization (it indexes the concrete matrix instead of W) and mul! (which is lazy).

Sounds good

ChrisRackauckas · 2018-07-29T12:58:09Z

src/derivative_utils.jl

+      else
+        J = ForwardDiff.derivative(uf,uprev)
+      end
+      W = mass_matrix - dtgamma*J


out of place isn't using WOperator yet?

It's used in the branch where DiffEqBase.has_jac(f) == true (https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/pull/443/files/21b1a5b7fae562c8f514cc445002453fa916da6c#diff-35ab888388acf2e7239e17f33e5eefe7R257 and https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/pull/443/files/21b1a5b7fae562c8f514cc445002453fa916da6c#diff-35ab888388acf2e7239e17f33e5eefe7R272).

Why not all of the time?

I guess it's not necessary all of the time, and if we don't allow a linsolve option on it then it's pointless for now.

1. Fix global `fun` being polluted by complex_tests.jl 2. Add Compat to REQUIRE (used by iterator_tests.jl)

MSeeker1340 · 2018-07-29T13:48:49Z

Finally got everything to pass. @ChrisRackauckas If everything so far looks good to you, then we can close this after the CI tests finish. I'll be working on StochasticDiffEq in the meantime.

ChrisRackauckas · 2018-07-29T14:16:06Z

test/differentiation_traits_tests.jl

@@ -26,7 +26,7 @@ using OrdinaryDiffEq, Test
    nothing
  end

-  Lotka_f = ODEFunction(Lotka,jac=Lotka_jac,tgrad=Lotka_tgrad)
+  Lotka_f = ODEFunction(Lotka,jac=Lotka_jac,jac_prototype=zeros(2,2),tgrad=Lotka_tgrad)


but nothing still works?

Unfortunately it doesn't now.

Basically, previously we're assuming that jac_prototype is a dense matrix (allocated in alg_cache). However, now that it can in theory be everything but can't be simply determined by the constructor of ODEFunction, we need to explicitly provide one.

we need to have this work to keep the user interface clean. There should just be one spot where there's a hook that nothing means the same as dense.

ChrisRackauckas · 2018-07-29T14:16:31Z

test/REQUIRE

@@ -4,3 +4,4 @@ ParameterizedFunctions 3.0.0
 StaticArrays
 DiffEqCallbacks
 DiffEqOperators
+Compat


where is this used?

iterator_test.jl

ChrisRackauckas · 2018-07-29T14:16:48Z

src/derivative_utils.jl

@@ -240,6 +240,32 @@ function calc_W!(integrator, cache::OrdinaryDiffEqMutableCache, dtgamma, repeat_
      else
        new_W = false
      end
+    else # concrete W using jacobian from `calc_J!`


This branch is taken when?

It's when f neither has invW nor jac, in which case the jacobian is computed using ForwardDiff/finite difference via calc_J!.

The WOperator implementation doesn't yet include this type of jacobian, which is something we could improve on in the future (might be difficult though because jacobian via differencing uses u and uf, and is different than the (J,u,p,t) interface for diffeq operators.

ChrisRackauckas · 2018-07-29T14:19:02Z

An integration test where an ODE is integrated with the mass matrix being a sparse matrix is required to show that this all captures the new functionality.

MSeeker1340 · 2018-07-29T14:26:37Z

An integration test where an ODE is integrated with the mass matrix being a sparse matrix is required to show that this all captures the new functionality.

You're right. I can modify the temporary testing script I wrote a while ago for this purpose.

ChrisRackauckas · 2018-07-29T14:34:56Z

Fixes SciML/DifferentialEquations.jl#259

ChrisRackauckas · 2018-07-29T16:46:37Z

test/utility_tests.jl

+    fun2 = ODEFunction(_f; mass_matrix=mm, jac=(u,p,t) -> t*A)
+    fun1_ip = ODEFunction(_f_ip; mass_matrix=mm)
+    fun2_ip = ODEFunction(_f_ip; mass_matrix=mm,
+      jac_prototype=DiffEqArrayOperator(similar(A); update_func=(J,u,p,t) -> (J .= t .* A; J)))


ChrisRackauckas · 2018-07-29T16:49:15Z

Looks great. Merge this when the tests pass.

MSeeker1340 · 2018-07-29T16:53:04Z

About the integration test: In the temporary test script, I compared solution using 1) lazy W (by providing a jac function/jac prototype) and 2) ForwardDiff (by providing an ODEFunction without jac) against the analytical solution. I have changed it so that the two solutions are compared against each other, which would be more suitable for checking the lazy W implementation indeed work.

I have also used a time-dependent jacobian so that the update functionality can be checked. As per your suggestion, sparse mass matrix and jacobian are used. One caveat is that I have to manually set the algorithm's linsolve to lu instead of the default lu!, because the latter does not support sparse arrays yet.

ChrisRackauckas · 2018-07-29T16:58:33Z

I have also used a time-dependent jacobian so that the update functionality can be checked. As per your suggestion, sparse mass matrix and jacobian are used. One caveat is that I have to manually set the algorithm's linsolve to lu instead of the default lu!, because the latter does not support sparse arrays yet.

https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/linear_nonlinear.jl#L14-L17

We can make the default be nothing and have that smartly check f and use that to determine lu vs lu!. I think that default can handle most cases well (until we get to DAEs, when we will want to make it smarter and automatically switch to qr, and on BandedMatrices it should default to the banded QR, etc). Separate issue.

WOperator and test

71147ee

MSeeker1340 added 7 commits July 20, 2018 12:49

Merge remote-tracking branch 'origin/master' into lazy-W

42ce8c5

Adjust position for WOperator

a4addeb

Update calc_W! for constant cache

197fc40

Scalar WOperator compatibility

1a1f5ad

Mass matrix for finite differece calc_W!

2f17963

Fix non-scalar adaptive ImplicitEuler

df84f70

Merge remote-tracking branch 'origin/master' into lazy-W

c57604b

MSeeker1340 added 2 commits July 20, 2018 15:57

Support transformed W

9f40d55

Merge remote-tracking branch 'origin/master' into lazy-W

66c62aa

MSeeker1340 added 4 commits July 22, 2018 15:11

Merge remote-tracking branch 'origin/master' into lazy-W

2c5850e

Add constructor with uninitiated mass_matrix

f7e425e

Add WOperator support to in-place calc_W!

4bd3617

Update J and W allocation for ImplicitEulerCache

bb185eb

MSeeker1340 added 2 commits July 24, 2018 09:49

Use mass_matrix from ODEFunction

d175f9b

Merge remote-tracking branch 'origin/master' into lazy-W

7b7a1cd

ChrisRackauckas reviewed Jul 24, 2018

View reviewed changes

ChrisRackauckas approved these changes Jul 24, 2018

View reviewed changes

Merge remote-tracking branch 'origin/master' into lazy-W

08f460d

MSeeker1340 added 4 commits July 28, 2018 21:20

Add back concrete W for functions without jac

b8ae169

The previously removed concrete W calculation section in `calc_W!` is needed for problems whose function does not have a jacobian but instead calculates it via forward diff/finite difference using `calc_J!`

Bugfix

e0eea48

Fix invW

72cde28

Building lazy W conflicts with how a *DEFunction with `invW` is handled in `calc_W!`. Thus the cache building criteria in `alg_cache` is changed accordingly.

Fix differentiation traits tests

80d072b

ChrisRackauckas reviewed Jul 29, 2018

View reviewed changes

Fix test scripts

87e174d

1. Fix global `fun` being polluted by complex_tests.jl 2. Add Compat to REQUIRE (used by iterator_tests.jl)

MSeeker1340 changed the title ~~WIP: Lazy W operator support~~ Lazy W operator support Jul 29, 2018

ChrisRackauckas reviewed Jul 29, 2018

View reviewed changes

MSeeker1340 added 3 commits July 29, 2018 12:24

Address jac_prototype=nothing case

90f377b

Merge remote-tracking branch 'origin/master' into lazy-W

e2beccb

Integration test

b520fe0

ChrisRackauckas reviewed Jul 29, 2018

View reviewed changes

ChrisRackauckas added 2 commits July 29, 2018 09:47

Update REQUIRE

c1e53da

Update utility_tests.jl

284fa70

ChrisRackauckas approved these changes Jul 29, 2018

View reviewed changes

ChrisRackauckas mentioned this pull request Jul 29, 2018

Default linear solver problem-adaptive behavior SciML/DifferentialEquations.jl#323

Closed

ChrisRackauckas merged commit 30d3749 into SciML:master Jul 29, 2018

ChrisRackauckas mentioned this pull request Jul 29, 2018

Jacobian types SciML/DifferentialEquations.jl#220

Closed

MSeeker1340 deleted the lazy-W branch July 30, 2018 16:10

Lazy W operator support #443

Lazy W operator support #443

Conversation

MSeeker1340 commented Jul 19, 2018

coveralls commented Jul 19, 2018 • edited

codecov bot commented Jul 20, 2018 • edited

Codecov Report

MSeeker1340 commented Jul 20, 2018

MSeeker1340 commented Jul 20, 2018

ChrisRackauckas commented Jul 20, 2018 • edited

ChrisRackauckas commented Jul 20, 2018

ChrisRackauckas commented Jul 20, 2018 • edited

MSeeker1340 commented Jul 22, 2018

MSeeker1340 commented Jul 22, 2018

ChrisRackauckas commented Jul 22, 2018

MSeeker1340 commented Jul 24, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

MSeeker1340 commented Jul 29, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ChrisRackauckas commented Jul 29, 2018

MSeeker1340 commented Jul 29, 2018

ChrisRackauckas commented Jul 29, 2018

Choose a reason for hiding this comment

ChrisRackauckas commented Jul 29, 2018

MSeeker1340 commented Jul 29, 2018

ChrisRackauckas commented Jul 29, 2018

coveralls commented Jul 19, 2018 •

edited

codecov bot commented Jul 20, 2018 •

edited

ChrisRackauckas commented Jul 20, 2018 •

edited

ChrisRackauckas commented Jul 20, 2018 •

edited