From d31b7224e3ac2671666d53ceccb910b8c2f2d7fc Mon Sep 17 00:00:00 2001 From: Maxime Bouyges Date: Mon, 12 Jan 2026 17:57:07 +0100 Subject: [PATCH 1/3] Add DirectLinearSolver support for length-1 vectors --- examples/4_gallery.jl | 123 ++++++++++++++++++++++++++++++++++++++++++ src/settings.jl | 5 +- test/examples.jl | 4 ++ 3 files changed, 131 insertions(+), 1 deletion(-) create mode 100644 examples/4_gallery.jl diff --git a/examples/4_gallery.jl b/examples/4_gallery.jl new file mode 100644 index 00000000..ee455b14 --- /dev/null +++ b/examples/4_gallery.jl @@ -0,0 +1,123 @@ +# # Gallery + +# ## Distance to the unit circle + +#= +We show how to compute a differentiable distance field to the unit circle. Note that, for this specific case, +- an analytical solution is available ($$d(x) = ||x||-1$$), allowing some tests; +- we could directly write an `Optim` function returning the distance function, and this function is differentiable. +=# +using ImplicitDifferentiation +using LinearAlgebra +using Optim +using StaticArrays +using ForwardDiff + +""" Find `θ0` such that the squared distance between `x` and the point `x_circle(θ0)` on the unit circle is minimal. """ +function forward_unit_circle(x) + function f(θ) + _θ = first(θ) + u = x .- SA[cos(_θ), sin(_θ)] + return sum(abs2, u) + end + θ_init = -π # pretend we don't know the solution + results = optimize(f, [θ_init]) + + y = Optim.minimizer(results) # vector of length 1 + z = nothing + + return y, z +end + +""" Gradient, with respect to `θ0` (represented here by `y`) of the squared distance between `x` and `x_circle(θ0)`. """ +function conditions(x, y, _z) + x1, x2 = x + θ = first(y) + c = cos(θ) + s = sin(θ) + ∇₂f = 2 .* [s * (x1 - c) + (-c) * (x2 - s)] + return ∇₂f +end + +#= +Now, let's define three different functions for our distance field: +- one using directly the forward function : `d_explicit`; +- one using an `ImplicitFunction` with an `IterativeLinearSolver`: `d_impl_iter`; +- one using an `ImplicitFunction` with the `DirectLinearSolver`: `d_impl_direct`. +=# +function d_forward(x) + θ = first(first(forward_unit_circle(x))) + dist = norm(x .- SA[cos(θ), sin(θ)]) + return dist +end + +impl_iter = ImplicitFunction(forward_unit_circle, conditions;) +function d_implicit_iter(x) + θ = first(first(impl_iter(x))) + dist = norm(x .- SA[cos(θ), sin(θ)]) + return dist +end + +impl_direct = ImplicitFunction( + forward_unit_circle, + conditions; + representation=MatrixRepresentation(), + linear_solver=DirectLinearSolver(), +) +function d_impl_direct(x) + θ = first(first(impl_direct(x))) + dist = norm(x .- SA[cos(θ), sin(θ)]) + return dist +end + +#= +Evaluated on a point, let's say, `x = [1.1, 0.0]`, they all return the same result (`d(x) ≈ 0.1`): +=# +p = [1.1, 0.0] +@show d_forward(p) +@show d_implicit_iter(p) +@show d_impl_direct(p) +@test d_forward(p) ≈ d_implicit_iter(p) ≈ d_impl_direct(p) #src + +#= +The most interesting part is the differentiation of two implicit functions. Note that the gradient of the +distance to the unit circle is the following unit vector: +```math +\nabla d(x) = \dfrac{1}{\lVert x \rVert} x. +``` +=# +ν_iter = ForwardDiff.gradient(d_implicit_iter, p) +@show ν_iter, norm(ν_iter) +@test abs(ν_iter[1] - 1.0) < eps(1.0) # src +@test abs(ν_iter[2]) < eps(0.0) # src +ν_direct = ForwardDiff.gradient(d_impl_direct, p) +@show ν_direct, norm(ν_direct) +@test abs(ν_direct[1] - 1.0) < eps(1.0) # src +@test abs(ν_direct[2]) < eps(0.0) # src + +#= +Check for another point, `x = [1.1, 1.1]`: +=# +p = [1.1, 1.1] +ν_iter = ForwardDiff.gradient(d_implicit_iter, p) +@show ν_iter, norm(ν_iter) +@test all(abs.(ν_iter .- √(2) / 2) .< 1e-9) # src +ν_direct = ForwardDiff.gradient(d_impl_direct, p) +@show ν_direct, norm(ν_direct) +@test all(abs.(ν_direct .- √(2) / 2) .< 1e-9) # src + +#= +Using the `DirectLinearSolver`, we can even access the hessian of our distance field. Analytically, +it reads: +```math +\nabla(\nabla d)(x) = \dfrac{1}{\lVert x \rVert^3} \begin{pmatrix} + x_2^ & -x_1 x_2 \\ + -x_1 x_2 & x_2^2 \\ +\end{pmatrix}. +``` +=# +p = [1.1, 1.2] +ℋ_iter = ForwardDiff.hessian(d_implicit_iter, p) +ℋ_ref = [(p[2]^2) (-p[1]*p[2]); (-p[1]*p[2]) (p[1]^2)] ./ norm(p)^3 +@show abs.(ℋ_iter .- ℋ_ref) +@test all(abs.(ℋ_iter .- ℋ_ref) .< 1e-9) # src diff --git a/src/settings.jl b/src/settings.jl index c5a3dcf1..3637e17d 100644 --- a/src/settings.jl +++ b/src/settings.jl @@ -19,7 +19,10 @@ Specify that linear systems `Ax = b` should be solved with a direct method. struct DirectLinearSolver <: AbstractSolver end function (solver::DirectLinearSolver)( - A::Union{AbstractMatrix,Factorization}, _Aᵀ, b::AbstractVector, x0::AbstractVector + A::Union{AbstractMatrix,Factorization,Number}, + _Aᵀ, + b::AbstractVector, + x0::AbstractVector, ) return A \ b end diff --git a/test/examples.jl b/test/examples.jl index d62129ba..ae402fa3 100644 --- a/test/examples.jl +++ b/test/examples.jl @@ -13,3 +13,7 @@ end @testitem "Tricks" begin include(joinpath(dirname(@__DIR__), "examples", "3_tricks.jl")) end + +@testitem "Gallery" begin + include(joinpath(dirname(@__DIR__), "examples", "4_gallery.jl")) +end From 09b1edd86be8896bbcc5b64ab4e194741ea47ec7 Mon Sep 17 00:00:00 2001 From: mbouyges Date: Tue, 13 Jan 2026 15:30:47 +0100 Subject: [PATCH 2/3] Fix test --- examples/4_gallery.jl | 37 ++++++++++++++++--------------------- 1 file changed, 16 insertions(+), 21 deletions(-) diff --git a/examples/4_gallery.jl b/examples/4_gallery.jl index ee455b14..dd0a637b 100644 --- a/examples/4_gallery.jl +++ b/examples/4_gallery.jl @@ -45,18 +45,16 @@ Now, let's define three different functions for our distance field: - one using an `ImplicitFunction` with an `IterativeLinearSolver`: `d_impl_iter`; - one using an `ImplicitFunction` with the `DirectLinearSolver`: `d_impl_direct`. =# -function d_forward(x) - θ = first(first(forward_unit_circle(x))) +function d_generic(x, method) + θ = first(first(method(x))) dist = norm(x .- SA[cos(θ), sin(θ)]) return dist end +d_forward(x) = d_generic(x, forward_unit_circle) + impl_iter = ImplicitFunction(forward_unit_circle, conditions;) -function d_implicit_iter(x) - θ = first(first(impl_iter(x))) - dist = norm(x .- SA[cos(θ), sin(θ)]) - return dist -end +d_implicit_iter(x) = d_generic(x, impl_iter) impl_direct = ImplicitFunction( forward_unit_circle, @@ -64,11 +62,7 @@ impl_direct = ImplicitFunction( representation=MatrixRepresentation(), linear_solver=DirectLinearSolver(), ) -function d_impl_direct(x) - θ = first(first(impl_direct(x))) - dist = norm(x .- SA[cos(θ), sin(θ)]) - return dist -end +d_implicit_direct(x) = d_generic(x, impl_direct) #= Evaluated on a point, let's say, `x = [1.1, 0.0]`, they all return the same result (`d(x) ≈ 0.1`): @@ -76,8 +70,8 @@ Evaluated on a point, let's say, `x = [1.1, 0.0]`, they all return the same resu p = [1.1, 0.0] @show d_forward(p) @show d_implicit_iter(p) -@show d_impl_direct(p) -@test d_forward(p) ≈ d_implicit_iter(p) ≈ d_impl_direct(p) #src +@show d_implicit_direct(p) +@test d_forward(p) ≈ d_implicit_iter(p) ≈ d_implicit_direct(p) #src #= The most interesting part is the differentiation of two implicit functions. Note that the gradient of the @@ -88,12 +82,13 @@ distance to the unit circle is the following unit vector: =# ν_iter = ForwardDiff.gradient(d_implicit_iter, p) @show ν_iter, norm(ν_iter) -@test abs(ν_iter[1] - 1.0) < eps(1.0) # src -@test abs(ν_iter[2]) < eps(0.0) # src -ν_direct = ForwardDiff.gradient(d_impl_direct, p) +@show abs(ν_iter[1] - 1.0) +@test abs(ν_iter[1] - 1.0) < 3e-6 # src +@test abs(ν_iter[2]) < 2e-19 # src +ν_direct = ForwardDiff.gradient(d_implicit_direct, p) @show ν_direct, norm(ν_direct) -@test abs(ν_direct[1] - 1.0) < eps(1.0) # src -@test abs(ν_direct[2]) < eps(0.0) # src +@test abs(ν_direct[1] - 1.0) < 3e-6 # src +@test abs(ν_direct[2]) < 2e-19 # src #= Check for another point, `x = [1.1, 1.1]`: @@ -102,7 +97,7 @@ p = [1.1, 1.1] ν_iter = ForwardDiff.gradient(d_implicit_iter, p) @show ν_iter, norm(ν_iter) @test all(abs.(ν_iter .- √(2) / 2) .< 1e-9) # src -ν_direct = ForwardDiff.gradient(d_impl_direct, p) +ν_direct = ForwardDiff.gradient(d_implicit_direct, p) @show ν_direct, norm(ν_direct) @test all(abs.(ν_direct .- √(2) / 2) .< 1e-9) # src @@ -117,7 +112,7 @@ it reads: ``` =# p = [1.1, 1.2] -ℋ_iter = ForwardDiff.hessian(d_implicit_iter, p) +ℋ_iter = ForwardDiff.hessian(d_implicit_direct, p) ℋ_ref = [(p[2]^2) (-p[1]*p[2]); (-p[1]*p[2]) (p[1]^2)] ./ norm(p)^3 @show abs.(ℋ_iter .- ℋ_ref) @test all(abs.(ℋ_iter .- ℋ_ref) .< 1e-9) # src From f3b81f7e3b73ff019150274f1f25d78b712393b5 Mon Sep 17 00:00:00 2001 From: mbouyges Date: Wed, 14 Jan 2026 09:04:20 +0100 Subject: [PATCH 3/3] remove example from gallery --- examples/4_gallery.jl | 118 ------------------------------------------ test/examples.jl | 4 -- test/systematic.jl | 17 ++++++ 3 files changed, 17 insertions(+), 122 deletions(-) delete mode 100644 examples/4_gallery.jl diff --git a/examples/4_gallery.jl b/examples/4_gallery.jl deleted file mode 100644 index dd0a637b..00000000 --- a/examples/4_gallery.jl +++ /dev/null @@ -1,118 +0,0 @@ -# # Gallery - -# ## Distance to the unit circle - -#= -We show how to compute a differentiable distance field to the unit circle. Note that, for this specific case, -- an analytical solution is available ($$d(x) = ||x||-1$$), allowing some tests; -- we could directly write an `Optim` function returning the distance function, and this function is differentiable. -=# -using ImplicitDifferentiation -using LinearAlgebra -using Optim -using StaticArrays -using ForwardDiff - -""" Find `θ0` such that the squared distance between `x` and the point `x_circle(θ0)` on the unit circle is minimal. """ -function forward_unit_circle(x) - function f(θ) - _θ = first(θ) - u = x .- SA[cos(_θ), sin(_θ)] - return sum(abs2, u) - end - θ_init = -π # pretend we don't know the solution - results = optimize(f, [θ_init]) - - y = Optim.minimizer(results) # vector of length 1 - z = nothing - - return y, z -end - -""" Gradient, with respect to `θ0` (represented here by `y`) of the squared distance between `x` and `x_circle(θ0)`. """ -function conditions(x, y, _z) - x1, x2 = x - θ = first(y) - c = cos(θ) - s = sin(θ) - ∇₂f = 2 .* [s * (x1 - c) + (-c) * (x2 - s)] - return ∇₂f -end - -#= -Now, let's define three different functions for our distance field: -- one using directly the forward function : `d_explicit`; -- one using an `ImplicitFunction` with an `IterativeLinearSolver`: `d_impl_iter`; -- one using an `ImplicitFunction` with the `DirectLinearSolver`: `d_impl_direct`. -=# -function d_generic(x, method) - θ = first(first(method(x))) - dist = norm(x .- SA[cos(θ), sin(θ)]) - return dist -end - -d_forward(x) = d_generic(x, forward_unit_circle) - -impl_iter = ImplicitFunction(forward_unit_circle, conditions;) -d_implicit_iter(x) = d_generic(x, impl_iter) - -impl_direct = ImplicitFunction( - forward_unit_circle, - conditions; - representation=MatrixRepresentation(), - linear_solver=DirectLinearSolver(), -) -d_implicit_direct(x) = d_generic(x, impl_direct) - -#= -Evaluated on a point, let's say, `x = [1.1, 0.0]`, they all return the same result (`d(x) ≈ 0.1`): -=# -p = [1.1, 0.0] -@show d_forward(p) -@show d_implicit_iter(p) -@show d_implicit_direct(p) -@test d_forward(p) ≈ d_implicit_iter(p) ≈ d_implicit_direct(p) #src - -#= -The most interesting part is the differentiation of two implicit functions. Note that the gradient of the -distance to the unit circle is the following unit vector: -```math -\nabla d(x) = \dfrac{1}{\lVert x \rVert} x. -``` -=# -ν_iter = ForwardDiff.gradient(d_implicit_iter, p) -@show ν_iter, norm(ν_iter) -@show abs(ν_iter[1] - 1.0) -@test abs(ν_iter[1] - 1.0) < 3e-6 # src -@test abs(ν_iter[2]) < 2e-19 # src -ν_direct = ForwardDiff.gradient(d_implicit_direct, p) -@show ν_direct, norm(ν_direct) -@test abs(ν_direct[1] - 1.0) < 3e-6 # src -@test abs(ν_direct[2]) < 2e-19 # src - -#= -Check for another point, `x = [1.1, 1.1]`: -=# -p = [1.1, 1.1] -ν_iter = ForwardDiff.gradient(d_implicit_iter, p) -@show ν_iter, norm(ν_iter) -@test all(abs.(ν_iter .- √(2) / 2) .< 1e-9) # src -ν_direct = ForwardDiff.gradient(d_implicit_direct, p) -@show ν_direct, norm(ν_direct) -@test all(abs.(ν_direct .- √(2) / 2) .< 1e-9) # src - -#= -Using the `DirectLinearSolver`, we can even access the hessian of our distance field. Analytically, -it reads: -```math -\nabla(\nabla d)(x) = \dfrac{1}{\lVert x \rVert^3} \begin{pmatrix} - x_2^ & -x_1 x_2 \\ - -x_1 x_2 & x_2^2 \\ -\end{pmatrix}. -``` -=# -p = [1.1, 1.2] -ℋ_iter = ForwardDiff.hessian(d_implicit_direct, p) -ℋ_ref = [(p[2]^2) (-p[1]*p[2]); (-p[1]*p[2]) (p[1]^2)] ./ norm(p)^3 -@show abs.(ℋ_iter .- ℋ_ref) -@test all(abs.(ℋ_iter .- ℋ_ref) .< 1e-9) # src diff --git a/test/examples.jl b/test/examples.jl index ae402fa3..d62129ba 100644 --- a/test/examples.jl +++ b/test/examples.jl @@ -13,7 +13,3 @@ end @testitem "Tricks" begin include(joinpath(dirname(@__DIR__), "examples", "3_tricks.jl")) end - -@testitem "Gallery" begin - include(joinpath(dirname(@__DIR__), "examples", "4_gallery.jl")) -end diff --git a/test/systematic.jl b/test/systematic.jl index cbb70055..cfc8abc9 100644 --- a/test/systematic.jl +++ b/test/systematic.jl @@ -19,6 +19,23 @@ using TestItems test_implicit(scen) test_implicit(scen2) end + + # Test for output vector of length 1 + for (linear_solver, backends) in Iterators.product( + [DirectLinearSolver(), IterativeLinearSolver()], + [nothing, (; x=AutoForwardDiff(), y=AutoZygote())], + ) + yield() + scen = Scenario(; + solver=x -> (sqrt.(x), nothing), + conditions=(x, y, z) -> y .^ 2 .- x, + x=[1.0], + implicit_kwargs=(; representation, linear_solver, backends), + ) + scen2 = add_arg_mult(scen) + test_implicit(scen) + test_implicit(scen2) + end end; @testitem "Operator" setup = [TestUtils] begin