Existence and Uniqueness Verification for Singular Zeros

of Nonlinear Systems

by R. Baker Kearfott

Department of Mathematics, University of

Southwestern Louisiana

            rbk@usl.edu

and Jianwei Dian

Department of Mathematics, University of

Southwestern Louisiana and Arnold Neumaier

Institute of Mathematics, University of Vienna

The General Question

Given F : x → ℝ^{n} and x ∈ 𝕀ℝ^{n}, rigorously verify:

• there exists a unique x^{∗} ∈ x such that F (x^{∗}) = 0,

(1)

Computer arithmetic can be used to verify the assertion in Problem (1), with the aid of interval extensions and computational fixed point theorems.

The General Question

Uses

• Producing rigorous bounds on approximate solutions to linear and nonlinear systems (The approximate solutions can be computed with traditional techniques.)

– in analysis of stability of structures, where one wants to prove that all eigenvalues have negative real parts

– in robust computational geometry (surface intersection problems, etc.)

• As a tool in branch and bound algorithms in global optimization.

• As a tool in the verification that all zeros of a nonlinear system have been found in a region of ℝ^{n}.

The Nonsingular Case

Traditional Interval Newton Methods

Assumptions (roughly stated):

1. The Jacobi matrix F′(x^{∗}) is nonsingular.

2. x^{∗} is near the center of x.

3. The component widths of x are small.

4. N(F; x, x̌) is the image of x under an appropriate, preconditioned interval Newton method, with x̌ the center of x.

Then:

1. The preconditioned F′(x) is approximately the identity matrix.

2. Thus, N(F; x, x̌) ⊂ x. This proves that there is a unique solution of F (x) = 0 in x.

singular system verification SCAN’98–4

An Example

Example 1 Take

f_{1}(x_{1}, x_{2}) = x_{1}^{2} − x_{2},
f_{2}(x_{1}, x_{2}) = x_{1} − x_{2}^{2},

and x = (x_{1}, x_{2})^{T} = ([−0.1, 0.1], [−0.1, 0.3])^{T}.

Take x̌ = (0, 0.1)^{T}, so F (x̌) = (−0.1, −0.01)^{T}, and an interval extension of the Jacobi matrix is

F′(x) = ( \begin{matrix}2x_{1} & −1 \\ 1 & −2x_{2} \end{matrix} ) = ( \begin{matrix}[−.2, .2] & −1 \\ 1 & [−0.6, 0.2] \end{matrix} ) .

Precondition by the inverse of the midpoint matrix

Y = \{m(F′(X))\}^{−1} = ( \begin{matrix}−0.2 & 1 \\ −1 & 0 \end{matrix} ) ,

An Example (continued)

so the corresponding linear interval system is

Y F′(x)(x − x̌) = −Y F (x̌),

i.e.,

( \begin{matrix}[0.96, 1.04] & [−0.4, 0.4] \\ [−0.2, 0.2] & 1 \end{matrix} ) (x−x̌) = ( \begin{matrix}−0.01 \\ −0.1 \end{matrix} ) .

The interval Gauss–Seidel method applied to this system proves a unique solution in x:

x̃_{1} = 0 − \frac{−0.01 + [−0.4, 0.4][−0.2, 0.2]}{[0.96, 1.04]}
⊆ [−0.094, 0.073] ⊂ int([−0.1, 0.1]).

Similarly, x̃_{2} ⊂ int(x_{2}).

Singularities

When the Jacobi matrix F′(x^{∗}) is singular, computations as above cannot possibly prove existence and uniqueness. Example 2 Take

f_{1}(x_{1}, x_{2}) = x_{1}^{2} − x_{2},
f_{2}(x_{1}, x_{2}) = x_{1}^{2} + x_{2},

and x = (x_{1}, x_{2})^{T} = ([−0.1, 0.1], [−0.1, 0.3])^{T}. For such systems, the best that a preconditioner can do is reduce the Jacobi matrix to approximately the form

( \begin{matrix}∗ & 0 & … & 0 & \overbrace{∗ … ∗}^{n - rank} \\ 0 & ∗ & 0 … & 0 & ∗ … ∗ \\ ⋮ & ⋮ & ... & ⋮ & ⋮ \\ 0 & … & 0 & ∗ & ∗ … ∗ \\ 0 & … & 0 & 0 & 0 … 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & … & 0 & 0 & 0 … 0 \end{matrix} ) .

Singularities

Philosophical Considerations

Uniqueness verification within the original x is not possible. Alternatives are:

1. (Easy but perhaps not always adequate) Verify the system has an ϵ–approximate solution within x.

2. (Only possible in special cases) Verify the system has at least one solution in x.

3. (More difficult computationally) Verify the system has an exact number of solutions, counting multiplicities, in a complex extension of x.

Which alternative is appropriate in particular contexts?

ϵ–Approximate Solution

1. Simply evaluate the interval extension F(x).

2. If ∥F(x)∥ < ϵ, then each component of F is less than ϵ everywhere within x. (But all components of f may not simultaneously vanish.)

3. This technique can be used, along with scaling, etc., to handle non-isolated solutions in branch and bound algorithms.

Verification of at Least One Solution

1. The topological degree (to be explained shortly) may be computed over x.

2. If the topological degree is non-zero, there is at least one solution of F (x) = 0 in x.

3. No conclusion can be reached if the topological degree is zero.

Verification of the Exact Multiplicity

1. If F : ℂ^{n} → ℂ^{n}, then the topological degree of F over x gives the exact number of solutions, counting multiplicities.

2. If F : ℝ^{n} → ℝ^{n}, and F can be extended analytically into ℂ^{n}, then computations can verify existence of an exact solution or solutions (with multiplicity computed by the algorithm) within a small region of complex space containing x.

3. Sometimes the value of the topological degree (i.e. the multiplicity of the solution) is of interest beyond the existence/uniqueness question.

The Topological Degree

What is it?

• If F : x ⊂ ℝ^{n} → ℝ^{n}, F′(x^{∗}) ≠ 0 wherever F (x^{∗}) = 0, x^{∗} ∈ x, and F (x) ≠ 0 when x ∈ ∂x, then the degree d(F, x, 0) is the number of x^{∗} ∈ x, F (x^{∗}) = 0 with det(F′(x^{∗})) > 0, minus the number of such x^{∗} ∈ x with det(F′(x^{∗})) < 0.

• d(F, x, 0) is a continuous function of F, and is defined even if det(F′(x^{∗})) = 0, as long as there are no solutions to F (x) = 0 on ∂x.

• If F is extended to ℂ^{n} and is thought of as mapping ℝ^{2n} to ℝ^{2n}, and x is embedded in a box z ∈ ℂ^{2n}, then d(F, z, 0) is equal to the exact number of z ∈ z, F (z) = 0, counting multiplicities.

singular system verification SCAN’98–12

An Example

• If

f_{1}(x, y) = x^{2} − y^{2} − ϵ^{2}
f_{2}(x, y) = 2xy,

If ϵ ≠ 0, then F has solutions at (x, y) = (ϵ, 0) and (x, y) = (−ϵ, 0). Since det(F′(x)) = 4(x^{2} + y^{2}) = 4ϵ^{2} at each of these solutions, d(F, z, 0) = 2, where

z = {(x, y) | x ∈ [−0.1, 0.1], y ∈ [−δ, δ]}

for any δ > 0.

• If ϵ = 0, then d(F, z, 0) is still equal to 2, even though the Jacobi matrix vanishes at the only solution (x, y) = (0, 0).

How is it Computed?

• d(F, x, 0) depends only on values of F on ∂x.

• Define

F_{¬k}(x) = (f_{1}(x), …, f_{k−1}(x),
f_{k+1}(x), …, f_{n}(x)),

and select s ∈ {−1, 1}. Then d(F, x, 0) is equal to the number of zeros of F_{¬k} on ∂x with positive orientation at which sgn(f_{k}) = s, minus the number of zeros of F_{¬k} on ∂x with negative orientation at which sgn(f_{k}) = s.

• The orientation is computed by computing the sign of the determinant of the Jacobian of F_{¬k} and by taking account of which face.

singular system verification SCAN’98–14

Computation of the Degree

[figure: f_{1}(x) = 0; f_{2} > 0; f_{2} > 0; f_{1} = x_{1}^{2} − x_{2}^{2}; f_{2} = 2x_{1}x_{2}]

Computational Cost

1. Directly finding all zeros of F_{¬k} on ∂x can be done in a straightforward branch and bound algorithm. However, that is perhaps too expensive for mere verification purposes.

2. The structure of the preconditioned system can be used to greatly simplify the computations.

3. The widths of the box x constructed about the approximate solution can be chosen so that only several one-dimensional searches need be done to compute d(F, z, 0), where F : ℂ^{n} → ℂ^{n}.

Notation and Assumptions

• For F : ℝ^{n} → ℝ^{n}, extend F to complex space: z = x + iy, u_{k}(x, y) = ℜ(f_{k}(z)) and v_{k}(x, y) = ℑ(f_{k}(z)).

• Define F̃(x, y) = (u_{1}(x, y), v_{1}(x, y), …, u_{n}(x, y), v_{n}(x, y)) : ℝ^{2n} → ℝ^{2n}.

• Assume F (x^{∗}) ≈ 0.

• Assume F has been preconditioned (say, through an incomplete LU factorization). Also assume F′(x^{∗}) has null space of dimension 1, so

f_{k}(x) ≈ (x_{k} − x^{∗}_{k}) + \frac{∂f_{k}}{∂x_{n}}(x^{∗})(x_{n} − x^{∗}_{n})

for 1 ≤ k ≤ n − 1.

Structure of the System

One-Dimensional Null Space For 1 ≤ k ≤ (n − 1),

u_{k}(x, y) = (x_{k} − x̌_{k}) + \frac{∂f_{k}}{∂x_{n}}(x̌)(x_{n} − x̌_{n})
+O (∥(x − x̌, y)∥^{2})
v_{k}(x, y) = y_{k} + \frac{∂f_{k}}{∂x_{n}}(x̌)y_{n}
+O (∥(x − x̌, y)∥^{2}) ,

and

u_{n}(x, y) = O (∥(x − x̌, y)∥^{2}) ,
v_{n}(x, y) = O (∥(x − x̌, y)∥^{2}) .

Structure of the System

Consequences

1. Mean-value interval extensions u_{k} and v_{k} can be formed, 1 ≤ k ≤ n − 1.

2. If x_{n} is known precisely, formally solving u_{k}(x, y) = 0 for x_{k} gives x_{k} with w(x_{k}) = O (∥(x − x̌, y)∥^{2}), 1 ≤ k ≤ n − 1.

3. If y_{n} is known precisely, formally solving v_{k}(x, y) = 0 for y_{k} gives y_{k} with w(y_{k}) = O (∥(x − x̌, y)∥^{2}), 1 ≤ k ≤ n − 1.

A Degree-Computation Algorithm

Construction of the Box z

1. Define x = ([\underline{x}_{1}, \overline{x}_{1}], …, [\underline{x}_{n}, \overline{x}_{n}]) and y = ([\underline{y}_{1}, \overline{y}_{1}], …, [\underline{y}_{n}, \overline{y}_{n}]).

2. Define x\underline{k} as (x, y) with [\underline{x}_{k}, \overline{x}_{k}] replaced by \underline{x}_{k}, and define x\overline{k} as (x, y) with [\underline{x}_{k}, \overline{x}_{k}] replaced by \overline{x}_{k}. Similarly define y\underline{k} and y\overline{k}.

3. u_{k}(x, y) = 0 on x\underline{k}, 1 ≤ k ≤ n − 1, at approximately

x_{n} = x̌_{n} + \frac{x̌_{k} − \underline{x}_{k}}{∂f_{k}/∂x_{n}(x̌)},

u_{k}(x, y) = 0 on x\overline{k}, 1 ≤ k ≤ n − 1, at approximately

x_{n} = x̌_{n} + \frac{x̌_{k} − \overline{x}_{k}}{∂f_{k}/∂x_{n}(x̌)}.

singular system verification SCAN’98–20

A Degree-Computation Algorithm

Construction of z (continued)

4. Similarly, v_{k}(x, y) = 0 on y\underline{k}, 1 ≤ k ≤ n − 1, at approximately y_{n} = \frac{−\underline{y}_{k}}{∂f_{k}/∂x_{n}(x̌)}, and v_{k}(x, y) = 0 on y\overline{k}, 1 ≤ k ≤ n − 1, at approximately y_{n} = \frac{−\overline{y}_{k}}{∂f_{k}/∂x_{n}(x̌)}.

5. Thus, if x_{n} is chosen so

w(x_{n}) ≤ \frac{1}{2} min_{1≤k≤n−1} \{\frac{w(x_{k})}{|∂f_{k}/∂x_{n}(x̌)|}\} ,

then it is unlikely that u_{k}(x, y) = 0 on either x\underline{k} or x\overline{k}.

6. Similarly, if y_{n} is chosen so that

w(y_{n}) ≤ \frac{1}{2} min_{1≤k≤n−1} \{\frac{w(y_{k})}{|∂f_{k}/∂x_{n}(x̌)|}\} ,

then it is unlikely that v_{k}(x, y) = 0 on either y\underline{k} or y\overline{k}.

singular system verification SCAN’98–21

Degree Computation

An Actual Algorithm

1. For k = 1 to n − 1:

(a) Do mean-value interval evaluations of u_{k}(x, y) over x\underline{k} and x\overline{k} to show that u_{k}(x, y) ≠ 0 on these faces of z.

(b) Similarly do second-order interval evaluations of v_{k}(x, y) over y\underline{k} and y\overline{k} to show that v_{k}(x, y) ≠ 0 on these faces of z.

2. On x\underline{n} and x\overline{n}:

(a) Use mean-value extensions u_{k}(x, y) = 0 to solve for x_{k} with width O (∥(x − x̌, y)∥^{2}), 1 ≤ k ≤ n − 1.

(b) Perform a binary search on y_{n} to find verified intervals where F̃_{¬v_{n}} = 0 and v_{n}(x, y) > 0.

singular system verification SCAN’98–22

A Degree Computation Algorithm

(continued)

3. On y\underline{n} and y\overline{n} (similar to Step 2):

(a) Use mean-value extensions v_{k}(x, y) = 0 to solve for x_{k} with width O (∥(x − x̌, y)∥^{2}), 1 ≤ k ≤ n − 1.

(b) Perform a binary search on x_{n} to find verified intervals where F̃_{¬v_{n}} = 0 and v_{n}(x, y) > 0.

4. For each solution to F_{¬v_{n}} = 0 found in Steps 2b and 3b, compute an orientation, to sum to find the degree.

The Degree Computation Algorithm

Some Notes

• For each interval on y_{n} produced in the search in Step 2b, narrower intervals on y_{k} can be produced with a mean-value extension v_{k}(x) = 0, 1 ≤ k ≤ n − 1.

• Similarly, for each interval on x_{n} produced in the search in Step 3b, narrower intervals on x_{k} can be produced with a mean-value extension u_{k}(x) = 0, 1 ≤ k ≤ n − 1.

• An interval Newton method can be set up for F̃_{¬v_{n}} to verify existence and uniqueness.