Skip to content

Commit

Permalink
Merge branch 'master' of https://github.com/briochemc/F1Method.jl
Browse files Browse the repository at this point in the history
* 'master' of https://github.com/briochemc/F1Method.jl:
  Update README.md
  Update README.md
  • Loading branch information
briochemc committed May 6, 2019
2 parents 64c43bd + 33413da commit 7e7a899
Showing 1 changed file with 58 additions and 25 deletions.
83 changes: 58 additions & 25 deletions README.md
Original file line number Diff line number Diff line change
@@ -1,8 +1,7 @@

<img src="https://user-images.githubusercontent.com/4486578/57202054-3d1c4400-6fe4-11e9-97d7-9a1ffbfcb2fc.png" alt="logo" title="F1method" align="right" height="200"/>

F-1 Method
==========
# F-1 Method

<p>
<a href="https://doi.org/10.5281/zenodo.2667835">
Expand Down Expand Up @@ -34,51 +33,85 @@ F-1 Method
</a>
</p>

This package implements the F-1 method descibed in Pasquier et al. (2019).
It allows for efficient quasi-auto-differentiation of an objective function defined implicitly by the solution of a steady-state problem represented by a discretized system of nonlinear partial differential equations that takes the form
This package implements the F-1 method described in Pasquier et al. (in preparation).
It allows for efficient quasi-auto-differentiation of an objective function defined implicitly by the solution of a steady-state problem.

<img src="https://latex.codecogs.com/svg.latex?&space;\boldsymbol{F}(\boldsymbol{x},&space;\boldsymbol{p})&space;=&space;0" title="Eq1"/>
Consider a discretized system of nonlinear partial differential equations that takes the form

where
<a href="https://www.codecogs.com/eqnedit.php?latex=\boldsymbol{x}" target="_blank"><img src="https://latex.codecogs.com/gif.latex?\boldsymbol{x}" title="\boldsymbol{x}" /></a>
is a column vector of the model state variables and
<a href="https://www.codecogs.com/eqnedit.php?latex=\boldsymbol{p}" target="_blank"><img src="https://latex.codecogs.com/gif.latex?\boldsymbol{p}" title="\boldsymbol{p}" /></a>
is a vector of parameters.
Specifically, the F-1 method allows for an efficient computation of both the gradient vector and the Hessian matrix of a generic objective function defined by
<a href="https://www.codecogs.com/eqnedit.php?latex=\hat{f}(\boldsymbol{p})&space;\equiv&space;f(\boldsymbol{s},\boldsymbol{p})" target="_blank"><img src="https://latex.codecogs.com/gif.latex?\hat{f}(\boldsymbol{p})&space;\equiv&space;f(\boldsymbol{s},\boldsymbol{p})" title="\hat{f}(\boldsymbol{p}) \equiv f(\boldsymbol{s},\boldsymbol{p})" /></a>
where
<a href="https://www.codecogs.com/eqnedit.php?latex=\boldsymbol{s}(\boldsymbol{p})" target="_blank"><img src="https://latex.codecogs.com/gif.latex?\boldsymbol{s}(\boldsymbol{p})" title="\boldsymbol{s}(\boldsymbol{p})" /></a>
is the steady-state solution of the system, i.e., such that
<a href="https://www.codecogs.com/eqnedit.php?latex=\boldsymbol{F}(\boldsymbol{s},&space;\boldsymbol{p})&space;=&space;0" target="_blank"><img src="https://latex.codecogs.com/gif.latex?\boldsymbol{F}(\boldsymbol{s},&space;\boldsymbol{p})&space;=&space;0" title="\boldsymbol{F}(\boldsymbol{s}, \boldsymbol{p}) = 0" /></a>.
```
F(x,p) = 0
```

where `x` is a column vector of the model state variables and `p` is a vector of parameters.
The F-1 method then allows for an efficient computation of both the gradient vector and the Hessian matrix of a generic objective function defined by

```
objective(p) = f(s(p),p)
```

where `s(p)` is the steady-state solution of the system, i.e., such that `F(s(p),p) = 0` and where `f(x,p)` is for example a measure of the mismatch between observed state, parameters, and observations.
Optimizing the model is then simply done by minimizing `objective(p)`.
(See Pasquier et al., in prep., for more details.)

## Advantages of the F-1 method

The F-1 method is **easy** to use, gives **accurate** results, and is computationally **fast**:

- **Easy** — The F-1 method basically just needs the user to provide a solver (for finding the steady-state), the mismatch function, `f`, the state function, `F`, and their derivatives, `∇ₓf` and `∇ₓF` w.r.t. the state `x`.
(Note these derivatives can be computed numerically, via the [ForwardDiff](https://github.com/JuliaDiff/ForwardDiff.jl) package for example.)
- **Accurate** — Thanks to dual and hyperdual numbers, the accuracy of the gradient and Hessian, as computed by the F-1 method, are close to machine precision.
(The F-1 method uses the [DualNumbers](https://github.com/JuliaDiff/DualNumbers.jl) and [HyperDualNumbers](https://github.com/JuliaDiff/HyperDualNumbers.jl) packages.)
- **Fast** — The F-1 method is as fast as if you derived and an analytical formulas for every first and second derivatives in *and* used those in the most efficient way.
This is because the bottleneck of such computations is the number of matrix factorizations, and the F-1 method only requires a single one (compared to standard autodifferentiation methods that take the steady-state solver as a black box).

## What's needed?

A requirement of the F-1 method is that the Jacobian matrix <a href="https://www.codecogs.com/eqnedit.php?latex=\mathbf{A}&space;=&space;\nabla_{\boldsymbol{x}}\boldsymbol{F}(\boldsymbol{x},\boldsymbol{p})" target="_blank"><img src="https://latex.codecogs.com/gif.latex?\mathbf{A}&space;=&space;\nabla_{\boldsymbol{x}}\boldsymbol{F}(\boldsymbol{x},\boldsymbol{p})" title="\mathbf{A} = \nabla_{\boldsymbol{x}}\boldsymbol{F}(\boldsymbol{x},\boldsymbol{p})" /></a> can be created, stored, and factorized.
A requirement of the F-1 method is that the Jacobian matrix `A = ∇ₓf` can be created, stored, and factorized.

To use the F-1 method, the user must:

- Make sure that there is a suitable algorithm `alg` to solve the steady-state equation
- overload the `solve` function and the `SteadyStateProblem` constructor from [DiffEqBase](https://github.com/JuliaDiffEq/DiffEqBase.jl). (An example is given in the CI tests — see, e.g., the [`test/simple_setup.jl`](test/simple_setup.jl) file.)
- Provide the derivatives of `f` and `F` with respect to the state, `x`.

### Simple usage
## A concrete example

Make sure you have olverloaded `solve` from DiffEqBase.
Once an initial state, `x₀`, and some parameters, `p₀`, are chosen, simply evaluate the derivatives with
Make sure you have olverloaded `solve` from DiffEqBase
(an example of how to do this is given in the [documentation](https://briochemc.github.io/F1Method.jl/stable/)).
Once initial values for the state, `x₀`, and parameters, `p₀`, are chosen, simply initialize the required memory cache, `mem` via

```julia
# Initialize the cache for storing reusable objects
mem = F1Method.initialize_mem(x₀, p₀)
```

wrap the functions into functions of `p` only via

```julia

# Wrap the objective, gradient, and Hessian functions
objective(p) = F1Method.(f, F, ∇ₓF, mem, p, myAlg(); my_options...)
gradient(p) = F1Method.∇f̂(f, F, ∇ₓf, ∇ₓF, mem, p₀, myAlg(); my_options...)
hessian(p) = F1Method.∇²f̂(f, F, ∇ₓf, ∇ₓF, mem, p₀, myAlg(); my_options...)
gradient(p) = F1Method.∇f̂(f, F, ∇ₓf, ∇ₓF, mem, p, myAlg(); my_options...)
hessian(p) = F1Method.∇²f̂(f, F, ∇ₓf, ∇ₓF, mem, p, myAlg(); my_options...)
```

# Compute the objective function, 𝑓̂(𝒑)
and compute the objective, gradient, or Hessian via either of

```julia
objective(p₀)

# Compute the gradient, ∇𝑓̂(𝒑)
gradient(p₀)

# Compute the Hessian matrix, ∇²𝑓̂(𝒑)
hessian(p₀)
```

That's it.
You were told it was simple, weren't you?
Now you can test how fast and accurate it is!
(Or trust our published work, Pasquier et al., in prep.)

## Citing the software

If you use this package, or implement your own package based on the F-1 method please cite us.
If you use the F-1 method, please cite our Pasquier et al. (in prep.) publication.
If you also use this package directly, please cite us using the [CITATION.bib](./CITATION.bib), which contains a bibtex entry for the software.

3 comments on commit 7e7a899

@briochemc
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register()

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request updated: JuliaRegistries/General/571

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if Julia TagBot is installed, or can be done manually through the github interface, or via:

git tag -a v0.1.0 -m "<description of version>" 7e7a8997f58b141c8469bd4f32728f6cdb68eaca
git push origin v0.1.0

@briochemc
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

TagBot tag

Please sign in to comment.