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

Add Halley's method via descent API #404

Open
wants to merge 8 commits into
base: master
Choose a base branch
from

Conversation

tansongchen
Copy link
Contributor

@tansongchen tansongchen commented Apr 8, 2024

Add Halley's method by computing Hessian-vector-vector product with TaylorDiff.jl in $O(n)$ time. Note that this isn't structured as something like "HessianCache", but instead structured as a descent method, which can be viewed as "adding some correction term to the Newton's method". This is mainly because Hessian, either implicitly or explicitly, isn't likely to be used elsewhere.

Tests haven't been added because TaylorDiff.jl should be added to an extension. Once that is done, it is straightforward to add tests.

Checklist

  • Appropriate tests were added
  • Any code changes were done in a way that does not break public API
  • All documentation related to code changes were updated
  • The new code follows the
    contributor guidelines, in particular the SciML Style Guide and
    COLPRAC.
  • Any new documentation only uses public API

Additional context

Add any other context about the problem here.

end
b = cache.b
# compute the hessian-vector-vector product
hvvp = derivative(x -> cache.f(x, cache.p), u, δu, 2)
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
hvvp = derivative(x -> cache.f(x, cache.p), u, δu, 2)
hvvp = derivative(Base.Fix2(cache.f, cache.p), u, δu, 2)

Copy link
Member

Choose a reason for hiding this comment

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

This will also need a case for inplace cache.f

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Got it. TaylorDiff.jl haven't got a in-place mode yet, but will do this week

src/descent/halley.jl Outdated Show resolved Hide resolved
src/algorithms/halley.jl Outdated Show resolved Hide resolved
Copy link

codecov bot commented Apr 9, 2024

Codecov Report

Attention: Patch coverage is 66.66667% with 17 lines in your changes are missing coverage. Please review.

Project coverage is 77.45%. Comparing base (7deb1ed) to head (92264c0).

Files Patch % Lines
src/descent/halley.jl 64.58% 17 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##           master     #404      +/-   ##
==========================================
- Coverage   86.45%   77.45%   -9.01%     
==========================================
  Files          47       49       +2     
  Lines        2872     2918      +46     
==========================================
- Hits         2483     2260     -223     
- Misses        389      658     +269     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@tansongchen
Copy link
Contributor Author

tansongchen commented Apr 26, 2024

@avik-pal Support for in-place f and tests are added. Let me know if you have any questions!

I think there isn't an eqivalent for Fix2 that fixes the third argument, which is needed by evaluating the derivatives given p. Is it ok to use a closure?

@avik-pal
Copy link
Member

I think there isn't an eqivalent for Fix2 that fixes the third argument, which is needed by evaluating the derivatives given p. Is it ok to use a closure?

Yes mark it as @closure and do it

@avik-pal avik-pal linked an issue Apr 30, 2024 that may be closed by this pull request
src/descent/halley.jl Outdated Show resolved Hide resolved
Comment on lines 1 to 16
@testitem "Halley method" begin
f(u, p) = u .* u .- p
f!(fu, u, p) = fu .= u .* u .- p
u0 = [1.0, 1.0]
p = 2.0

# out-of-place
prob1 = NonlinearProblem(f, u0, p)
sol1 = solve(prob1, Halley())
@test sol1.u ≈ [sqrt(2.0), sqrt(2.0)]

# in-place
prob2 = NonlinearProblem(f!, u0, p)
sol2 = solve(prob2, Halley())
@test sol2.u ≈ [sqrt(2.0), sqrt(2.0)]
end
Copy link
Member

Choose a reason for hiding this comment

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

Let's add a bit more comprehensive tests here similar to SimpleNonlinearSolve

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Sure, will do

@tansongchen
Copy link
Contributor Author

I am trying to add this as a block in 23_test_problems_tests.jl so I can make use of these problems. Do you know how to run a @testitem independently? With the vscode built-in test runner, I am getting the following error:

image

@avik-pal
Copy link
Member

avik-pal commented May 2, 2024

Yeah ReTestItems doesn't yet work with VSCode Test Runner. In a julia REPL activate the test using TestEnv then do

using ReTestItems

ReTestItems.runtests(<path to file or test directory>; name="<name of the testitem>")

@tansongchen
Copy link
Contributor Author

Ok I added it to this file near NewtonRaphson. Some cases haven't work mainly due to limitations of TaylorDiff.jl

@tansongchen
Copy link
Contributor Author

@avik-pal Do you have an out-of-place version of 23 problems? When testing with SimpleHallely, I hit

12:48:34 | START (1/1) test item "SimpleHalley" at test/core/23_test_problems_tests.jl:53
┌ Error: ErrorException("SimpleHalley currently only supports out-of-place nonlinear problems")
└ @ Main.var"##RobustnessTesting#247" ~/Public/SimpleNonlinearSolve.jl/test/core/23_test_problems_tests.jl:28

@avik-pal
Copy link
Member

Don't the tests in simplenonlinearsolve use the oop version?

@tansongchen
Copy link
Contributor Author

So currently, the tests at https://github.com/SciML/SimpleNonlinearSolve.jl/blob/main/test/core/23_test_problems_tests.jl doesn't include Halley's method, and when I tried to add Halley in the same way as Newton, I get the error above.

And the file at https://github.com/SciML/DiffEqProblemLibrary.jl/blob/master/lib/NonlinearProblemLibrary/src/NonlinearProblemLibrary.jl doesn't have out-of-place f's

@tansongchen
Copy link
Contributor Author

Ok so I manually ran [1, 5, 8, 15, 16] against SimpleHalley, and 1, 5, 15, 16 also failed. 8 is due to an error on TaylorDiff.jl and has been fixed. I need to release a new version of TaylorDiff to make the CIs pass though.

@tansongchen
Copy link
Contributor Author

Should be fine now, @avik-pal could you review that again?

@avik-pal
Copy link
Member

Handle the test failures, the Misc ones are real

@tansongchen
Copy link
Contributor Author

Ok, tried to fix the precompilation issue and added compat

@tansongchen
Copy link
Contributor Author

Now all fails are unrelated

@avik-pal
Copy link
Member

Are you sure? There are no failures on master. Some weird interaction of Taylor diff with zygote?

@tansongchen
Copy link
Contributor Author

hmm, that's weird. The Zygote stuff in TaylorDiff doesn't have type piracy so shouldn't have weird interactions

@avik-pal avik-pal mentioned this pull request Jun 1, 2024
@avik-pal
Copy link
Member

avik-pal commented Jun 1, 2024

@tansongchen TaylorDiff is definitely causing problems. See #441

@avik-pal avik-pal closed this Jun 22, 2024
@avik-pal avik-pal reopened this Jun 22, 2024
@@ -28,6 +28,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
TaylorDiff = "b36ab563-344f-407b-a36a-4f200bebf99c"
Copy link
Member

Choose a reason for hiding this comment

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

TaylorDiff brings in Symbolics, which seems like quite a heavy dep. We have quite a few downstream users notably OrdinaryDiffEq.

Won't this affect load times? @ChrisRackauckas @oscardssmith thoughts on this? Should we make it a weak dep

Copy link
Member

Choose a reason for hiding this comment

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

Can we make this be an extension alg in a lib like what we're doing in OrdinaryDiffEq.jl? NonlinearSolve.jl is going to have the same issue of just having too many algorithms at some point, and this seems like something that can be a lib extension

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.

Jacobian and Hessian Free Halley's Method
3 participants