Existence Verification for Singular Zeros of Nonlinear

Systems

by R. Baker Kearfott, rbk@usl.edu

Department of Mathematics, University of

Southwestern Louisiana and

      \underline{Jianwei Dian}, dian@usl.edu

Department of Mathematics, University of Southwestern Louisiana

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.

Singularities

When the Jacobi matrix F′(x^{∗}) is singular, computations as above cannot possibly prove existence and uniqueness. For such systems, the best that a preconditioner can do is reduce the Jacobi matrix to approximately the form

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

The Topological Degree

Uses : 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.

Uses : 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.

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 \underline{zeros of F_{¬k} on ∂x} with positive orientation at which sgn(f_{k}) = s, minus the number of \underline{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.

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}.

Structure of the System

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.

Structure of the System

Consequences

F′(x̌) is approximately the form

( \begin{matrix}1 & 0 & … & 0 & ∗ \\ 0 & 1 & 0 … & 0 & ∗ \\ ⋮ & ⋮ & ... & ⋮ & ⋮ \\ 0 & … & 0 & 1 & ∗ \\ 0 & … & 0 & 0 & 0 \end{matrix} ) .

So,

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

for 1 ≤ k ≤ n − 1.

f_{n}(x) = \frac{1}{2} ∑_{k,l=1}^{n} \frac{∂^{2}f_{n}}{∂x_{k}∂x_{l}}(x̌)(x_{k} − x̌_{k})(x_{l} − x̌_{l})
+O (∥x − x̌∥^{3})

singular system verification May, 1999 SIAM–11

Structure of the System

Consequences (continued) 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}) ,

So, for 1 ≤ k ≤ (n − 1),

u_{k}(x, y) ≈ (x_{k} − x̌_{k}) + \frac{∂f_{k}}{∂x_{n}}(x̌)(x_{n} − x̌_{n})
v_{k}(x, y) ≈ y_{k} + \frac{∂f_{k}}{∂x_{n}}(x̌)y_{n}

Structure of the System

Consequences (continued)

The following things are useful in the search phase (to be explained shortly) of the topological degree approach.

1. 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.

2. 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.

Computation of d(F̃, z, 0)

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

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

z = (x_{1}, y_{1}, …, x_{n}, y_{n})

2. The center of x_{k} is x̌_{k}.

The center of y_{k} is 0.

3. 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}.

4. Consider

\begin{aligned}F̃_{¬u_{n}}(x, y) & = & (u_{1}(x, y), v_{1}(x, y), …, \\ & & u_{n−1}(x, y), v_{n−1}(x, y), v_{n}(x, y)) \end{aligned}

on the boundary of (x, y).

For 1 ≤ k ≤ (n − 1), on x\underline{k}, F̃_{¬u_{n}}(x, y) = 0

⟹ u_{k}(x, y) ≈ (\underline{x}_{k} − x̌_{k}) + \frac{∂f_{k}}{∂x_{n}}(x̌)(x_{n} − x̌_{n}) = 0
⟹ \frac{|\underline{x}_{k} − x̌_{k}|}{|∂f_{k}/∂x_{n}(x̌)|} = |x_{n} − x̌_{n}|
⟹ \frac{w(x_{k})}{|∂f_{k}/∂x_{n}(x̌)|} ≤ w(x_{n})

Similarly, on x\overline{k},

F̃_{¬u_{n}}(x, y) = 0 ⟹ \frac{w(x_{k})}{|∂f_{k}/∂x_{n}(x̌)|} ≤ w(x_{n})

singular system verification May, 1999 SIAM–15

For 1 ≤ k ≤ (n − 1), on y\underline{k}, F̃_{¬u_{n}}(x, y) = 0

⟹ u_{k}(x, y) ≈ \underline{y}_{k} + \frac{∂f_{k}}{∂x_{n}}(x̌)y_{n} = 0
⟹ \frac{|\underline{y}_{k}|}{|∂f_{k}/∂x_{n}(x̌)|} = |y_{n}|
⟹ \frac{w(y_{k})}{|∂f_{k}/∂x_{n}(x̌)|} ≤ w(y_{n})

Similarly, on y\overline{k},

F̃_{¬u_{n}}(x, y) = 0 ⟹ \frac{w(y_{k})}{|∂f_{k}/∂x_{n}(x̌)|} ≤ w(y_{n})

singular system verification May, 1999 SIAM–16

1. Thus, if x_{n} is chosen so that

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}.

2. 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}.

3. So, we can eliminate 4n − 4 of the 4n faces of the boundary of (x, y), since we have arranged to verify F̃_{¬u_{n}}(x, y) ≠ 0 on each of these faces. Thus, we only need to search the four faces x\underline{n}, x\overline{n}, y\underline{n} and y\overline{n}, regardless of how large n is.

singular system verification May, 1999 SIAM–17

Some Hints

If

1. the approximations

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

for 1 ≤ k ≤ n − 1.

f_{n}(x) ≈ \frac{1}{2} ∑_{k,l=1}^{n} \frac{∂^{2}f_{n}}{∂x_{k}∂x_{l}}(x̌)(x_{k} − x̌_{k})(x_{l} − x̌_{l})

are exact; and 2.

∑_{k,l=1}^{n} \frac{∂^{2}f_{n}}{∂x_{k}∂x_{l}}(x̌)α_{k}α_{l} ≠ 0,

where α_{k} = ∂f_{k}/∂x_{n}(x̌) for 1 ≤ k ≤ n − 1, and α_{n} = −1,

then, d(F̃, z, 0) = 2.

singular system verification May, 1999 SIAM–18

Some Hints

Also, under the above assumptions,

• F̃_{¬u_{n}}(x, y) = 0 has no solutions on each of the 4n − 4 faces x\underline{k}, x\overline{k}, y\underline{k} and y\overline{k}, 1 ≤ k ≤ n − 1.

• F̃_{¬u_{n}}(x, y) = 0 has a unique solution on each of the 4 faces x\underline{n}, x\overline{n}, y\underline{n} and y\overline{n}.

We also know where the solution is, for example, on x\underline{k} the solution is (x̌_{1}, 0, x̌_{2}, 0, …, x̌_{n−1}, 0, \underline{x}_{n}, 0). This knowledge helps in the search phase which handles the faces x\underline{n}, x\overline{n}, y\underline{n} and y\overline{n} when the approximations are not exact.

The Actual Algorithm

Construction of the Box z

1. For 1 ≤ k ≤ n − 1, x_{k} is chosen to be centered at x̌_{k}. We need to take into consideration the accuracy of the approximate solver which finds the approximate solution of F (x) = 0.

2. For 1 ≤ k ≤ n − 1, y_{k} is chosen to be centered at 0. The width of y_{k} is chosen to be small, with some freedom.

3. x_{n} is chosen so that

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

4. 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̌)|}\} .

singular system verification May, 1999 SIAM–20

Elimination Phase

For k = 1 to n − 1

1. Use 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. This implies F̃_{¬u_{n}}(x, y) = 0 has no solutions on x\underline{k} and x\overline{k}

2. Use 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. This implies F̃_{¬u_{n}}(x, y) = 0 has no solutions on y\underline{k} and y\overline{k}

The Actual Algorithm

Search Phase

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

(a) Narrow x_{k}: 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.

Search Phase (continued)

(b) Search y_{n}:

i.A. A small subinterval y_{n}^{0} of y_{n} is

chosen centered at 0. B. The mean-value extensions for

v_{k}(x, y) = 0 are used to solve for y_{k} with width O (max(∥(x − x̌, y)∥^{2}, ∥y_{n}^{0}∥)), 1 ≤ k ≤ n − 1.

C. An interval Newton method can be set up for F̃_{¬u_{n}} to verify existence and uniqueness of a zero in y_{n}^{0}.

ii. y_{n}^{0} is inflated as much as possible provided the existence and uniqueness of the zero can be verified.

iii. Verify F̃_{¬u_{n}} = 0 has no solutions in the rest of y_{n}

The Actual Algorithm

Search Phase (continued)

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

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

Search Phase (continued)

(b) Search x_{n}:

i.A. A small subinterval x_{n}^{0} of x_{n} is

chosen centered at x̌_{n}. B. The mean-value extensions for

u_{k}(x, y) = 0 are used to solve for x_{k} with width O (max(∥(x − x̌, y)∥^{2}, ∥x_{n}^{0}∥)), 1 ≤ k ≤ n − 1.

C. An interval Newton method can be set up for F̃_{¬u_{n}} to verify existence and uniqueness of a zero in x_{n}^{0}.

ii. x_{n}^{0} is inflated as much as possible provided the existence and uniqueness of the zero can be verified.

iii. Verify F̃_{¬u_{n}} = 0 has no solutions in the rest of y_{n}

Search Phase (continued)

3. For each solution to F̃_{¬u_{n}} = 0 found in Steps 1b and 2b, compute the sign of u_{n}. Eliminate the solutions with negative u_{n}.

4. For each solution to F̃_{¬u_{n}} = 0 left in Steps 3, compute the sign of the determinant to get an orientation. Sum to get the degree.

Preliminary Experimental Results

1. Elimination phase: for all experiments that had been done, it was always successful to eliminate the 4n − 4 faces x\underline{k}, x\overline{k}, y\underline{k} and y\overline{k}, 1 ≤ k ≤ n − 1.

2. Search phase: for some of the experiments, it successfully searched the four faces x\underline{n}, x\overline{n}, y\underline{n} and y\overline{n}, and verified the degree was 2. For others, different problems occurred, indicating the need to adjust some parameters.