Equation and Variable Scaling
=============================

A polynomial system may have coefficients which vary widely in magnitude.
This may lead to inaccurate solutions.   Equation scaling multiplies 
all coefficients in the same equation by the same constant.  
Variable scaling replaces the original variables by new variables times
some constant.  Both equation and variable scaling can reduce 
the variability of the magnitudes of the coefficients of the system.

In this section, the following functions are applied:

::

    from phcpy.solutions import verify
    from phcpy.solver import solve
    from phcpy.scaling import double_scale_system
    from phcpy.scaling import double_scale_solutions

solving without scaling
-----------------------

Consider the polynomials below.

::

    p = ['0.000001*x^2 + 0.000004*y^2 - 4;', '0.000002*y^2 - 0.001*x;']

Observe the variation in magnitude between the coefficients
of the polynomials in ``p``.  First, the system is solved,
without scaling the coefficients.

::

    psols = solve(p)
    for (idx, sol) in enumerate(psols):
        print('Solution', idx+1, ':')
        print(sol)

Then the solutions are

::

    Solution 1 :
    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x :  1.23606797749980E+03  -4.34288187660045E-10
     y : -7.86151377757428E+02   2.63023405572965E-10
    == err :  2.598E-03 = rco :  4.053E-04 = res :  8.002E-10 =
    Solution 2 :
    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x :  1.23606797749981E+03   1.58919085844426E-11
     y :  7.86151377757436E+02   9.62482268770893E-12
    == err :  5.733E-04 = rco :  4.053E-04 = res :  6.661E-11 =
    Solution 3 :
    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x : -3.23606797750004E+03  -1.73527763271358E-11
     y :  4.58877485485452E-12  -1.27201964951414E+03
    == err :  1.310E-03 = rco :  2.761E-04 = res :  1.855E-10 =
    Solution 4 :
    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x : -3.23606797749979E+03   5.48229606785362E-14
     y :  1.44974035789780E-14   1.27201964951407E+03
    == err :  6.525E-05 = rco :  2.761E-04 = res :  4.063E-14 =

Observe the magnitude of the values for the coordinates of the solutions
and the estimates for the inverse condition numbers in ``rco``.
The forward errors ``err`` are rather large and the residual ``res``
not so close to the machine precision.

The function ``verify`` evaluates the polynomials in the system 
at the solutions and returns the sum of the absolute values.

::

    verify(p, psols)

yields the number

::

    4.3339840153217635e-12

With scaling this error can be reduced.

solving after scaling
---------------------

Equation and variable scaling is applied in double precision.

::

    (q, c) = double_scale_system(p)

Consider now the coefficients of the scaled system:

::

    for pol in q:
        print(pol)

and here are the scaled polynomials:

::

    9.99999999999994E-01*x^2 + 9.99999999999996E-01*y^2 - 1.00000000000000E+00;
      " + 9.99999999999996E-01*y^2 - 9.99999999999997E-01*x;


and scaling coefficients in ``c`` are

::

    [3.30102999566398,
     0.0,
     2.999999999999999,
     0.0,
     -0.6020599913279623,
     0.0,
     -0.30102999566398136,
     0.0,
     0.04028876017114153,
     0.0]

We solve the scaled system:

::

    qsols = solve(q)
    for (idx, sol) in enumerate(qsols):
        print('Solution', idx+1, ':')
        print(sol)

which gives the solutions

::

    Solution 1 :
    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x : -1.61803398874990E+00   2.94545607917864E-90
     y :  1.47272803958932E-90   1.27201964951407E+00
    == err :  1.475E-16 = rco :  2.268E-01 = res :  6.661E-16 =
    Solution 2 :
    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x : -1.61803398874990E+00   2.94545607917864E-90
     y : -1.47272803958932E-90  -1.27201964951407E+00
    == err :  1.475E-16 = rco :  2.268E-01 = res :  6.661E-16 =
    Solution 3 :
    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x :  6.18033988749897E-01  -2.92604772168262E-98
     y :  7.86151377757425E-01   0.00000000000000E+00
    == err :  8.868E-17 = rco :  4.601E-01 = res :  1.110E-16 =
    Solution 4 :
    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x :  6.18033988749897E-01  -2.92604772168262E-98
     y : -7.86151377757425E-01   0.00000000000000E+00
    == err :  8.868E-17 = rco :  4.601E-01 = res :  1.110E-16 =

All solutions of the scaled system are well conditioned 
with small forward and backward errors.

The scaling coefficients in ``c`` are used to bring the solutions 
of the scaled problem to the original coordinates.

::

    ssols = double_scale_solutions(len(q), qsols, c)

The output of

::

    for (idx, sol) in enumerate(ssols):
        print('Solution', idx+1, ':')
        print(sol)

is

::

    Solution 1 :
    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x : -3.23606797749979E+03   5.89091215835726E-87
     y :  1.47272803958932E-87   1.27201964951407E+03
    == err :  1.475E-16 = rco :  2.268E-01 = res :  6.661E-16 =
    Solution 2 :
    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x : -3.23606797749979E+03   5.89091215835726E-87
     y : -1.47272803958932E-87  -1.27201964951407E+03
    == err :  1.475E-16 = rco :  2.268E-01 = res :  6.661E-16 =
    Solution 3 :
    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x :  1.23606797749979E+03  -5.85209544336522E-95
     y :  7.86151377757424E+02   0.00000000000000E+00
    == err :  8.868E-17 = rco :  4.601E-01 = res :  1.110E-16 =
    Solution 4 :
    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x :  1.23606797749979E+03  -5.85209544336522E-95
     y : -7.86151377757424E+02   0.00000000000000E+00
    == err :  8.868E-17 = rco :  4.601E-01 = res :  1.110E-16 =

Let us look at the sum of the backward errors.
Executing

::

    verify(p, ssols)

yields the error 

::

    4.4853010194856324e-14

Observe that the error is about one hundred times 
smaller than without scaling."