Example D1: Deflation to Locate Multiple roots
It is sometimes the case that we are trying to solve for the solution of a system of equations that has multiple zeros where one is (or a few are) already known. Case in point in dynamical systems: computing bifurcated branches emanating from a "main" branch of solutions.
We use the method of deflation to discourage the solver from converging to the previously known solution(s) by penalizing the residue by the distance to that solution. Let us suppose that $u_0$ is the previously known solution point and $R(u)$ is the residue function that needs to be solved ($u\in \mathbb{R}^n$, $R: \mathbb{R}\to \mathbb{R}$). We write the deflated residue
\[R_d(u) = \frac{1}{||u-u_0||_2^2} R(u),\]
where $||x||_2$ denotes the 2-norm. The limit $u\to u_0$ may or may not exist for $R_d(u)$ based on the speed at which $R(u)$ goes to zero at $u_0$, and this may affect the behavior of $R_d(u)$. Notwithstanding this, we find using the 2-norm as above to be quite effective for most problems.
\[R_d(u)\]
discourages the solver from approaching the known solution $u_0$, allowing us to converge to another solution if it exists and if our initial guess is close enough.
juliajim implements DEFLATEDRES! which is a deflation wrapper routine that can be used on top of any existing residue function to achieve deflation.
This example demonstrates this through a simple 2-state system of polynomial equations with two known roots:
\[R(u; p) = \begin{bmatrix} u_2\\-p u_1-u_2+\frac{u_1^3}{4} \end{bmatrix},\]
which is solved by $(u_1,u_2)=(-2\sqrt{p},0), (0,0), (2\sqrt{p},0)$.
Load Packages
using NonlinearSolve
using LinearAlgebra
using juliajim.MDOFUTILSSetup Residue function
function fres!(up, r=nothing, j=nothing, jp=nothing)
if !(r === nothing)
r[:] = [up[2];-up[3]*up[1]-up[2]+up[1]^3/4];
end
if !(j === nothing)
j[:, :] = [0 1;-up[3]+3up[1]^2/4 -1];
end
if !(jp === nothing)
jp[:] = [0;-1];
end
endfres! (generic function with 4 methods)Obtain the first solution
par = 3.0;
u0 = [3.0,0];
fun1 = NonlinearFunction((r,u,p)->fres!([u;p], r),
jac=(J,u,p)->fres!([u;p], nothing,J),
paramjac=(Jp,u,p)->fres!([u;p], nothing,nothing,Jp));
prob1 = NonlinearProblem(fun1, u0, par);
sol1 = solve(prob1, show_trace=Val(true));
Algorithm: NewtonRaphson(
descent = NewtonDescent(),
autodiff = AutoForwardDiff(),
vjp_autodiff = AutoFiniteDiff(
fdtype = Val{:forward}(),
fdjtype = Val{:forward}(),
fdhtype = Val{:hcentral}(),
dir = true
),
jvp_autodiff = AutoForwardDiff(),
concrete_jac = Val{false}()
)
---- ------------- -----------
Iter f(u) inf-norm Step 2-norm
---- ------------- -----------
0 2.25000000e+00 0.00000000e+00
1 8.64000000e-01 6.00000000e-01
2 4.41013120e-02 1.28571429e-01
3 1.38792329e-04 7.30382447e-03
4 1.39016088e-09 2.31315914e-05
5 1.77635684e-15 2.31693479e-10
Final 1.77635684e-15
----------------------This returns the solution
sol1retcode: Success
u: 2-element Vector{Float64}:
3.464101615137755
0.0which is exactly $(2\sqrt{p},0)$ as expected.
Deflated Solution
uC = [sol1.u; par];
fun2 = NonlinearFunction((r,u,p)-> DEFLATEDRES!([u;p], uC, fres!; R=r),
jac=(J,u,p)-> DEFLATEDRES!([u;p], uC, fres!; J=J),
paramjac=(Jp,u,p)-> DEFLATEDRES!([u;p], uC, fres!; Jp=Jp));
prob2 = NonlinearProblem(fun2, u0, par);
sol2 = solve(prob2, show_trace=Val(true));
Algorithm: NewtonRaphson(
descent = NewtonDescent(),
autodiff = AutoForwardDiff(),
vjp_autodiff = AutoFiniteDiff(
fdtype = Val{:forward}(),
fdjtype = Val{:forward}(),
fdhtype = Val{:hcentral}(),
dir = true
),
jvp_autodiff = AutoForwardDiff(),
concrete_jac = Val{false}()
)
---- ------------- -----------
Iter f(u) inf-norm Step 2-norm
---- ------------- -----------
0 1.04461524e+01 0.00000000e+00
1 4.73423048e+00 3.78395951e-01
2 1.98329177e+00 5.77129827e-01
3 7.33414117e-01 7.27236432e-01
4 2.22548715e-01 6.97298275e-01
5 4.60853796e-02 4.52588317e-01
6 3.72572734e-03 1.52574887e-01
7 3.12508123e-05 1.46513168e-02
8 2.25490352e-09 1.24985209e-04
9 1.17423574e-17 9.01961399e-09
Final 1.17423574e-17
----------------------This returns the trivial solution $(0,0)$ for the above initial guess:
sol2retcode: Success
u: 2-element Vector{Float64}:
4.6969429454406597e-17
0.0Adding multiple deflators
DEFLATEDRES! allows specifying multiple solutions (in a list) for deflation. Let us specify both the solutions found above.
uCs = [[sol1.u;par], [sol2.u;par]];
fun3 = NonlinearFunction((r,u,p)-> DEFLATEDRES!([u;p], uCs, fres!, R=r),
jac=(J,u,p)-> DEFLATEDRES!([u;p], uCs, fres!, J=J),
paramjac=(Jp,u,p)-> DEFLATEDRES!([u;p], uCs, fres!, J=J));
prob3 = NonlinearProblem(fun3, u0, par);
sol3 = solve(prob3, show_trace=Val(true));
Algorithm: NewtonRaphson(
descent = NewtonDescent(),
autodiff = AutoForwardDiff(),
vjp_autodiff = AutoFiniteDiff(
fdtype = Val{:forward}(),
fdjtype = Val{:forward}(),
fdhtype = Val{:hcentral}(),
dir = true
),
jvp_autodiff = AutoForwardDiff(),
concrete_jac = Val{false}()
)
---- ------------- -----------
Iter f(u) inf-norm Step 2-norm
---- ------------- -----------
0 1.06961524e+01 0.00000000e+00
1 5.09526538e+00 3.95815545e-01
2 2.72599526e+00 7.10468839e-01
3 5.01974668e+00 2.48624726e+00
4 2.41439058e+00 5.78303591e-01
5 9.61732005e-01 1.00218788e+00
6 1.76376748e-01 1.01608731e+00
7 4.70383197e-03 2.67469974e-01
8 3.07342153e-06 7.51630224e-03
9 1.30876791e-12 4.91747026e-06
10 1.85037171e-16 2.09402865e-12
Final 1.85037171e-16
----------------------Now we recover the solution $(-2\sqrt{p},0)$!
sol3retcode: Success
u: 2-element Vector{Float64}:
-3.4641016151377544
0.0Now isn't that cool? We've used exactly the same initial guess and have progressively recovered all the solutions of the system.
What if we kept going?
In the practical setting we are sometimes not destined to know a priori how many solutions to expect. So let us keep going to see what happens.
uCs = [[s.u;par] for s in [sol1,sol2,sol3]];
fun4 = NonlinearFunction((r,u,p)-> DEFLATEDRES!([u;p], uCs, fres!, R=r),
jac=(J,u,p)-> DEFLATEDRES!([u;p], uCs, fres!, J=J),
paramjac=(Jp,u,p)-> DEFLATEDRES!([u;p], uCs, fres!, J=J));
prob4 = NonlinearProblem(fun4, u0, par);
sol4 = solve(prob4, show_trace=Val(true))retcode: Stalled
u: 2-element Vector{Float64}:
-24810.65311212503
-3.81241509304851e12NonlinearSolve will try out its list of methods and then finally get stalled. While this implies the absence of another solution in this case (because we know this), in the general case this merely implies the absence of another solution in the basin of attraction of the provided initial guess. The general problem of numerically eliminating possible solutions is quite complex!
This page was generated using Literate.jl.