Newton’s Method, Deflation and Multiplicity¶
The instructions in the code cells in this chapter require the following functions to be imported:
from phcpy.solutions import make_solution, strsol2dict
from phcpy.deflation import double_newton_step
from phcpy.deflation import double_double_newton_step
from phcpy.deflation import quad_double_newton_step
from phcpy.deflation import double_deflate
from phcpy.deflation import double_multiplicity
Newton’s method¶
Let us start with a simple system of two polynomials in two unknowns.
pols = ['x^2 + 4*y^2 - 4;', '2*y^2 - x;']
At a regular solution, Newton’s method doubles the accuracy in each step. For the example, we start at a solution where the coordinates are given with only three decimal places.
sols = [make_solution(['x', 'y'], [1.23, -0.786])]
print(sols[0])
Here is then the solution with three decimal places of accuracy:
t : 0.000000000000000E+00 0.000000000000000E+00
m : 1
the solution for t :
x : 1.230000000000000E+00 0.0
y : -7.860000000000000E-01 0.0
== err : 0.000E+00 = rco : 1.000E+00 = res : 0.000E+00 =
The function double_newton_step expects a list of solution strings
as its second argument.
sols = double_newton_step(pols, sols)
print(sols[0])
The outcome of one step with Newton’s method in double precision is
t : 0.00000000000000E+00 0.00000000000000E+00
m : 0
the solution for t :
x : 1.23607623318386E+00 0.00000000000000E+00
y : -7.86154018188250E-01 0.00000000000000E+00
== err : 6.230E-03 = rco : 1.998E-01 = res : 3.706E-05 =
Observe the values for err (forward error),
rco (estimated inverse of the condition number), and
res (the residual or backward error).
The multiplicity field m turned 0 because the default
tolerance was not reached and the solution could not be counted
as a proper solution. Let us reset the multiplicity,
making a new solution with the values from the sols[0].
The string representation of the solution is converted
into a dictionary first:
sold = strsol2dict(sols[0])
sold
which yields
{'t': 0j,
'm': 0,
'err': 0.00623,
'rco': 0.1998,
'res': 3.706e-05,
'x': (1.23607623318386+0j),
'y': (-0.78615401818825+0j)}
and then we make a new solution:
sol = make_solution(['x', 'y'], [sold['x'], sold['y']])
print(sol)
As input for the next Newton step we will use
t : 0.000000000000000E+00 0.000000000000000E+00
m : 1
the solution for t :
x : 1.236076233183860E+00 0.000000000000000E+00
y : -7.861540181882500E-01 0.000000000000000E+00
== err : 0.000E+00 = rco : 1.000E+00 = res : 0.000E+00 =
Calling double_newton_step for the second time
sols = [sol]
sols = double_newton_step(pols, sols)
print(sols[0])
gives
t : 0.00000000000000E+00 0.00000000000000E+00
m : 0
the solution for t :
x : 1.23606797751503E+00 0.00000000000000E+00
y : -7.86151377766704E-01 0.00000000000000E+00
== err : 1.090E-05 = rco : 1.998E-01 = res : 1.100E-10 =
Observe that the value of res dropped from magnitude 1.0e-5
down to 1.0e-10, corresponding to the well conditioning of the root.
However, the multiplicy field is still zero because the estimate
for the forward error is still too high.
The third application of the Newton step
sold = strsol2dict(sols[0])
sol = make_solution(['x', 'y'], [sold['x'], sold['y']])
sols = [sol]
sols = double_newton_step(pols, sols)
print(sols[0])
then yields the solution
t : 0.00000000000000E+00 0.00000000000000E+00
m : 1
the solution for t :
x : 1.23606797749979E+00 0.00000000000000E+00
y : -7.86151377757423E-01 0.00000000000000E+00
== err : 2.452E-11 = rco : 1.998E-01 = res : 4.441E-16 =
The value for the residual res is very close to machine precision
and the solution is considered a proper regular solution.
We can double the precision to double double, and continue with the previous approximate solution:
sols = double_double_newton_step(pols, sols)
print(sols[0])
what is printed is below
t : 0.00000000000000000000000000000000E+00 0.00000000000000000000000000000000E+00
m : 1
the solution for t :
x : 1.23606797749978969640917366873130E+00 0.00000000000000000000000000000000E+00
y : -7.86151377757423286069558585843026E-01 0.00000000000000000000000000000000E+00
== err : 3.036E-16 = rco : 1.998E-01 = res : 3.944E-31 =
and doing this once more
sols = double_double_newton_step(pols, sols)\n",
print(sols[0])"
shows
t : 0.00000000000000000000000000000000E+00 0.00000000000000000000000000000000E+00
m : 1
the solution for t :
x : 1.23606797749978969640917366873128E+00 0.00000000000000000000000000000000E+00
y : -7.86151377757423286069558585842966E-01 0.00000000000000000000000000000000E+00
== err : 6.272E-32 = rco : 1.998E-01 = res : 0.000E+00 =
and then on to quad double precision:
sols = quad_double_newton_step(pols, sols)
print(sols[0])
which then gives the solution:
t : 0.0000000000000000000000000000000000000000000000000000000000000000E+00 0.0000000000000000000000000000000000000000000000000000000000000000E+00
m : 1
the solution for t :
x : 1.2360679774997896964091736687312762354406183596115257242708972454E+00 0.0000000000000000000000000000000000000000000000000000000000000000E+00
y : -7.8615137775742328606955858584295892952312205783772323766490197015E-01 0.0000000000000000000000000000000000000000000000000000000000000000E+00
== err : 7.070E-33 = rco : 1.998E-01 = res : 2.101E-64 =
deflation¶
At an isolated singular solution, deflation is a method to restore the quadratic convergence of Newton’s method.
pols = ['(29/16)*x^3 - 2*x*y;', 'x^2 - y;']
Obviously, (0, 0) is a root of the pols but even if we start
very close to (0, 0), Newton’s method converges extremely slowly.
sols = [make_solution(['x', 'y'],[1.0e-6, 1.0e-6])]
print(sols[0])
The initial solution is then
t : 0.000000000000000E+00 0.000000000000000E+00
m : 1
the solution for t :
x : 1.000000000000000E-06 0.0
y : 1.000000000000000E-06 0.0
== err : 0.000E+00 = rco : 1.000E+00 = res : 0.000E+00 =
Running one step with Newton’s method:
sols = double_newton_step(pols, sols)
print(sols[0])
gives
t : 0.00000000000000E+00 0.00000000000000E+00
m : 0
the solution for t :
x : 9.99999906191101E-07 0.00000000000000E+00
y : 9.99999812409806E-13 0.00000000000000E+00
== err : 1.000E-06 = rco : 5.625E-13 = res : 1.875E-19 =
Observe the tiny value for rco, the estimate of the inverse
of the condition number.
The next step with Newton’s method
sols = double_newton_step(pols, sols)
print(sols[0])
shows
t : 0.00000000000000E+00 0.00000000000000E+00
m : 0
the solution for t :\n",
x : 6.66666604160106E-07 0.00000000000000E+00
y : 3.33333270859482E-13 0.00000000000000E+00
== err : 3.333E-07 = rco : 2.778E-14 = res : 1.111E-13 =
which does not brings us much closer to (0, 0).
Now we apply deflation in double precison:
solsd = double_deflate(pols, sols)
print(solsd[0])
which then yields
t : 0.00000000000000E+00 0.00000000000000E+00
m : 1
the solution for t :
x : 9.46532112069346E-24 4.09228221015004E-24
y : 1.02357542351685E-24 -2.03442589046821E-24
== err : 5.292E-12 = rco : 5.608E-03 = res : 1.885E-15 =
Deflation also works on systems with more equations than unknowns.
pols = ['x^2;', 'x*y;', 'y^2;']
Once again, (0, 0) is the obvious root of pols.
sols = [make_solution(['x', 'y'], [1.0e-6, 1.0e-6])]
print(sols[0]
So, we start at
t : 0.000000000000000E+00 0.000000000000000E+00
m : 1
the solution for t :
x : 1.000000000000000E-06 0.0
y : 1.000000000000000E-06 0.0
== err : 0.000E+00 = rco : 1.000E+00 = res : 0.000E+00 =
and apply deflation
sols = double_deflate(pols, sols, tolrnk=1.0e-8)
print(sols[0])
which shows
t : 0.00000000000000E+00 0.00000000000000E+00
m : 1
the solution for t :
x : 1.25000000000000E-07 0.00000000000000E+00
y : 1.25000000000000E-07 0.00000000000000E+00
== err : 1.250E-07 = rco : 8.165E-01 = res : 1.562E-14 =
This is not an improvement, because the tolerance on the numerical
rank, given by tolrnk=1.0e-8 is too severe.
If we relax this tolerance to 1.0e-4,
sols = double_deflate(pols, sols, tolrnk=1.0e-4)
print(sols[0])
then we obtain
t : 0.00000000000000E+00 0.00000000000000E+00
m : 1
the solution for t :
x : -5.87747175411144E-39 -5.14278778484751E-39
y : 2.93873587705572E-39 -1.83670992315982E-40
== err : 2.757E-23 = rco : 4.082E-01 = res : 1.984E-38 =
and the coordinates of this numerical approximation are well below the double precision.
multiplicity structure¶
The multiplicity can be computed locally starting at the solution. Consider the following example:
pols = [ 'x**2+y-3;', 'x+0.125*y**2-1.5;']
The solution (1, 2) is a multiple root of pols.
Let us prepare the input:
sol = make_solution(['x', 'y'], [1, 2])
print(sol)
so, we work with
t : 0.000000000000000E+00 0.000000000000000E+00
m : 1
the solution for t :
x : 1.000000000000000E+00 0.0
y : 2.000000000000000E+00 0.0
== err : 0.000E+00 = rco : 1.000E+00 = res : 0.000E+00 =
as input to double_multiplicity
multiplicity, hilbert_function = double_multiplicity(pols, sol)
The value of multiplicity equals 3
which is confirmed by the content of hilbert_function:
[1, 1, 1, 0, 0, 0]
Thus, the multiplicity of (1, 2) as a root of pols equals three.