An intrinsic homotopy for intersecting algebraic varieties

Andrew J. Sommese^{∗}  Jan Verschelde^{†}

Charles W. Wampler^{‡}

September 2, 2004

Abstract

Recently we developed a diagonal homotopy method to compute a numerical representation of all positive dimensional components in the intersection of two irreducible algebraic sets. In this paper, we rewrite this diagonal homotopy in intrinsic coordinates, which reduces the number of variables, typically in half. This has the potential to save a significant amount of computation, especially in the iterative solving portion of the homotopy path tracker. Three numerical experiments all show a speedup of about a factor two. 2000 Mathematics Subject Classification. Primary 65H10; Secondary 13P05, 14Q99, 68W30. Key words and phrases. Components of solutions, embedding, generic points, homotopy continuation, irreducible components, numerical algebraic geometry, polynomial system.

∗Department of Mathematics, University of Notre Dame, Notre Dame, IN 46556-4618, U.S.A.; Email: sommese@nd.edu URL: http://www.nd.edu/˜ sommese. This material is based upon work supported by the National Science Foundation under Grant No. 0105653; and the Duncan Chair of the University of Notre Dame.

†Department of Mathematics, Statistics, and Computer Science, University of Illinois at Chicago, 851 South Morgan (M/C 249), Chicago, IL 60607-7045, U.S.A.; Email: jan@math.uic.edu or jan.verschelde@na-net.ornl.gov URL: http://www.math.uic.edu/˜ jan. This material is based upon work supported by the National Science Foundation under Grant No. 0105739 and Grant No. 0134611.

‡General Motors Research and Development, Mail Code 480-106-359, 30500 Mound Road, Warren, MI 48090-9055, U.S.A.; Email: Charles.W.Wampler@gm.com.

Our goal is to compute the irreducible decomposition of A ∩ B ⊂ ℂ^{k}, where A and B are irreducible algebraic sets. In particular, suppose that

• A is an irreducible component of the solution set of a polynomial system f_{A}(u) = 0 defined on ℂ^{k}, and similarly

• B is an irreducible component of the solution set of a polynomial system f_{B}(u) = 0 defined on ℂ^{k}.

This includes the important special case when f_{A} and f_{B} are the same system, but A and B are distinct irreducible components.

Casting this problem into the framework of numerical algebraic geometry, we assume that all components are represented as witness sets. For an irreducible component A ⊂ ℂ^{k} of dimension dim(A) and degree deg(A), a witness set consists of a generic k − dim(A) dimensional linear subspace L ⊂ ℂ^{k} and the deg(A) points of intersection A ∩ L. We assume that at the outset we are given such sets for A and B, and our goal is to compute witness sets for the irreducible components of A ∩ B. The intersection may break into several such components, and the components may have various dimensions. Our methods proceed in two phases: we first find a witness superset guaranteed to contain witness points for all the components, then we break this set into its irreducible components. We recently reported on an algorithm [15], herein called the extrinsic^{1} homotopy method, for computing a witness superset for A ∩ B. This can then be decomposed into irreducible components using the methods in [14] and its references.

Abstracting away the details, which are discussed more fully in §1, the extrinsic method consists of a cascade of homotopies in unknowns x ∈ ℂ^{N} and path parameter t ∈ [0, 1], each of the form

H(x, t) := [ \begin{matrix}f(x) \\ t(P x + p) + (1 − t)(Qx + q) \end{matrix} ] = 0  (1)

where f : ℂ^{N} → ℂ^{m} is a system of polynomial equations, P, Q are (N − m) × N full-rank matrices, and p, q ∈ ℂ^{(N−m)} are column vectors. There is a homotopy of this form for each dimension where A ∩ B could have one or more solution components. We know solution values for x at t = 1 and wish to track solution paths x(t) implicitly defined by (1) as t → 0 to get x(0).

1The terminology extrinsic/intrinsic is in analogy with the homotopies of [4].

At any specific value of t, this looks like

Ĥ(x, t) = [ \begin{matrix}f(x) \\ R(t)x + r(t) \end{matrix} ] = 0,  (2)

where R = tP + (1 − t)Q and r = tp + (1 − t)q. The homotopy is constructed such that we are assured that R(t) is full rank for all t ∈ [0, 1]. Thus, the linear subspace of solutions of R(t)x + r(t) = 0 can be parameterized by u ∈ ℂ^{m} in the form

x(u, t) = R^{⊥}(t)u + x_{p}(t),  (3)

where x_{p}(t) is any particular solution and R^{⊥}(t) is the right null space of R(t), that is, R^{⊥} is a full-rank N × m matrix with RR^{⊥} = 0. We may restrict Ĥ to this linear subspace to obtain

H̃(u, t) := Ĥ(x(u, t), t) = f(R^{⊥}(t)u + x_{p}(t)) = 0,  (4)

where we have dropped the linear equations because by construction, they are identically zero for all t. We refer to this as the intrinsic form of the equations.

The problem with (4) is that it requires computing R^{⊥} and x_{p} at each new value of t as we follow the homotopy paths. Because of this, H̃(x) offers little, if any, computational advantage over the extrinsic Ĥ(x).

Although not generally possible, for some P, Q, p, q, one can convert the extrinsic homotopy (1) into an intrinsic homotopy of the form

H̃(u, t) = f(t(Cu + c) + (1 − t)(Du + d)) = 0,  (5)

in which the path parameter t appears linearly. This means that the linear algebra to compute C, D ∈ ℂ^{N×m} and c, d ∈ ℂ^{N} is done just once at the outset, rather than being repeated at each value of t. This can save a significant amount of computation and is also simpler to implement.

This paper is organized as follows. In §1, we review the extrinsic homotopies formulated in [15] for intersecting algebraic varieties, and in §2.1 and §2.2, we show how to convert these to the linear intrinsic form. A comparison of the numerical behavior of the extrinsic homotopies and intrinsic homotopies is presented in §3. Acknowledgements. The authors thank the referees for their careful reading and suggestions to improve the paper.

1 Extrinsic Diagonal Homotopies

Let A ⊂ ℂ^{k} and B ⊂ ℂ^{k} be as in the opening paragraph, having dimensions a and b respectively. We have bounds on the dimension of components of A ∩ B as follows. After renaming if necessary, we may assume a ≥ b. The largest possible dimension of A ∩ B is therefore b, which happens if and only if B is contained in A. We can check this possibility using a homotopy membership [12] test to see if a generic point of B is in A. If so, we have A ∩ B = B and no further computation is needed. Otherwise, we know that the largest possible dimension of A ∩ B is b − 1. On the other hand, because the codimension of A ∩ B is at most the sum of the codimensions of A and B, the smallest possible dimension of any component of A ∩ B is max(a + b − k, 0). For a particular problem, one might have available some tighter bounds on dim(A ∩ B), and if so, one can take advantage of that knowledge in the algorithm to follow. Accordingly, we introduce the symbols h^{∗} and h_{0} as follows:

b ≥ h^{∗} > dim(A ∩ B),  (6)
max(a + b − k, 0) ≤ h_{0} ≤ min(dim(any component of A ∩ B)). (7)

Unless we have other knowledge, we use the defaults h^{∗} = b and h_{0} = max(a+ b − k, 0).

Instead of working directly in ℂ^{k}, we find the intersection A∩B by casting the problem into (u, v) ∈ ℂ^{k+k} and restricting to the diagonal u − v = 0. More precisely, the product X := A × B ⊂ ℂ^{k+k} is an affine variety of dimension a + b, i.e., an irreducible affine algebraic set of dimension a + b. The intersection of A and B can be identified, e.g., [2, Ex. 13.15] or [10, pg. 122ff], with X ∩ Δ where Δ is the diagonal of ℂ^{k+k} defined by the system

δ(u, v) := [ u_{1} − v_{1} ; ⋮ ; u_{k} − v_{k} ] = 0  (8)

with (u, v) giving the coordinates of ℂ^{k+k}.

The initial data consists of witness sets for A and B. That is, our data for A consists of a generic system L_{A}(u) = 0 of a linear equations and the deg(A) solutions {α_{1}, …, α_{deg(A)}} ⊂ ℂ^{k} of the system

[ \begin{matrix}f_{A}(u) \\ L_{A}(u) \end{matrix} ] = 0,  (9)

and similarly the data for B consists of a generic system L_{B}(v) = 0 of b linear equations and the deg(B) solutions {β_{1}, …, β_{deg(B)}} ⊂ ℂ^{m} of the system

[ \begin{matrix}f_{B}(v) \\ L_{B}(v) \end{matrix} ] = 0.  (10)

Remark 1.1 We are not assuming that A and B occur with multiplicity one in the solution sets of their respective systems f_{A}(u) = 0 and f_{B}(v) = 0. If the multiplicity is greater than one, we must use a singular path tracker [13].

The extrinsic algorithm can be summarized concisely by introducing a bit of matrix notation. First, let

w = [ \begin{matrix}u \\ v \end{matrix} ] ∈ ℂ^{2k},  (11)

and introduce a column vector of “slack” variables z ∈ ℂ^{k}. Also, define the k × k projection matrix

P_{h} = diag (\underbrace{1, …, 1}_{h}, \underbrace{0, …, 0}_{k−h})  (12)

Left multiplication by P_{h} picks out the first h rows of its multiplicand and right multiplication picks out the first h columns of its multiplier. Note also that P_{h}^{2} = P_{h}. Similarly, let P_{ji} be the k × k matrix

P_{ji} = diag (\underbrace{0, …, 0}_{j}, \underbrace{1, …, 1}_{i−j}, \underbrace{0, …, 0}_{k−i}),  (13)

which picks out rows (or columns) j + 1, …, i. It is useful to note that P_{j} + P_{ji} = P_{i}.

The formulation of the homotopy requires several random matrices as follows. First, we choose generic matrices

M ∈ ℂ^{(k−a)×#(f_{A})},  N ∈ ℂ^{(k−b)×#(f_{B})},  (14)

where #(f_{A}) is the number of functions in the system f_{A}(x) associated to component A, and similarly for #(f_{B}). These are used to define

F(w) := [ \begin{matrix}Mf_{A}(u) \\ Nf_{B}(v) \end{matrix} ] .  (15)

Note that A×B is an irreducible component of the solution set of the system F(w) = 0. Next, we choose 𝔸 a generic (a + b) × k matrix, and let

A = [ 𝔸 −𝔸 ] ∈ ℂ^{(a+b)×2k}  (16)

so Aw = 𝔸(u − v). Finally, we choose generic matrices

B ∈ ℂ^{(a+b)×k},  C ∈ ℂ^{k×2k}  d ∈ ℂ^{k×1}.  (17)

In all these, a matrix with random complex elements will be generic with probability one.

Since the smallest dimensional nonempty component of A ∩ B is of dimension at least max{0, a + b − k}, it follows from [15, Lemma (3.1)] that we can find the irreducible decomposition of A ∩ B by finding the irreducible decomposition of Aw = 0 on X = A × B. For this purpose, we consider a cascade of homotopies of the form

E_{h}(w, z) = [ \begin{matrix}F(w) \\ Aw + BP_{h}z \\ z − P_{h}(Cw + d) \end{matrix} ] = 0,  (18)

which is well-defined for any integer 0 ≤ h ≤ k. Denoting the entries of z as z_{1}, …, z_{k}, note that the last row of this matrix equation implies that (z_{h+1}, …, z_{k}) = 0. The method for generating a witness superset consists of solving E_{h^{∗}}(w, z) = 0 and then descending sequentially down the cascade to solve E_{j}(w, z) = 0 for j = h^{∗} − 1, …, h_{0}.

The rationale behind the cascade is that the linear system P_{h}(Cw +d) = 0 is a linear slice that cuts out witness points for solution components of dimension h. The vector z is a set of slack variables. A solution point of E_{h}(w, z) = 0 for which z = 0 is on the slice and thus gives a witness point. Solution points with z ≠ 0 are not on the slice, and we call these “nonsolutions.” These become the starting points for the next step of the cascade. (We state this more formally below, after giving more details of the algorithm.) For each step down the cascade, one more slack variable is set to zero and a corresponding hyperplane is removed from the slice. The recycling of nonsolutions as starting points for the next step of the cascade is valid due to the fact that for j < i, E_{j}(w, z) is just E_{i}(w, z) with certain elements of B, C, and d set to zero. This is justified in [15].

The following steps of the algorithm still need to be described:

• how to solve E_{h^{∗}}(w, z) = 0,

• how to descend the cascade, and

• how to reap the witness points from the solutions at each level of the cascade.

The homotopy to solve E_{h^{∗}}(w, z) = 0 is

[ \begin{matrix}F(w) \\ (1 − t) [ \begin{matrix}Aw + BP_{h^{∗}}z \\ z − P_{h^{∗}}(Cw + d) \end{matrix} ] + tγ [ \begin{matrix}L_{A}(u) \\ L_{B}(v) \\ z \end{matrix} ] \end{matrix} ] = 0,  (19)

where γ is a random complex number. At t = 1, solution paths start at the deg(A)×deg(B) nonsingular solutions {(α_{1}, β_{1}), …, (α_{deg(A)}, β_{deg(B)})} ⊂ ℂ^{2k} obtained by combining the witness points for A and B. At t = 0, the solution paths terminate at the desired start solutions for E_{h^{∗}}(w, z) = 0. In [15] we ended the homotopy at E_{b}(w, z) = 0, but the argument works equally well with h^{∗} in place of b.

The homotopy connecting E_{i} to E_{j} for j < i is

H_{i,j}(τ, w, z) := [ \begin{matrix}F(w) \\ Aw + BP_{i}z \\ z − (P_{j} + τP_{ji})(Cw + d) \end{matrix} ] = 0,  (20)

where τ goes from 1 to 0 along a sufficiently general 1-real-dimensional curve. For example, for all but finitely many γ ∈ ℂ of absolute value 1, τ = r + γr(1−r) as r goes from 1 to 0 on the real interval suffices. Another possibility, relevant in what comes below, is

τ = t/(t + γ(1 − t))  (21)

as t goes from 1 to 0 on the real interval.

In the cascade of homotopies from [15] (based on [11]), we start out with the finite set G_{i} of nonsingular solutions of E_{i} with z_{i} ≠ 0. Tracking these start solutions we end up with a set of solutions G_{i,j}^{E} of E_{j} with z_{h} = 0 for h > j. In [15], j = i − 1, but the argument there works immediately for any j < i. The key points about the set G_{i,j}^{E} is that 1. the set G_{j} equals the set of points in G_{i,j}^{E} for which z_{j} ≠ 0;

2. the set of points Ŵ_{j} ⊂ G_{i,j}^{E} for which z_{h} = 0 for all h ≤ j contains a witness point set W_{j} for the j-dimensional components of the solution set of the intersection of A and B.

We also know that the set of points in G_{i,j}^{E} for which z_{h} = 0 for all h ≤ j equals the set of points in G_{i,j}^{E} for which z_{j} = 0. We wish to set up an intrinsic homotopy such that analogs of the above key facts hold true.

2 Setting Up Intrinsic Homotopies

The extrinsic homotopies of (19) and (20) use the variables (w, z) ∈ ℂ^{2k} ×ℂ^{k}. Each has a + b + k linear equations which we wish to eliminate by converting to an intrinsic homotopy. The result will be homotopies in intrinsic variables y ∈ ℂ^{2k−a−b}. Note that 2k − (a + b) is the codimension of A × B in ℂ^{2k}. It is also the sum ā + b̄ of the codimension ā = k − a of A in ℂ^{k} and the codimension b̄ = k − b of B in ℂ^{k}. Since this quantity appears frequently in the expressions below, we define

m = 2k − a − b.  (22)

Accordingly, our intrinsic homotopy variables are y ∈ ℂ^{m}.

2.1 Intrinsic Start Homotopy

In this section, we replace the extrinsic start homotopy of (19) with one having the intrinsic form of (5). Fixing a particular solution

w_{1} = [ \begin{matrix}u_{p} \\ v_{p} \end{matrix} ]  (23)

of

[ \begin{matrix}L_{A}(u) \\ L_{B}(v) \end{matrix} ] = 0,  (24)

choose a basis W_{1} ∈ ℂ^{2k×m} of the null space N_{1} of

[ \begin{matrix}L_{A}(u) − L_{A}(0) \\ L_{B}(v) − L_{B}(0) \end{matrix} ] = 0.  (25)

The solutions (α_{i}, β_{j}) of (24) arising from (9) and (10) correspond to (w_{1} + N_{1}) ∩ (A × B).

Fixing a particular solution w_{2} of

Aw + BP_{h^{∗}}(Cw + d) = 0,  (26)

choose a basis W_{2} ∈ ℂ^{2k×m} of the null space N_{2} of

Aw + BP_{h^{∗}}Cw = 0.  (27)

We have the intrinsic homotopy with variable y ∈ ℂ^{m}

F (τ [w_{1} + W_{1}y] + (1 − τ) [w_{2} + W_{2}y]) = 0.  (28)

An irreducible analytic set A_{1} ⊂ ℂ^{2k} is said to be transverse to an irreducible analytic set A_{2} ⊂ ℂ^{2k} if

1. the intersection A_{1} ∩ A_{2} is contained in the set of smooth points of both A_{1} and A_{2}; and

2. at any point x ∈ A_{1} ∩ A_{2},

dim (T_{A_{1},x} ∩ T_{A_{2},x}) = dim T_{A_{1},x} + dim T_{A_{2},x} − 2k.

Note that this implies that A_{1}∩A_{2} is a manifold whose connected components have dimension dim A_{1} + dim A_{2} − 2k. Since w_{1} + N_{1} is transverse to A × B, the (2k − a − b)-dimensional affine subspace given by

{τ_{1} [w_{1} + W_{1}y] + τ_{2} [w_{2} + W_{2}y] | y ∈ ℂ^{m}}  (29)

is transverse to A × B for all but a finite set of [τ_{1}, τ_{2}] ∈ ℙ^{1}. In particular for all but a finite number of γ ∈ ℂ of absolute value one, with the relation between τ and t as in (21), the m-dimensional affine subspace given by

{τ [w_{1} + W_{1}y] + (1 − τ) [w_{2} + W_{2}y] | y ∈ ℂ^{m}}  (30)

is transverse to A × B for all t ∈ (0, 1]. By genericity in the choices of A, B, C, d, this is true for t = 0 also. Thus using the homotopy (28) to track the paths starting with the (α_{i}, β_{j}) at t = 1, we get the start solutions of the cascade at t = 0.

In practice it will be convenient to go directly from solutions (α_{i}, β_{j}) of (24) arising from (9) and (10) to E_{h^{∗}−1} or any E_{j} with j < h^{∗}. Doing this we want to know that the limits of the paths of the intrinsic homotopy starting with the solutions (α_{i}, β_{j}) contain the subset G_{j} for which z_{j} ≠ 0 and a set of points Ŵ_{j} which contains a set of witness points W_{j}. This is true for both the intrinsic and the earlier extrinsic homotopy of [15]. The reason why this is so is that the solutions G_{j} ∪Ŵ_{j} are contained in the set of isolated solutions of E_{j} restricted to A × B. Therefore by [15, Lemma A.1], there is a Zariski open set of t ∈ ℂ such that except for a finite choice of γ of absolute value one in (21), G_{j} ∪ Ŵ_{j} are limits of isolated solutions of the homotopy (28) restricted to A×B. Since the solutions at t = 1 of the homotopy (28) on A×B are the transversal intersection with the m-dimensional affine subspace given by Eq.(30), it follows that for the t near 1 this is still true. Thus the isolated solutions of the homotopy (28) for a Zariski open set of the t are continuations from solutions (α_{i}, β_{j}) of (24) arising from (9) and (10), and in consequence G_{j} ∪ Ŵ_{j} are contained in limits of isolated solutions of the homotopy (28) restricted to A × B starting at these points.

The current default is to go directly from solutions (α_{i}, β_{j}) of (24) arising from (9) and (10) to E_{h^{∗}−1}.

2.2 Intrinsic Cascade Homotopies

In this section, we convert the extrinsic cascade homotopies of (20) into intrinsic homotopies of the form of (5). This must be done a bit more delicately than what was done for the start homotopy, because we must preserve the containment of H_{i,j} inside the parameter space of E_{i} so that we retain the properties stated at the end of §1. We do this by deriving an intrinsic homotopy whose path is exactly the same as a generic real path from τ = 1 to τ = 0 in (20).

We start by eliminating z by substitution from the last block row of (20) into the middle row. We use the facts that for i > j, P_{i}P_{j} = P_{j} and P_{i}P_{ji} = P_{ji} to obtain

H_{i,j}(t, w) := [ \begin{matrix}F(w) \\ Aw + B(P_{j} + τP_{ji})(Cw + d) \end{matrix} ] = 0,  (31)

which, abusing notation, we still call H_{i,j}. By similar abuse of notation, we use E_{h}(w) in place of E_{h}(w, z) after eliminating z from (18).

Our first observation concerns the existence of a constant particular solution throughout the cascade.

Lemma 2.1 The inhomogeneous linear system

[ \begin{matrix}I_{k} & −I_{k} \\ & C \end{matrix} ] w = [ \begin{matrix}0 \\ −d \end{matrix} ]  (32)

has a unique nonzero solution ϵ.

Proof. The genericity of C implies the invertibility of [ \begin{matrix}I_{k} & −I_{k} \\ & C \end{matrix} ]. □

Notice that this implies that both Aϵ = 0 and Cϵ + d = 0, and therefore w = ϵ is a solution of

Aw + B(P_{j} + τP_{ji})(Cw + d) = 0  (33)

for any i, j, τ.

Let Y_{h} be the homogeneous linear system

Y_{h} := (A + BP_{h}C)w = 0.  (34)

The following lemma concerning the null space of Y_{h} is crucial for the conversion to an intrinsic form.

Lemma 2.2 For any j and i such that h_{0} ≤ j < i ≤ h^{∗}, there exist matrices E_{i,j} ∈ ℂ^{2k×(m−i+j)} and F_{i,j}, G_{i,j} ∈ ℂ^{2k×(i−j)} such that

1. [E_{i,j} F_{i,j}] = Null Y_{i}

2. [E_{i,j} G_{i,j}] = Null Y_{j}

3. P_{ji}CF_{i,j} = P_{ji}CG_{i,j} = [ \begin{matrix}0 \\ I_{i−j} \\ 0 \end{matrix} ],

where the (i − j) × (i − j) identity matrix I_{i−j} appears in rows j + 1, …, i.

Proof. We must first establish that Y_{i} and Y_{j} are full row rank a + b so that m = 2k − a − b is the correct dimension of their null spaces. Since A depends on generic 𝔸 (see (16)) and B and C are generic, it suffices to show that there is at least one choice of 𝔸, B, C such that Y_{h} is full rank for h_{0} ≤ h ≤ h^{∗}. For a + b < k, it suffices to choose B = 0, C = 0 and choose 𝔸 to make Y_{h} = [I_{a+b} 0 − I_{a+b} 0]w. For a + b > k, choose 𝔸 = [I_{k} 0]^{T}, choose B with I_{a+b−k} in the lower left and C with I_{a+b−k} in the upper left. Since h ≥ a + b − k, this suffices to make Y_{h} full rank, as one may check by direct substitution.

Next, we establish that Y_{i} and Y_{j} share a null subspace of dimension m − i + j. Note that

Y_{i} = (A + B(P_{j} + P_{ji})C)w = (Y_{j} + BP_{ji}C)w.  (35)

The matrix BP_{ji}C is independent of BP_{j}C because the projection matrices pick out different rows and columns of generic matrices B and C. Accordingly, the subspace Null Y_{i} ∩Null Y_{j} = Null Y_{j} ∩Null (BP_{ji}C). These have dimension m and codimension (i − j), respectively, and they meet transversely, so the intersection has dimension m − i + j. Let E_{i,j} be any basis for this subspace.

Now, suppose F̂_{i,j} completes a basis [E_{i,j} F̂_{i,j}] for Null Y_{i}. It must be independent of Null (BP_{ji}C), and since B is generic, this implies that P_{ji}CF̂_{i,j} must be full rank. Since P_{ji} zeros out all but rows j + 1, …, i, this implies that

P_{ji}CF̂_{i,j} = [ \begin{matrix}0 \\ Q \\ 0 \end{matrix} ]  (36)

must have a full-rank (i − j) × (i − j) matrix Q in rows j + 1, …, i. Then, F_{i,j} = F̂_{i,j}Q^{−1} completes the basis of Y_{i} while also satisfying Condition 3 of the lemma. Similar reasoning shows the existence of G_{i,j}. □

Choosing a random γ ∈ ℂ, we form the linear system

W_{i,j}(t, y) = ϵ + [ E_{i,j} tF_{i,j} + γ(1 − t)G_{i,j} ] y  (37)

where y ∈ ℂ^{m} and ϵ is as in Lemma 2.1. From this, we form the intrinsic homotopy

H_{i,j}(t, y) = F (W_{i,j}(t, y)) = 0,  (38)

and track y as t goes from 1 to 0 on the real interval.

The crucial fact behind the equivalence of the intrinsic and extrinsic homotopies is that the space intrinsically parameterized in (37) is the same for appropriate choices of parameters as the space that we extrinsically cut out with linear equations before.

Lemma 2.3 For all but a finite number of γ ∈ ℂ of absolute value one, it follows that for any t ∈ [0, 1], the kernel of the linear system

Aw + B(P_{j} + τP_{ji})(Cw + d) = 0.  (39)

on ℂ^{2k} with τ = t/(t + γ(1 − t)) is parameterized by W_{i,j}(t, y) where y ∈ ℂ^{m}.

Proof. This follows immediately for t = 0 and 1 with no restriction on γ of absolute value 1 by taking τ equal to 0 and 1 respectively.

Combining this with the dimension of the kernel of (39) being at least m, we conclude that the dimension of the kernel of (39) is exactly m except for finitely many 0 ≠ γ ∈ ℂ. In particular, for all but a finite number γ of absolute value 1, the dimension of the kernel of (39) for τ = t/(t + γ(1 − t)) with t ∈ (0, 1) is of dimension m. Since ϵ satisfies both Aϵ = 0 and Cϵ + d = 0, it is therefore enough to show that for all (t, y)

(A + B(P_{j} + τP_{ji})C) [ E_{i,j} tF_{i,j} + γ(1 − t)G_{i,j} ] y = 0.  (40)

Since the columns of E_{i,j} are in Null Y_{j} ∩Null (BP_{ji}C), it is annihilated. Since y is arbitrary, we must have

(A + B(P_{j} + τP_{ji})C) [tF_{i,j} + γ(1 − t)G_{i,j}] = 0.  (41)

Since F_{i,j} is in Null Y_{i} and G_{i,j} is in Null Y_{j}, this is the same as

B ((τ − 1)tP_{ji}CF_{i,j} + τγ(1 − t)P_{ji}CG_{i,j}) = 0.  (42)

By Condition 3 of Lemma 2.2, this becomes

((τ − 1)t + τγ(1 − t)) B [ \begin{matrix}0 \\ I_{i−j} \\ 0 \end{matrix} ] ,  (43)

which equals zero since τ = t/(t + γ(1 − t)). □

We rephrase Lemma 2.3.

Lemma 2.4 For all but a finite number of γ ∈ ℂ of absolute value one, it follows that for any t ∈ [0, 1], the system

F (W_{i,j}(t, y)) = 0  (44)

on ℂ^{m} is the intrinsic system associated to the system

[ \begin{matrix}F(w) \\ Aw + BP_{i}z \\ z − (P_{j} + τP_{ji})(Cw + d) \end{matrix} ] = 0  (45)

with τ = t/(t + γ(1 − t)).

We define G_{i} as the set of nonsingular solutions of H_{i,j}(1, ω) on which P_{i}(Cw + d) is nonzero and which correspond to points of A × B; G_{j} as the set of nonsingular solutions of H_{i,j}(0, ω) on which P_{j}(Cw +d) is nonzero and which correspond to points of A × B; and G_{i,j} as the set of limits obtained by tracking G_{i} from t = 1 to t = 0 using the homotopy H_{i,j}(t, ω).

Theorem 2.5 The subset Ŵ_{j} ⊂ G_{i,j} on which P_{j}(Cw+d) is zero contains a set of witness points for the j-dimensional components of A ∩ B. These witness points include deg(Z) distinct points for each irreducible j-dimensional component Z of A ∩ B. Moreover G_{j} ⊂ G_{i,j}.

Proof. The sets G_{i}, G_{j} considered as sets of solutions of the extrinsic systems E_{i}, E_{j} on ℂ^{2k} are the same as the sets occurring in the homotopy of [15]. The extrinsic homotopy from [15] that we discussed in §1 is simply a differentiable path P parameterized by t ∈ [0, 1] on a complex line ℓ in the parameter space of the systems E_{i}(w, z) joining a general point E_{i} to a general point E_{j} of the linear subspace of systems of the from E_{j}(w, z). The only fact about the path P used in [15] is that it depends on a choice of γ ∈ ℂ of absolute value 1, which can be chosen, except for a finite number of complex numbers of absolute value 1, so that P avoids a certain finite subset B of ℓ. In Lemma 2.4 we show that the intrinsic homotopy leads to systems on the same complex line ℓ. What changed is that the path P′ on ℓ is not linearly related to the original path P. But since the path P′ depends on a choice of γ ∈ ℂ of absolute value 1, which can still be chosen, except for a finite number of complex numbers of absolute value 1, so that P′ avoids the finite subset B of ℓ, the same conclusions of [15] still hold. □

2.3 Algorithm Summary

The homotopy algorithm to intersect two positive dimensional varieties in intrinsic coordinates is described below. After the initialization, there are three stages. First is the homotopy to start the cascade, followed by the homotopy to find a witness set for the top dimensional part of A∩B. Thirdly, all lower dimensional parts of A ∩ B are computed in a loop from b − 2 down to h_{0}. The second and third stage are separate because we can avoid a coordinate transformation. Also, in many cases – such as the important application of the intersection with a hypersurface – the loop will never be executed.

Some subroutines used in the algorithm below are just implementations of one formula in the paper, e.g.: Combine implements (15). Next we describe briefly the other subroutines.

The linear algebra operations to deal with solutions in intrinsic coordinates are provided in the subroutines Start Plane, Project, Initialize, Basis, and Transform. Given the equations for L_{A} and L_{B}, Start Plane first computes a basis for the null space of L_{A}^{−1}(0) and L_{B}^{−1}(0) before doubling the coordinates into a corresponding basis in ℂ^{2k}. After orthonormalization of the basis, Project computes the intrinsic coordinates for the product of the given witness sets of A and B. The subroutine Initialize first generates the random matrices A, B, C, and d before computing the ϵ of Lemma 2.1. In addition, Initialize returns the operator Y, which returns for any h the corresponding Y_{h} of (34). Lemma 2.2 is implemented by Basis, while Transform converts the coordinates for the solutions from one basis into another.

The path tracking is done by the procedure Track. On input are the homotopy and start solutions. Except from the set up of the homotopy in intrinsic coordinates, one can implement Track along the lines of general path following methods, see [1], [6, 7], or [9].

The subroutine Filter takes on input the witness sets W for higher dimensional components and the list Z. On return is W, augmented with a witness set for the solution set at the current dimension, and a filtered list Z of nonsolutions. The list Z given to Filter may contain points on higher dimensional solution sets. To remove such points, a homotopy membership test as proposed in [12] can be applied. Recently, an interesting alternative was proposed by Li and Zeng in [8]. The nonsolutions serve as start solutions in the cascade to find witness sets for the lower dimensional solution sets. If Z becomes empty after Filter, the algorithm terminates.

Algorithm 2.6 Intersecting two Positive Dimensional Varieties A and B.

Input: k, a, b (a ≥ b);  dim(A) = a, dim(B) = b, A, B ⊂ ℂ^{k}
f_{A}(u) = 0, f_{B}(v) = 0;  polynomial systems in u, v ∈ ℂ^{k}
L_{A}(u) = 0, L_{B}(v) = 0; dim(L_{A}^{−1}(0)) = k−a, dim(L_{B}^{−1}(0)) = k−b
W_{A}, W_{B}.  solutions in witness sets for A and B
Output: F(x) = 0;  system combined from f_{A}, f_{B} in x ∈ ℂ^{k}
L = [L_{h_{0}}, …, L_{b−1}];  list of linear spaces, dim(L_{i}^{−1}(0)) = i
W = [W_{h_{0}}, …, W_{b−1}].  solutions W_{i} in i-dim witness sets
F := Combine(f_{A}, f_{B});  combine systems f_{A} and f_{B} as in (15)
S := Start Plane(L_{A}, L_{B});  basis for plane defining W_{A} × W_{B}
Z := Project(W_{A} × W_{B}, S);  solutions to start the cascade
[Y, ϵ] := Initialize(k, a, b);  linear space Aw + BP_{h}C(w + d) = 0
[E, F, G] := Basis(Y_{b}, Y_{b−1});  basis for Null Y_{b} and Null Y_{b−1}
W (t, y) := [tS + (1 − t)[ϵ + [E F]];  deform start plane S into [E F]

with t using formula (21)

Z := Track(F(W (t, y)), Z);  homotopy to start the cascade
Z := Track(F, [E, F, G], Z);  find top dimensional component
[W_{b−1}, Z] := Filter(W, Z);  keep witness sets and nonsolutions
h_{0} := max(a + b − k, 0);  minimal dim(A ∩ B)
for j from b − 2 down to h_{0} do  compute witness set at dimension j
[E, F, G] := Basis(Y_{j+1}, Y_{j});  W (t, y) = ϵ + [E tF + γ(1 − t)G]y
Z := Transform(Z, [E, F]);  coordinates into new basis [E F]
Z := Track(F, [E, F, G], Z);  homotopy F(W_{j+1,j}(t, y)) = 0
[W_{j}, Z] := Filter(W, Z);  keep witness sets and nonsolutions

end for.

We will present the details of the algorithm on the first example in the next section below.

3 Numerical Experiments

The algorithms in this paper have been implemented and tested with PHCpack [16], for public release in version 2.3. To compare with our implementation in extrinsic coordinates, we use the same examples as in [15]. All computations were done on a 2.4 GHz Linux machine.

(1) An Example from Calculus. In this example, we intersect a cylinder A with a sphere B. More precisely, A = { (x, y, z) | x^{2}+y^{2}−1 = 0 } and B = {(x, y, z) | (x + 0.5)^{2} + y^{2} + z^{2} − 1 = 0 }. The intersection A ∩ B is a curve of degree four. Since k = 3, a = 2, and b = 2: h_{0} = 1, so there are only two homotopies, each defining four solution paths.

We now illustrate the details of the algorithm on this example.

F := Combine(f_{A}, f_{B}): Since we have the right number of equations for A and B, i.e.: f_{A}(x, y, z) = x^{2} + y^{2} − 1 = 0 and f_{B}(x, y, z) = (x + 0.5)^{2} + y^{2} + z^{2} − 1 = 0, the matrices M and N are not really needed, and F(w) consists of f_{A}(u) = 0 and f_{B}(v) = 0, with w = [u v]^{T}.

S := Start Plane(L_{A}, L_{B}): The random hyperplanes L_{A}(u) = 0 and L_{B}(v) = 0 each define a random line in 3-space. If L_{A}(u) = 0 is represented by u(y) = ϵ_{A} + ρ_{A}y, and if L_{B}(v) = 0 is represented by v(y) = ϵ_{B} + ρ_{B}y, then S is generated by

w = [ \begin{matrix}u \\ w \end{matrix} ] = [ \begin{matrix}ϵ_{A} \\ ϵ_{B} \end{matrix} ] + [ \begin{matrix}ρ_{A} \\ 0 \end{matrix} ] y_{1} + [ \begin{matrix}0 \\ ρ_{B} \end{matrix} ] y_{2}.  (46)

Z := Project(W_{A} × W_{B}, S): The solutions to start the cascade are obtained from appending the coordinates of the two points in W_{B} to every point in the witness set W_{A}. As #W_{A} = 2, there are four start solutions in the cascade. The second argument S of Project is used to project the (x, y, z)-coordinates of the witness sets into the (y_{1}, y_{2})-coordinates in the representation of S in (46).

[Y, ϵ] := Initialize(3, 2, 2): At this stage, random matrices A, B, C, and d are generated as in (16) and in (17), a solution ϵ of Aw + BP_{h}C(w + d) = 0 is computed, and the operator Y is returned.

[E, F, G] := Basis(Y_{2}, Y_{1}): After the execution of Basis, we have [E F] = Null Y_{2} and [E G] = Null Y_{1}, as in Lemma 2.2.

Z := Track(F(W (t, y)), Z): With this homotopy defined by W (t, y) we track four paths to start the cascade.

Z := Track(F, [E, F, G], Z): The four paths tracked here end at the four points which form a witness set for the curve of intersection A ∩ B of degree four.

As b − 2 = 0, h_{0} = 1, the range for j is empty, and the loop in the algorithm is not executed.

(2) An Illustration of the Cascade. In this example we need to execute the cascade to find the point of intersection. We consider the components A = { x = 0, y = 0 } and B = { z = 0, w = 0 } as solution sets of the same system f(x, y, z, w) = [xz, xw, yz, yw]^{T} = 0. We have k = 4, a = 2, and b = 2.

(3) Adding an Extra Leg to a Moving Platform. In this example we cut a hypersurface A in ℂ^{8} with a curve B, i.e.: a = 7 and b = 1. The application concerns a Griffis-Duffy platform [3] (analyzed by Husty and Karger in [5] and subsequently in [14]) where A ∩ B can be interpreted as adding a seventh leg to the platform so it no longer moves. As deg(A) = 2 and deg(B) = 28 (ignoring the mechanically irrelevant components), there are 56 paths to trace, by two homotopies.

In the Table 1 below we list all important dimensions of the three example applications. A summary of the execution times is reported in Table 2.

example dimensions and degrees of A and B  deg(A)
k dim(A) deg(A) dim(B) deg(B) m M × deg(B)
(1) 3  2  2  2  2  2 7  4
(2) 4  2  1  2  1  4 10  1
(3) 8  7  2  1  28  8 17  56

Table 1: Dimension and degrees of the two irreducible sets A and B for the three examples, followed by #variables m = 2k − dim(A) − dim(B), M = 3k −a (which is the #variables in the extrinsic homotopy), and number of paths deg(A) × deg(B) at the start of the cascade.

In these numerical experiments, we save about half of the computational time when working in intrinsic coordinates. Comparing the number of variables of the original extrinsic method, M = 3k − a for the examples tested, with the number for the intrinsic method, m = 2k − dim(A) − dim(B), we have in these experiments 3k − a = 7, 10, 17 variables reduced to 2, 4, 8, or more than half. Since the cost of linear solving is O(n^{3}), this implies about a eight-fold reduction in the cost of linear solving. Linear solving can be

Homotopies  Total CPU Time
0  1  2 intrinsic extrinsic
(1) 0.03 0.01 –  0.04  0.07
(2) 0.01 0.02 0.01  0.04  0.11
(3) 9.90 5.94 –  15.84  34.70

Table 2: Timings in CPU user seconds on 2.4GHz Linux machine. The second column concerns the homotopy to start the cascade, in the third column are the timings for the top dimensional components, followed by the eventual next homotopy in the cascade.

a significant portion of the total cost, as it is used in Newton’s method for tracking the homotopy paths. The experimental results suggest that this was accounting for about half of the total cost in the extrinsic method, but accounts for a much less significant fraction of the computational cost of the intrinsic method. The other 50% or so of the cost remains, which is attributable to function evaluation, data transfer, and other overhead. The cost of function evaluation can vary dramatically from one polynomial system to another, so we cannot definitively expect the same percentage savings for all systems, but we can say that the intrinsic formulation seems to give a substantial reduction in computational time.

References

[1] E.L. Allgower and K. Georg. Numerical Continuation Methods, an Introduction, volume 13 of Springer Ser. in Comput. Math. Springer– Verlag, Berlin Heidelberg New York, 1990. Reprinted in 2003 by SIAM as volume 45 in the Classics in Applied Mathematics series.

[2] D. Eisenbud. Commutative Algebra with a View Toward Algebraic Geometry, volume 150 of Graduate Texts in Mathematics. Springer– Verlag, New York, 1995.

[3] M. Griffis and J. Duffy. Method and apparatus for controlling geometrically simple parallel mechanisms with distinctive connections. US Patent 5,179,525, 1993.

[4] B. Huber, F. Sottile, and B. Sturmfels. Numerical Schubert calculus. J. of Symbolic Computation, 26(6):767–788, 1998.

[5] M.L. Husty and A. Karger. Self-motions of Griffis-Duffy type parallel manipulators Proc. 2000 IEEE Int. Conf. Robotics and Automation, CDROM, San Francisco, CA, April 24–28, 2000.

[6] T.Y. Li. Numerical solution of multivariate polynomial systems by homotopy continuation methods. Acta Numerica 6:399–436, 1997.

[7] T.Y. Li. Numerical solution of polynomial systems by homotopy continuation methods. In Handbook of Numerical Analysis. Volume XI. Special Volume: Foundations of Computational Mathematics, edited by F. Cucker, pages 209–304, 2003.

[8] T.Y. Li and Z. Zeng. A rank revealing method and its applications. Preprint, 2003.

[9] A. Morgan. Solving polynomial systems using continuation for engineering and scientific problems. Prentice-Hall, Englewood Cliffs, N.J., 1987.

[10] D. Mumford. Algebraic Geometry I: Complex Projective Varieties, volume 221 of Grundlehren der mathematischen Wissenschaften 221. Springer–Verlag, Berlin, 1976.

[11] A.J. Sommese and J. Verschelde. Numerical homotopies to compute generic points on positive dimensional algebraic sets. Journal of Complexity 16(3):572–602, 2000.

[12] A.J. Sommese, J. Verschelde, and C.W. Wampler. Numerical irreducible decomposition using projections from points on the components. In Symbolic Computation: Solving Equations in Algebra, Geometry, and Engineering, volume 286 of Contemporary Mathematics, edited by E.L. Green, S. Hoşten, R.C. Laubenbacher, and V. Powers, pages 37–51. AMS 2001.

[13] A.J. Sommese, J. Verschelde, and C.W. Wampler. A method for tracking singular paths with application to the numerical irreducible decomposition. In Algebraic Geometry, a Volume in Memory of Paolo Francia, edited by M.C. Beltrametti, F. Catanese, C. Ciliberto, A. Lanteri, C. Pedrini. W. de Gruyter, pages 329-345, W. de Gruyter, 2002.

[14] A.J. Sommese, J. Verschelde, and C.W. Wampler. Symmetric functions applied to decomposing solution sets of polynomial systems. SIAM J. Numer. Anal. 40(6):2026–2046, 2002.

[15] A.J. Sommese, J. Verschelde, and C.W. Wampler. Homotopies for intersecting solution components of polynomial systems. To appear SIAM J. Numer. Anal.

[16] J. Verschelde. Algorithm 795: PHCpack: A general-purpose solver for polynomial systems by homotopy continuation. ACM Transactions on Mathematical Software 25(2): 251–276, 1999. Software available at

http://www.math.uic.edu/~jan.