Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

New NLsolve #125

Merged
merged 9 commits into from
Dec 25, 2017
Merged

New NLsolve #125

merged 9 commits into from
Dec 25, 2017

Conversation

pkofod
Copy link
Member

@pkofod pkofod commented Oct 27, 2017

This pr closes #110 (now or never) and at the same time allows for fewer allocations in repeated calls to nlsolve without going all the way. That will come later (soon).

The Jacobian (and "fvec" or F(x)) is now in df, so no n-by-n array is allocated each time. This also means that there is no need for DifferentiableSparseMultivariateFunction and DifferentiableGivenSparseMultivariateFunction (... kind of a mouthful), as df is simply parameterized by the type of the Jacobian. The name for the objective type has changed as well to DifferentiableVector.

The choice of keeping the Jacobian and F(x) vector in df means that it has to be seeded, so it complicates matters a bit. I think it's worth it to avoid the Sparse and GivenSparse types, though.

I plan to major tag and put boundaries on dependent packages, and do a small write up on transitioning from one version to the other.

Up next: choose linear solver, dispatch based choice of method, even fewer allocations, docs, more.

@sglyon
Copy link
Collaborator

sglyon commented Oct 28, 2017

Wow! This is a very large effort.

I don't think I will have the time in the coming days to do a full review by myself.

@KristofferC would you be able to help with the review. Maybe @ChrisRackauckas also?

@ChrisRackauckas
Copy link
Contributor

This sounds perfect. I'll review it soon.

README.md Outdated
@@ -87,19 +87,19 @@ nlsolve(f!, initial_x)
```

Alternatively, you can construct an object of type
`DifferentiableMultivariateFunction` and pass it to `nlsolve`, as in:
`DifferentiableVector` and pass it to `nlsolve`, as in:
Copy link
Contributor

Choose a reason for hiding this comment

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

Why this name? It's supposed to be a function-ish object still?

Copy link
Member Author

Choose a reason for hiding this comment

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

Well, in Optim these are just called OnceDifferentiable for example. I think the Function makes in unnecessarily long to write. Here in NLsolve i figured... I'm not sure what I figure, I just found the "Multivariate" terminologically wrong (it's multivalued if anything). Could just call it Differentiable? I don't know. Happy to bikeshed.

README.md Outdated
nlsolve(df, initial_x)
```

If your function `f!` operates not on `AbstractVector` but on arbitrary
`AbstractArray` you have to specify `initial_x` in the constructor of the
`DifferentiableMultivariateFunction`:
`DifferentiableVector`:
Copy link
Contributor

Choose a reason for hiding this comment

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

this makes it a confusing name. So it can handle non-vector arrays?

Copy link
Member Author

Choose a reason for hiding this comment

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

I get your point, requesting bikeshed backup immediately ;)

README.md Outdated
@@ -116,90 +116,90 @@ nlsolve(f!, initial_x, autodiff = true)
### Jacobian available

If, in addition to `f!(x::AbstractVector, fx::AbstractVector)`, you have a
function `g!(x::AbstractVector, gx::AbstractMatrix)` for computing the Jacobian
function `j!(x::AbstractVector, gx::AbstractMatrix)` for computing the Jacobian
Copy link
Contributor

Choose a reason for hiding this comment

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

+1 I like this semantic change.

Copy link
Member Author

@pkofod pkofod Oct 29, 2017

Choose a reason for hiding this comment

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

I see I missed gx in the positional argument in the README, but yeah, g was always "confusing" to me.

src/NLsolve.jl Outdated
include("autodiff.jl")
include("mcp.jl")
df_path = "differentiable_vectors/"
include(df_path*"differentiable_functions.jl")
Copy link
Contributor

Choose a reason for hiding this comment

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

use joinpath instead

Copy link
Collaborator

@KristofferC KristofferC Oct 29, 2017

Choose a reason for hiding this comment

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

I rather just have inlucde("path/file.jl"). That way static analyzers (like the VSCode plugin) can follow the include. It is not powerful enough to "constant propagate" joinpaths (yet).

end
finite_difference_jacobian!(f, x, fx, jx, autodiff)
end
function j!(jx, x)
Copy link
Contributor

Choose a reason for hiding this comment

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

If they aren't going to be handled now, you should open an issue about handling these allocations better. Is there a reason you don't just pre-allocate and enclose the fx here instead of getting a new one every time??

Copy link
Member Author

Choose a reason for hiding this comment

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

I didn't really look at this code besides renaming stuff an making sure it didn't work. I was waiting for DiffEqDiffTools to clean this up.

function fj!(fx, jx, x)
f!(fx, x)
function f(x)
fx = similar(x)
Copy link
Contributor

Choose a reason for hiding this comment

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

pre-allocate an enclose?

Copy link
Member Author

Choose a reason for hiding this comment

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

see above

@ChrisRackauckas
Copy link
Contributor

Do you have a test for something like a BandedMatrix Jacobian?

@pkofod
Copy link
Member Author

pkofod commented Oct 29, 2017

Do you have a test for something like a BandedMatrix Jacobian?

I should add one.

@KristofferC
Copy link
Collaborator

The reason for DifferentiableGivenSparseMultivariateFunction is so that one can give a known sparsity pattern for a sparse Jacobian. Is this still possible?

@pkofod
Copy link
Member Author

pkofod commented Oct 29, 2017

Yes

Edit: to be a bit more specific: you just pass it using the J keyword

@pkofod
Copy link
Member Author

pkofod commented Oct 29, 2017

Instead of DifferentiableVector would something like NLEquations be better? Or DifferentiableEquations or something like that.

@pkofod
Copy link
Member Author

pkofod commented Oct 29, 2017

A simple (very simply) example showing that BandedMatrices work (after waiting forty five hours to precompile BandedMatrices)

julia> function f!(fstor, x)
       fstor.=(1/3).*x.^3
       end
f! (generic function with 1 method)

julia> function j!(jstor, x)
       for i = 1:3
       jstor[i,i] = x[i]^2
       end
       end
j! (generic function with 1 method)

julia> using NLsolve, BandedMatrices
 
julia> dv = DifferentiableVector(f!, j!, ones(3); J = bones(Float64, 3, 0, 0))
NLsolve.DifferentiableVector{#f!,#j!,NLsolve.#fj!#10{#f!,#j!},Array{Float64,1},BandedMatrices.BandedMatrix{Float64}}(f!, j!, NLsolve.fj!, [6.93768e-310, 6.9377e-310, 6.9377e-310], 3×3 BandedMatrices.BandedMatrix{Float64}:
nl 1.0          
      1.0     
           1.0, s[0], [0])

julia> NLsolve.jacobian(dv)
3×3 BandedMatrices.BandedMatrix{Float64}:
 1.0          
      1.0     
           1.0

julia> nlsolve(dv, ones(3))
Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [1.0, 1.0, 1.0]
 * Zero: [0.00228366, 0.00228366, 0.00228366]
 * Inf-norm of residuals: 0.000000
 * Iterations: 15
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 16
 * Jacobian Calls (df/dx): 16

julia> NLsolve.jacobian(dv)
3×3 BandedMatrices.BandedMatrix{Float64}:
 5.2151e-6                      
            5.2151e-6           
                       5.2151e-6

@ChrisRackauckas
Copy link
Contributor

(after waiting forty five hours to precompile BandedMatrices)

JuliaLinearAlgebra/BandedMatrices.jl#39

As of yesterday it was a bad idea lol

@pkofod
Copy link
Member Author

pkofod commented Oct 29, 2017

So, since DifferentiableMultivariateFunction is sort of terminologically wrong, I wanted to change it. Turns out DifferentiableVector is not too popular, and I get it. So what's the best alternative?

  • DifferentiableEquations (too close to DifferentialEquations?)
  • DifferentiableSystem
  • NLEquations
  • NonlinearEquations
  • ...

Happy to take suggestions.

@KristofferC
Copy link
Collaborator

KristofferC commented Oct 29, 2017

Is there a reason not to use the exact same function types that are used in Optim? Edit: Because they assume to return scalars?

@pkofod
Copy link
Member Author

pkofod commented Oct 30, 2017

Edit: Because they assume to return scalars?

Something like that makes is a bit difficult to pull off, but I guess it is doable if we (as I asked a bit about the other day on gitter) separate the cache variables from these types. Then OnceDifferentiable would play the role of DifferentiableMultivariateFunction. The name is fine with me. Just requires another round of updating NLSolversBase ^^

@pkofod
Copy link
Member Author

pkofod commented Oct 30, 2017

Something like that makes is a bit difficult to pull off, but I guess it is doable if we (as I asked a bit about the other day on gitter) separate the cache variables from these types. Then OnceDifferentiable would play the role of DifferentiableMultivariateFunction. The name is fine with me. Just requires another round of updating NLSolversBase ^^

It would requite some uniform way of specifying the Jacobian/gradient cache type (what's a good keyword name here?)

@pkofod
Copy link
Member Author

pkofod commented Oct 31, 2017

@KristofferC @ChrisRackauckas you both have reviewer's debt here, so I'm pinging you. What's your take on simply calling it OnceDifferentiable like in NLSolversBase, and for now you would simply have to specify a namespace/module if you load both?

using NLsolve
f! ...
j! ...
x0 ...
od = OnceDifferentiable(f!, j!, ...)
nlsolve(od, x0)
using Optim
# queue the warning
f = ...
g! = ...
optod = Optim.OnceDifferentiable(f, g!, ...)
optx0 = ...
optimize(optod, x0, ...)

Then I can work on merging the two approaches at a later point (and I think I know what I want to do)

@ChrisRackauckas
Copy link
Contributor

Why not go straight to this being the common approach for both, i.e. the new OnceDifferentiatible?

@pkofod
Copy link
Member Author

pkofod commented Oct 31, 2017

this being the common approach for both

sorry, but what is "this" here? :)

@pkofod
Copy link
Member Author

pkofod commented Oct 31, 2017

You mean fix NLSolversBase and Optim at the same time?

@ChrisRackauckas
Copy link
Contributor

ChrisRackauckas commented Oct 31, 2017

Would you be able to use this exact same type with NLsolve? What would need to be changed?

Edit: I meant Optim.

@pkofod
Copy link
Member Author

pkofod commented Oct 31, 2017

Would you be able to use this exact same type with NLsolve? What would need to be changed?

Some of the logic in NLSolversBase would be moved to Optim so as to make OnceDifferentiable more generic. Some of the "logic" I've baked into the type currently, would then be handled using a separate cache for Optim and one for NLsolve. This cache could then be made manually, or would be appropriately constructed based on the call being an optimize(df, x0) or nlsolve(df, x0) call.

@ChrisRackauckas
Copy link
Contributor

I think you should go straight to an architecture for the cache that works for both then, knowing that the soon after goal is to have both use the same type.

@pkofod
Copy link
Member Author

pkofod commented Oct 31, 2017

I may have asked this before on gitter, but would

OnceDifferentiable(f, g!; derivative = ...) # Optim
OnceDifferentiable(f!, j!; derivative = ...) # NLsolve

make sense if derivative is a keyword used to provide a gradient or Jacobian cache? It's the best word I can come up with that cover both the idea of a gradient (vector of partial derivatives) and Jacobian (matrix of partial derivatives). Then we could just have a common OnceDifferentiable.

TwiceDifferentiable might still have a specific hessian keyword. I don't think we'll do higher order methods here, as far as I can understand they're typically prohibitively expensive and awkward.

It could also just be

OnceDifferentiable(f, g!; df = ...)
OnceDifferentiable(f!, g!; df= ...)

That's much more concise as well. Neither of them incorporates "partial" ( or ∂ 😮 ) as I think it gets too long, so if the latter is preferred we'll have to settle with d :)

@pkofod
Copy link
Member Author

pkofod commented Dec 11, 2017

Okay, pushed a big update that works with JuliaNLSolvers/NLSolversBase.jl#27 . We no longer have the objective type in this package and that gives us a net loss of LOC currently. All tests pass locally although some tests have been turned off for now. Quite a few changes but I believe that

  • When this is done this will be a massive improvement
  • Everything can be deprecated if that's important (I'm not sure I'll have time to do it, so people will probably just have to accept that this is a massive overhaul that they need to understand the consequences off using NLsolve going forward).
  • This will help bring Optim and NLsolve closer together which was always my goal

@@ -13,7 +13,7 @@ function OnceDifferentiable(f!, F::AbstractArray, x::AbstractArray, autodiff = :
F = similar(x)
fj!(F, J, x)
end
return OnceDifferentiable(f!, j!, fj!, similar(x), similar(x))
return OnceDifferentiable(f!, j!, fj!, x, x)
Copy link
Contributor

Choose a reason for hiding this comment

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

is it fine that this aliases?

Copy link
Member Author

Choose a reason for hiding this comment

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

@pkofod pkofod closed this Dec 18, 2017
@pkofod pkofod reopened this Dec 18, 2017
@pkofod pkofod closed this Dec 19, 2017
@pkofod pkofod reopened this Dec 19, 2017
@pkofod pkofod merged commit 150f861 into master Dec 25, 2017
@pkofod pkofod deleted the pkm/df branch December 25, 2017 21:52
@KristofferC
Copy link
Collaborator

Lovely!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Change f!(x, fvec) to f!(fvec, x)
4 participants