Start the AD and mutation style extensions. Also document least squar…

Project.toml

@@ -10,13 +10,32 @@ PositiveFactorizations = "85a6dd25-e78a-55b7-8502-1745935b8125"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
+[weakdeps]
+AbstractDifferentiation = "c29ec348-61ec-40c8-8164-b8c60e9d9f3d"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
+SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
+
+[extensions]
+AbstractDifferentiationExt = "AbstractDifferentiation"
+ForwardDiffExt = "ForwardDiff"
+SparseForwardDiffExt = ["ForwardDiff", "SparseDiffTools", "Symbolics"]
+StaticArraysExt = "StaticArrays"
+
 [compat]
+AbstractDifferentiation = "0.6.0"
+ForwardDiff = "0.10.36"
 IterativeSolvers = "0.9"
 PositiveFactorizations = "0.2"
+ReverseDiff = "1.15.1"
 SparseDiffTools = "1.13.2"
+StaticArrays = "1.6.5"
 julia = "1"
 
 [extras]
+AbstractDifferentiation = "c29ec348-61ec-40c8-8164-b8c60e9d9f3d"
 ArbNumerics = "7e558dbc-694d-5a72-987c-6f4ebed21442"
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -26,6 +45,7 @@ IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

docs/mkdocs.yml

@@ -27,6 +27,9 @@ nav:
     - Tutorials:
         - Minimizing a function: 'optimization.md'
         - Solving non-linear equations: 'nonlineareq.md'
+        - Solving non-linear least squares problems: 'leastsquares.md'
+    - Extensions:
+        - Using extensions: 'extensions.md'
 #        - Gradients and Hessians: 'user/gradientsandhessians.md'
 #        - Configurable Options: 'user/config.md'
 #        - Linesearch: 'algo/linesearch.md'

docs/src/extensions.md

@@ -0,0 +1,275 @@
+# Extensions to NLSolvers
+Extensions are feature in Julia that allows package owners to provide functionality (through added methods to library functions) based on other packages. These packages can be useful for some users but cause increased load times for all users. Using extensions we've provided a few convenience functions that do not have to be loaded by all users, but only those who need it.
+
+## ForwardDiff.jl
+ForwardDiff.jl derivatives are available to users by using the `ScalarObjective(ForwardDiffAutoDiff(), x_seed; f=my_f)` constructor. Below is an example that shows that the hand-written and ForwardDiff derivatives agree.
+```
+using NLSolvers, ForwardDiff, Test
+function himmelblau!(x)
+    fx = (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
+    return fx
+end
+function himmelblau_g!(∇f, x)
+    ∇f[1] =
+        4.0 * x[1]^3 + 4.0 * x[1] * x[2] - 44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 - 14.0
+    ∇f[2] =
+        2.0 * x[1]^2 + 2.0 * x[2] - 22.0 + 4.0 * x[1] * x[2] + 4.0 * x[2]^3 - 28.0 * x[2]
+    ∇f
+end
+function himmelblau_fg!(∇f, x)
+    himmelblau!(x), himmelblau_g!(∇f, x)
+end
+function himmelblau_h!(∇²f, x)
+    ∇²f[1, 1] = 12.0 * x[1]^2 + 4.0 * x[2] - 44.0 + 2.0
+    ∇²f[1, 2] = 4.0 * x[1] + 4.0 * x[2]
+    ∇²f[2, 1] = ∇²f[1, 2]
+    ∇²f[2, 2] = 2.0 + 4.0 * x[1] + 12.0 * x[2]^2 - 28.0
+    return ∇²f
+end
+function himmelblau_fgh!(∇f, H, x)
+    himmelblau!(x), himmelblau_g!(∇f, x), himmelblau_h!(H, x)
+end
+
+objective = ScalarObjective(
+    f=himmelblau!,
+    g=himmelblau_g!,
+    fg=himmelblau_fg!,
+	h=himmelblau_h!,
+    fgh=himmelblau_fgh!,
+    )
+
+forward_objective = ScalarObjective(Experimental.ForwardDiffAutoDiff(), [3.0,3.0]; f=himmelblau!)
+
+# Test gradient
+G_forward = zeros(2)
+forward_objective.g(G_forward, [0.1,0.2])
+G = zeros(2)
+objective.g(G, [0.1,0.2])
+using Test
+@test G  ≈ G_forward
+
+# Test joint f and g evaluation
+G_forward = zeros(2)
+fG_forward, _ = forward_objective.fg(G_forward, [0.1,0.2])
+
+G = zeros(2)
+f, _ = objective.fg(G, [0.1,0.2])
+
+@test G  ≈ G_forward
+@test f  ≈ fG_forward
+
+# Test Hessian evaluation
+H_forward = zeros(2,2)
+forward_objective.h(H_forward, [0.1,0.2])
+
+H = zeros(2,2)
+objective.h(H, [0.1,0.2])
+
+@test H  ≈ H_forward
+
+# Test joint f, G, H evaluation
+H_forward = zeros(2,2)
+G_forward = zeros(2)
+f_forward, _, _ = forward_objective.fgh(G_forward, H_forward, [0.1,0.2])
+
+G = zeros(2)
+H = zeros(2,2)
+f, _, _= objective.fgh(G, H, [0.1,0.2])
+
+@test f  ≈ f_forward
+@test G  ≈ G_forward
+@test H  ≈ H_forward
+
+# Test hessian-vector calculations
+# test hv
+x0 = [0.1,0.2]
+hv = x0*0
+forward_objective.hv(hv, x0, 2*ones(2))
+
+objective.h(H, x0)
+@test H*(2.0*ones(2)) ≈ hv
+```
+
+## AbstractDifferentiation
+```
+using NLSolvers, ForwardDiff, Test
+import AbstractDifferentiation as AD
+function himmelblau!(x)
+    fx = (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
+    return fx
+end
+function himmelblau_g!(∇f, x)
+    ∇f[1] =
+        4.0 * x[1]^3 + 4.0 * x[1] * x[2] - 44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 - 14.0
+    ∇f[2] =
+        2.0 * x[1]^2 + 2.0 * x[2] - 22.0 + 4.0 * x[1] * x[2] + 4.0 * x[2]^3 - 28.0 * x[2]
+    ∇f
+end
+function himmelblau_fg!(∇f, x)
+    himmelblau!(x), himmelblau_g!(∇f, x)
+end
+function himmelblau_h!(∇²f, x)
+    ∇²f[1, 1] = 12.0 * x[1]^2 + 4.0 * x[2] - 44.0 + 2.0
+    ∇²f[1, 2] = 4.0 * x[1] + 4.0 * x[2]
+    ∇²f[2, 1] = ∇²f[1, 2]
+    ∇²f[2, 2] = 2.0 + 4.0 * x[1] + 12.0 * x[2]^2 - 28.0
+    return ∇²f
+end
+function himmelblau_fgh!(∇f, H, x)
+    himmelblau!(x), himmelblau_g!(∇f, x), himmelblau_h!(H, x)
+end
+
+objective = ScalarObjective(
+    f=himmelblau!,
+    g=himmelblau_g!,
+    fg=himmelblau_fg!,
+	h=himmelblau_h!,
+    fgh=himmelblau_fgh!,
+    )
+
+forward_objective = ScalarObjective(AD.ForwardDiffBackend(), [3.0,3.0]; f=himmelblau!)
+
+# Test gradient
+G_forward = zeros(2)
+forward_objective.g(G_forward, [0.1,0.2])
+G = zeros(2)
+objective.g(G, [0.1,0.2])
+using Test
+@test G  ≈ G_forward
+
+# Test joint f and g evaluation
+G_forward = zeros(2)
+fG_forward, _ = forward_objective.fg(G_forward, [0.1,0.2])
+
+G = zeros(2)
+f, _ = objective.fg(G, [0.1,0.2])
+
+@test G  ≈ G_forward
+@test f  ≈ fG_forward
+
+# Test Hessian evaluation
+H_forward = zeros(2,2)
+forward_objective.h(H_forward, [0.1,0.2])
+
+H = zeros(2,2)
+objective.h(H, [0.1,0.2])
+
+@test H  ≈ H_forward
+
+# Test joint f, G, H evaluation
+H_forward = zeros(2,2)
+G_forward = zeros(2)
+f_forward, _, _ = forward_objective.fgh(G_forward, H_forward, [0.1,0.2])
+
+G = zeros(2)
+H = zeros(2,2)
+f, _, _= objective.fgh(G, H, [0.1,0.2])
+
+@test f  ≈ f_forward
+@test G  ≈ G_forward
+@test H  ≈ H_forward
+
+# Test hessian-vector calculations
+# test hv
+x0 = [0.1,0.2]
+hv = x0*0
+forward_objective.hv(hv, x0, 2*ones(2))
+
+objective.h(H, x0)
+@test H*(2.0*ones(2)) ≈ hv
+
+
+
+using NLSolvers, Test
+import AbstractDifferentiation as AD
+
+using ReverseDiff
+
+function himmelblau!(x)
+    fx = (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
+    return fx
+end
+function himmelblau_g!(∇f, x)
+    ∇f[1] =
+        4.0 * x[1]^3 + 4.0 * x[1] * x[2] - 44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 - 14.0
+    ∇f[2] =
+        2.0 * x[1]^2 + 2.0 * x[2] - 22.0 + 4.0 * x[1] * x[2] + 4.0 * x[2]^3 - 28.0 * x[2]
+    ∇f
+end
+function himmelblau_fg!(∇f, x)
+    himmelblau!(x), himmelblau_g!(∇f, x)
+end
+function himmelblau_h!(∇²f, x)
+    ∇²f[1, 1] = 12.0 * x[1]^2 + 4.0 * x[2] - 44.0 + 2.0
+    ∇²f[1, 2] = 4.0 * x[1] + 4.0 * x[2]
+    ∇²f[2, 1] = ∇²f[1, 2]
+    ∇²f[2, 2] = 2.0 + 4.0 * x[1] + 12.0 * x[2]^2 - 28.0
+    return ∇²f
+end
+function himmelblau_fgh!(∇f, H, x)
+    himmelblau!(x), himmelblau_g!(∇f, x), himmelblau_h!(H, x)
+end
+
+objective = ScalarObjective(
+    f=himmelblau!,
+    g=himmelblau_g!,
+    fg=himmelblau_fg!,
+	h=himmelblau_h!,
+    fgh=himmelblau_fgh!,
+    )
+
+reverse_objective = ScalarObjective(AD.ReverseDiffBackend(), [3.0,3.0]; f=himmelblau!)
+
+# Test gradient
+G_reverse = zeros(2)
+reverse_objective.g(G_reverse, [0.1,0.2])
+G = zeros(2)
+objective.g(G, [0.1,0.2])
+using Test
+@test G  ≈ G_reverse
+
+# Test joint f and g evaluation
+G_reverse = zeros(2)
+fG_reverse, _ = reverse_objective.fg(G_reverse, [0.1,0.2])
+
+G = zeros(2)
+f, _ = objective.fg(G, [0.1,0.2])
+
+@test G  ≈ G_reverse
+@test f  ≈ fG_reverse
+
+# Test Hessian evaluation
+H_forward = zeros(2,2)
+reverse_objective.h(H_forward, [0.1,0.2])
+
+H = zeros(2,2)
+objective.h(H, [0.1,0.2])
+
+@test H  ≈ H_forward
+
+# Test joint f, G, H evaluation
+H_forward = zeros(2,2)
+G_reverse = zeros(2)
+f_forward, _, _ = reverse_objective.fgh(G_reverse, H_forward, [0.1,0.2])
+
+G = zeros(2)
+H = zeros(2,2)
+f, _, _= objective.fgh(G, H, [0.1,0.2])
+
+@test f  ≈ f_forward
+@test G  ≈ G_reverse
+@test H  ≈ H_forward
+
+# Test hessian-vector calculations
+# test hv
+x0 = [0.1,0.2]
+hv = x0*0
+reverse_objective.hv(hv, x0, 2*ones(2))
+
+objective.h(H, x0)
+@test H*(2.0*ones(2)) ≈ hv
+
+```
+
+## Future
+I plan to support SparseDiffTools.jl and more.