A Newton's Iteration Converges Quadratically to Nonisolated Solutions Too
aa r X i v : . [ m a t h . NA ] J a n A Newton’s Iteration Converges Quadratically toNonisolated Solutions Too
Zhonggang Zeng ∗ January 25, 2021
Abstract
The textbook Newton’s iteration is practically inapplicable on nonisolatedsolutions of unregularized nonlinear systems. With a simple modification, aversion of Newton’s iteration regains its local quadratic convergence to non-isolated zeros of smooth mappings assuming the solutions are semiregular asproperly defined regardless of whether the system is square, underdetermined oroverdetermined. Furthermore, the iteration serves as a de facto regularizationmechanism for computing singular zeros from empirical data. Even if the givensystem is perturbed so that the nonisolated solution disappears, the iterationstill locally converges to a stationary point that approximates a solution of theunderlying system with an error bound in the same order of the data accuracy.Geometrically, the iteration approximately converges to the nearest point onthe solution manifold. This extension simplifies nonlinear system modeling byeliminating the zero isolation process and enables a wide range of applicationsin algebraic computation.
Newton’s iteration as we know it loses its quadratic rate of convergence, if it convergesat all, to nonisolated solutions of equations with large errors in numerical computation.A subtle tweak of its formulation restores the fast rate of convergence and the optimalaccuracy as we shall elaborate in this paper.Perhaps there is no need to explain or even mention the importance of Newton’siteration in scientific computing, applied mathematics and numerical analysis. In ∗ Department of Mathematics, Northeastern Illinois University, Chicago, Illinois 60625, USA.email: [email protected] . Research is supported in part by NSF under grant DMS-1620337. x k +1 = x k − J ( x k ) − f ( x k ) for k = 0 , , . . . (1)is the standard method for solving systems of nonlinear equations in the form of f ( x ) = where J ( x ) is the Jacobian of the mapping f at x . It is welldocumented that Newton’s iteration quadratically converges to any isolated solutionunder natural conditions: The mapping is smooth and the initial iterate is near asolution at which the Jacobian is invertible.Solving systems of nonlinear equations is a standard topic in textbooks of numericalanalysis but the discussions have always been limited to isolated solutions of squaresystems. Models with nonisolated solutions frequently arise in scientific computingas we shall elaborate in § semiregular solutions, establish an extensionof Newton’s iteration for such solutions and prove its local quadratic convergence onexact equations along with local linear convergence on perturbed equations with em-pirical data. Furthermore, we provide a geometric interpretation of the convergencetendency, elaborate the modeling and applications involving nonisolated solutions anddemonstrate our software implementation with computing examples.An isolated zero is regular if the Jacobian is invertible, namely the nullity of theJacobian and the dimension of the zero are both zero. Nonisolated zeros of a smoothmapping can form a smooth submanifold of a positive dimension such as curvesand surfaces. We generalize the regularity of isolated zeros to the semiregularity ofnonisolated cases as the dimension being identical to the nullity of the Jacobian at thezero. Semiregular zeros of a positive dimension form branches with locally invariantdimensions (c.f. Lemma 3) and, near a semiregular zero, we prove a crucial propertythat every stationary point is a semiregular zero (c.f. Lemma 4).We extended Newton’s iteration to the form of x k +1 = x k − J rank- r ( x k ) † f ( x k ) for k = 0 , , . . . (2)for a smooth mapping f from an open domain in R m or C m to R n or C n where J rank- r ( x k ) † is the Moore-Penrose inverse of J rank- r ( x k ) that is the rank- r projection2f the Jacobian J ( x k ). We establish its local quadratic convergence (Theorem 1)under minimal natural assumptions: The mapping f is smooth and the initialiterate is near a semiregular zero at which the Jacobian is of rank r .Nonisolated solutions can be highly sensitive to data perturbations. When the sys-tem of equations is perturbed, represented with empirical data or solved using floatingpoint arithmetic, the nonisolated solution can be significantly altered or even disap-pears altogether. We prove that the proposed Newton’s iteration still converges to astationary point that approximates an exact solution of the underlying system withan accuracy in the same order of the data error (c.f. Theorem 2). In other words, theproposed extension of Newton’s iteration also serves as a regularization mechanismfor such an ill-posed zero-finding problem. A condition number can also be derivedfrom Theorem 2 for nonisolated zeros with respect to data perturbations.As a geometric interpretation, we shall illustrate the behavior of the Newton’s it-eration (2) for converging to roughly the point on the solution manifold nearest tothe initial iterate. We also elaborate case studies on mathematical modeling withnonisolated solutions, demonstrate our software implementation in solving for suchsolutions with step-by-step calling sequences and computing results.Extending Newton’s iteration beyond (1) by replacing the inverse with a certain kindof generalized inverses traces back to Gauss for solving least squares solutions ofoverdetermined systems by the Gauss-Newton iteration with an assumption that theJacobian is injective. For systems with rank-deficient Jacobians, Ben-Israel [3] is thefirst to propose using Moore-Penrose inverses to generalize Newton’s iteration as x k +1 = x k − J ( x k ) † f ( x k ) for k = 0 , , . . . (3)“but the conditions for the [convergence] theorem are somewhat restrictive and un-natural” [4]. Chu is the first to prove the local convergence of Newton’s iteration(3) with essentially minimal assumptions for underdetermined systems with surjec-tive Jacobians [6]. Applying the alpha theory, Dedieu and Kim [9] prove Newton’siteration (3) locally quadratically converges to a nonisolated solution under the as-sumption that the Jacobian has a constant deficient rank in a neighborhood of theinitial iterate.Nashed and Chen [21] extend Newton’s iteration further as x k +1 = x k − J ( x k ) f ( x k ) for k = 0 , , . . . (4)where, for any matrix A , the notation A stands for one of the outer inverses of A satisfying the identity A A A ≡ A , and prove its local quadratic convergenceto a stationary point ˆ x at which J ( x ) f (ˆ x ) = using a specific outer inverse J ( x k ) = (cid:0) I + J rank- r ( x ) † (cid:0) J ( x k ) − J ( x ) (cid:1)(cid:1) − J rank- r ( x ) † (5)for a proper r but it is unknown whether the limit ˆ x is a zero of f . Chen, Nashedand Qi [5] along with Levin and Ben-Israel [18] follow up with similar convergenceresults toward stationary points using outer inverses.3n comparison to those pioneer works, our extension (2) is suitable for any rank ofthe Jacobian at the zero. Furthermore, our convergence theorems require minimalassumptions and the iteration quadratically converges to a stationary point that isguaranteed to be a zero if the initial iterate is near a semiregular zero. The iterationalso permits data perturbations and floating point arithmetic while still converges toan approximate zero at linear rate.We loosely refer to our extension (2) as the rank- r Newton’s iteration or simply
New-ton’s iteration following a long-standing practice. The terminology is actually debat-able as it is fair to ask “Is Newton’s method really Newton’s method?” [10], and theterm may even be considered “an enduring myth” [16]. Heavily influenced by Fran¸coisVi`ete, Isaac Newton’s original method is not even an iteration and can be considereda special case of Joseph Raphson’s formulation restricted to univariate polynomialequations. The version (1) of Newton’s iteration can deservedly be credit to ThomasSimpson as well. As a myth or not, however, the term “Newton’s iteration” hasbeen used for all extensions of Newton’s original method and will likely to used as aconvention in the future.
Column vectors are denoted by boldface lower case letters such as b , x , y etc.with being a zero vector whose dimension can be derived from the context. Thevector spaces of n -dimensional real and complex column vectors are denoted by R n and C n respectively. The vector space of m × n complex matrices including realmatrices is denoted by C m × n . Matrices are denoted by upper case letters such as A , B , X , etc. with O and I denoting a zero matrix and an identity matrixrespectively. The range, kernel, rank and Hermitian transpose of a matrix A aredenoted by Range ( A ), Kernel ( A ), rank ( A ) and A H respectively. For any matrix A , its Moore-Penrose inverse [13, § A † is the unique matrix satisfying A A † A = A, A † A A † = A † , ( A A † ) H = A A † , ( A † A ) H = A † A. (6)The j -th largest singular value of A is denoted by σ j ( A ). Let U Σ V H be thesingular value decomposition of A where U = [ u , · · · , u m ] and V = [ v , · · · , v n ]are unitary matrices formed by the left singular vectors and the right singular vectorsrespectively. The rank- r projection A rank- r of A , also known as rank- r truncatedsingular value decomposition (TSVD) and rank- r approximation of A , is defined as A rank- r := σ ( A ) u v H + · · · + σ r ( A ) u r v H r Using singular values and singular vectors, the identity [13, § A † ≡ X σ j ( A ) > σ j ( A ) v j u H j A A † rank- r = A rank- r A † rank- r = [ u , · · · , u r ] [ u , · · · , u r ] H (7) A † rank- r A = A † rank- r A rank- r = [ v , · · · , v r ] [ v , · · · , v r ] H (8)that are orthogonal projections onto the subspaces spanned by { u , . . . , u r } and { v , . . . , v r } respectively. For any matrix A , denote A † rank- r as ( A rank- r ) † .We say f is a smooth mapping if f : Ω ⊂ R n → R m has continuous derivativesof second order, or f : Ω ⊂ C n → C m is holomorphic, where the domain Ω is anopen subset of C n or R n . We may designate a variable name, say x , for f anddenote f as x f ( x ). In that case the Jacobian of f at a particular x ∈ Ωis denoted by f x ( x ) or J ( x ) while J rank- r ( x ) or equivalently f x ( x ) rank- r isits rank- r projection. For a smooth mapping ( x , y ) f ( x , y ) at ( x , y ), thenotation f xy ( x , y ) represents its Jacobian with respect to ( x , y ) while f x ( x , y )and f y ( x , y ) denote the (partial) Jacobians with respect to x and y respectivelyat ( x , y ). Lemma 1
Let J ( x ) be the Jacobian of a smooth mapping f at any x in its opendomain Ω in C n or R n . Assume rank ( J ( x ∗ ) ) = r at x ∗ ∈ Ω . Then there isan open bounded convex subset Ω ∗ ∋ x ∗ of Ω and constants ζ , µ, η, ξ > suchthat, for every x , y ∈ Ω ∗ , the inequality rank ( J ( x ) ) ≥ r holds along with (cid:13)(cid:13) J rank- r ( x ) J rank- r ( x ) † − J rank- r ( y ) J rank- r ( y ) † (cid:13)(cid:13) ≤ ζ k x − y k (9) (cid:13)(cid:13) J rank- r ( x ) − J rank- r ( y ) (cid:13)(cid:13) ≤ µ k x − y k (10) (cid:13)(cid:13) J rank- r ( x ) † − J rank- r ( y ) † (cid:13)(cid:13) ≤ η k x − y k (11) (cid:13)(cid:13) J rank- r ( x ) † (cid:13)(cid:13) ≤ ξ (12) Proof.
From rank ( J ( x ∗ ) ) = r , we have σ r ( J ( x ∗ )) > σ r +1 ( J ( x ∗ )) = 0. Weyl’sTheorem [26, Corollary 4.31, p. 69] ensures the singular values to be continuous withrespect to the matrix entries. By the continuity of J ( x ) with respect to x in Ω,there is an open bounded convex neighborhood Ω ∗ of x ∗ with Ω ∗ ⊂ Ω such that σ r ( J ( x )) > σ r +1 ( J ( x )) for every x ∈ Ω ∗ . We can further assume Ω ∗ to besufficiently small so that (cid:13)(cid:13) J ( x ) − J ( y ) (cid:13)(cid:13) < (cid:0) σ r ( J ( x )) − σ r +1 ( J ( x )) (cid:1) (13)for all x , y ∈ Ω ∗ andmax x ∈ Ω ∗ (cid:13)(cid:13) J rank- r ( x ) † (cid:13)(cid:13) = max x ∈ Ω ∗ σ r ( J ( x )) ≤ (cid:13)(cid:13) J ( x ∗ ) † (cid:13)(cid:13) (14)so (12) is true. The left-hand side of (9) is the distance between the subspacesspanned by the first r left singular vectors of J ( x ) and J ( y ) respectively and,5y Wedin’s Theorem [30] (also see [27, Theorem 4.4]) and (13), there is a constant ζ > (cid:13)(cid:13) J rank- r ( x ) J rank- r ( x ) † − J rank- r ( y ) J rank- r ( y ) † (cid:13)(cid:13) ≤ (cid:13)(cid:13) J rank- r ( x ) † (cid:13)(cid:13) (cid:13)(cid:13) J ( x ) − J ( y ) (cid:13)(cid:13) ≤ ζ k x − y k for all x , y ∈ Ω ∗ . As a result, the inequality (10) follows from (9) and (cid:13)(cid:13) J rank- r ( x ) − J rank- r ( y ) (cid:13)(cid:13) = (cid:13)(cid:13) J rank- r ( x ) J rank- r ( x ) † J ( x ) − J rank- r ( y ) J rank- r ( y ) † J ( y ) (cid:13)(cid:13) (by (6) and (8)) ≤ (cid:13)(cid:13) J rank- r ( x ) J rank- r ( x ) † (cid:13)(cid:13) k J ( x ) − J ( y ) k + (cid:13)(cid:13) J rank- r ( x ) J rank- r ( x ) † − J rank- r ( y ) J rank- r ( y ) † (cid:13)(cid:13) k J ( y ) k . since (cid:13)(cid:13) J rank- r ( x ) J rank- r ( x ) † (cid:13)(cid:13) = 1 and k J ( y ) k is bounded on the compact set Ω ∗ .By [25, Theorem 3.3], there is a constant α > (cid:13)(cid:13) J rank- r ( x ) † − J rank- r ( y ) † (cid:13)(cid:13) ≤ α (cid:13)(cid:13) J rank- r ( x ) † (cid:13)(cid:13) (cid:13)(cid:13) J rank- r ( y ) † (cid:13)(cid:13) (cid:13)(cid:13) J rank- r ( x ) − J rank- r ( y ) (cid:13)(cid:13) leading to (11). ✷ Lemma 2
For F = C or R , let z φ ( z ) be a continuous injective mappingfrom an open set Ω in F n to F m . At any z ∈ Ω , there is an open neighborhood ∆ of φ ( z ) in F m such that, for every b ∈ ∆ , there exists a z b in Ω and anopen neighborhood Σ of z b with k b − φ ( z b ) k = min z ∈ Σ k b − φ ( z ) k . (15) Further assume φ is differentiable in Ω . Then φ z ( z b ) † ( b − φ ( z b )) = . (16) Proof.
Let Σ be an open bounded neighborhood of z such than Σ ⊂ Ω. Since φ is one-to-one and continuous, we have δ = min z ∈ Σ \ Σ k φ ( z ) − φ ( z ) k > . Let ∆ = (cid:8) y ∈ F m (cid:12)(cid:12) k y − φ ( z ) k < δ (cid:9) . Then, for every b ∈ ∆, there exists a z b ∈ Σ such that k φ ( z b ) − b k = min z ∈ Σ k φ ( z ) − b k . For every z ∈ Σ \ Σ ,however, k φ ( z ) − b k ≥ k φ ( z ) − φ ( z ) k − k φ ( z ) − b k > δ > k φ ( z ) − b k ≥ k φ ( z b ) − b k implying z b ∈ Σ and (15). Since a local minimum of k φ ( z ) − b k = ( φ ( z ) − b ) H ( φ ( z ) − b )occurs at the interior point z b ∈ Σ , it is straightforward to verify the equation φ z ( z b ) H ( φ ( z b ) − b ) = and thus (16) from Range ( φ z ( z b ) H ) = Range ( φ z ( z b ) † ). ✷ Semiregular zeros of smooth mappings
A point x ∗ is a zero of a mapping f if x = x ∗ is a solution of the equation f ( x ) = . A common notation f − ( ) stands for the set of all zeros of f . A zero x ∗ of f is isolated if there is an open neighborhood Λ of x ∗ in the domain of f such that f − ( ) ∩ Λ = { x ∗ } or x ∗ is nonisolated otherwise. A nonisolated zero x ∗ of a smooth mapping f may belong to a curve, a surface or a higher dimensionalsubset of f − ( ). We adopt a simple definition of the dimension of a zero as follows.For more in-depth elaboration on the dimension of zero sets, see [2, p. 17]. Definition 1 (Dimension of a Zero)
For F = C or R , let x ∗ be a zero ofa smooth mapping f : Ω ⊂ F m → F n . If there is an open neighborhood ∆ ⊂ Ω of x ∗ in F m such that ∆ ∩ f − ( ) = φ (Λ) where z φ ( z ) is a differentiableinjective mapping defined in a connected open set Λ in F k for a certain k > with φ ( z ∗ ) = x ∗ and rank ( φ z ( z ∗ ) ) = k , then the dimension of x ∗ as a zero of f is defined as dim f ( x ∗ ) := dim ( Range ( φ z ( z ∗ ))) ≡ rank ( φ z ( z ∗ ) ) = k. As a special case, an isolated zero is of dimension zero.
A so-defined k -dimensional zero x ∗ is on a smooth submanifold of dimension k in F m . If the dimension dim f ( x ∗ ) is well-defined, then φ z ( z ∗ ) in Definition 1is of full column rank and there is an open neighborhood Λ ∗ ⊂ Λ of z ∗ suchthat rank ( φ z (ˆ z ) ) ≡ k for all ˆ z ∈ Λ ∗ . Namely the dimension of a zero is locallyinvariant . We shall also say every x ∈ φ (Λ ∗ ) is in the same branch of zeros as x ∗ and, if a zero ˜ x is in the same branch of x ∗ , every zero ˆ x in the same branch of˜ x is in the same branch of x ∗ .An isolated zero x ∗ of f is regular if its dimension 0 is identical to the nullity ofthe Jacobian f x ( x ∗ ). This characteristic of regularity can naturally be generalizedto zeros of higher dimensions for being semiregular as defined below. There aretremendous advantages in computing when zeors are semiregular as we shall elaboratethroughout this paper. Definition 2 (Semiregular Zero)
A zero x ∗ of a smooth mapping x f ( x ) is semiregular if dim f ( x ∗ ) is well-defined and identical to nullity ( f x ( x ∗ ) ) . Namely dim f ( x ∗ ) + rank ( f x ( x ∗ ) ) = the dimension of the domain of f . (17) A zero is ultrasingular if it is not semiregular.
A system of equations f ( x ) = is said to be underdetermined if f is a mappingfrom Ω ⊂ F m to F n with m > n where F = C or R . A solution ofan underdetermined system is always semiregular if the Jacobian is surjective or,equivalently, of full row rank. For instance, let ( u ∗ , v ∗ ) be a zero of a smooth7apping ( u , v ) f ( u , v ) from R k × R m to R m and the partial Jacobian f v ( u ∗ , v ∗ ) is invertible. By the Implicit Mapping Theorem, there is a differentiablemapping u g ( u ) from a neighborhood Λ of u ∗ in R k to R m with g ( u ∗ ) = v ∗ and there is an open neighborhood ∆ of ( u ∗ , v ∗ ) in R k × R m suchthat ∆ ∩ f − ( ) = φ (Λ) where φ ( u ) = ( u , g ( u )) for u ∈ Λ. Furthermore, theJacobian φ u ( u ∗ ) is obviously of full column rank k . As a result, the dimension ofthe zero ( u ∗ , v ∗ ) is k that is identical to the nullity of Jacobian of f at ( u ∗ , v ∗ ),implying ( u ∗ , v ∗ ) is semiregular.Ultrasingular zeros can be isolated multiple zeros [7], isolated ultrasingularity embed-ded in nonisolated zero set (c.f. Example 8) or can form an entire branch of zeros(c.f. Example 9). Like the dimension, semiregularity is also invariant on a branch ofnonisolated zeros as asserted in the following lemma. Lemma 3 (Local Invariance of Semiregularity)
Let x ∗ be a semiregular zeroof a smooth mapping f . Then there is an open neighborhood ∆ ∗ of x ∗ such thatevery ˆ x ∈ ∆ ∗ ∩ f − ( ) is a semiregular zero of f in the same branch of x ∗ .Proof . The assertion is obviously true for isolated regular zeros. Let F be either C or R and the domain of f is an open subset in F n . Assume x ∗ is a semiregular k -dimensional zero of f . There is an open neighborhoods ∆ of x ∗ in F n andthere is an open connected set Λ of a certain z ∗ in F k along with a differentiableinjective mapping φ : Λ → F n such that φ ( z ∗ ) = x ∗ , ∆ ∩ f − ( ) = φ (Λ ), rank ( f x ( x ∗ ) ) = n − k and rank ( φ z ( z ∗ ) ) = k . By the continuity of singularvalues we can assume rank ( φ z ( z ) ) ≡ k for all z ∈ Λ . From f ( φ ( z )) ≡ and f x ( φ ( z )) φ z ( z ) ≡ O for all z ∈ Λ , we have nullity ( f x ( φ ( z )) ) ≥ k for all z ∈ Λ .By the continuity of singular values again, there is an open neighborhood ∆ ∗ ⊂ ∆of x ∗ such that rank ( f x ( x ) ) ≥ n − k for all x ∈ ∆ ∗ . Consequently, everyˆ x ∈ Λ ∩ f − ( ) is a semiregular zero where Λ = φ − (∆ ∗ ). ✷ We shall propose a new version of Newton’s iteration that, under proper conditions,converges to a stationary point ˆ x at which J rank- r (ˆ x ) † f (ˆ x ) = . The followingstationary point property of semiregular zeros ensures that, in a neighborhood of asemiregular zero, all stationary points are semiregular zeros in the same branch. Lemma 4 (Stationary Point Property)
Let x f ( x ) be a smooth mappingwith a semiregular zero x ∗ and r = rank ( f x ( x ∗ ) ) . Then there is an open neigh-borhood Ω ∗ of x ∗ such that, for any ˆ x ∈ Ω ∗ , the equality f x (ˆ x ) † rank- r f (ˆ x ) = holds if and only if ˆ x is a semiregular zero of f in the same branch of x ∗ .Proof. We claim there is a neighborhood Ω of x ∗ such that f (ˆ x ) = forevery ˆ x ∈ Ω with f x (ˆ x ) † rank- r f (ˆ x ) = . Assume this assertion is false. Namelythere is a sequence { x j } ∞ j =1 converging to x ∗ such that f x ( x j ) † rank- r f ( x j ) = but f ( x j ) = for all j = 1 , , . . . . Let z φ ( z ) be the parameterization of the8olution branch containing x ∗ as in Definition 1 with φ ( z ∗ ) = x ∗ . From Lemma 2with any sufficiently large j , there is a ˇ x j ∈ Ω ∗ ∩ f − ( ) = φ (∆) such that k x j − ˇ x j k = min z ∈ ∆ k x j − φ ( z ) k = k x j − φ ( z j ) k (18)at a certain z j with φ ( z j ) = ˇ x j , implying φ z ( z j ) φ z ( z j ) † x j − φ ( z j ) k x j − φ ( z j ) k = φ z ( z j ) k x j − φ ( z j ) k (cid:16) φ z ( z j ) † (cid:0) x j − φ ( z j ) (cid:1)(cid:17) = . (19)We claim ˇ x j → x ∗ as well when j → ∞ . Assume otherwise. Namely there isan ε > N >
0, there is a j > N with k ˇ x j − x ∗ k ≥ ε .However, we have k x j − x ∗ k < ε for all j larger than a certain N , implying k ˇ x j − x j k ≥ k ˇ x j − x ∗ k − k x ∗ − x j k > ε > k x j − x ∗ k that is a contradiction to (18).Since f ( x j ) = , we have x j = ˇ x j and we can assume ( x j − ˇ x j ) / k x j − ˇ x j k converges to a unit vector v for j → ∞ due to compactness. Then0 = lim j →∞ f x ( x j ) † rank- r (cid:0) f (ˇ x j ) − f ( x j ) (cid:1) k x j − ˇ x j k = lim j →∞ f x ( x j ) † rank- r f x ( x j ) (ˇ x j − x j ) k ˇ x j − x j k = f x ( x ∗ ) † f x ( x ∗ ) v (20)by (8) and (9), implying v ∈ Kernel ( f x ( x ∗ )). As a result, span { v } ⊕ Range ( φ z ( z ∗ )) ⊂ Kernel ( f x ( x ∗ ))since f x ( x ∗ ) φ z ( z ∗ ) = O due to f ( φ ( z )) ≡ in a neighborhood of z ∗ . From thelimit of (19) for j → ∞ , we have v ∈ Range ( φ z ( z ∗ )) ⊥ and thus nullity ( f x ( x ∗ ) ) ≥ rank ( φ z ( z ∗ ) ) + 1which is a contradiction to the semiregularity of x ∗ .For the special case where x ∗ is an isolated semiregular zero of f with dimension 0,the above proof applies with ˇ x j = x j so that (20) holds, implying a contradictionto nullity ( f x ( x ∗ ) ) = 0.By Lemma 3, there is a neighborhood Ω of x ∗ such that every x ∈ Ω ∩ f − ( )is a semiregular zero of f in the same branch of x ∗ . Thus the lemma holds forΩ ∗ = Ω ∩ Ω . ✷ Remark 1 (A note on terminology)
Keller [15] argues that the terminology “issomewhat unfortunate” on isolated and nonisolated zeros as he defines the later aszeros at which the Jacobian is not injective. The isolated zeros here are referred toas geometrically isolated in [15].
Nonisolated zeros defined by Keller include multiple0-dimensional zeros and zeros on positive dimensional branches. The term singularzero is broadly accepted to include multiple isolated zero and nonisolated zero (c.f. [2, § ultrasingular zero to distinguish those singularzeros from semiregular ones as defined in Definition 2.9 Convergence theorem on exact equations
Consider the system of equations in the form of f ( x ) = where f : Ω ⊂ F m → F n is a smooth mapping with F = C or F = R . The system can be square ( m = n ),underdetermined ( m > n ) or overdetermined ( m < n ). We propose the iteration x k +1 = x k − J rank- r ( x k ) † f ( x k ) for k = 0 , , . . . (21)for computing a zero x ∗ of f at which the Jacobian J ( x ∗ ) is of rank r particularlywhen x ∗ is on a branch of semiregular nonisolated zeros. We loosely refer to (21)as the rank- r Newton’s iteration or simply
Newton’s iteration since it is identical tothe commonly-known Newton’s iteration when r = m = n . Assume the mapping f is given with exact data. The following theorem establishes the local quadraticconvergence of the iteration (21). We shall consider the equation with empirical datain § Theorem 1 (Convergence Theorem)
Let f be a smooth mapping in an opendomain with a rank r Jacobian J ( x ∗ ) at a semiregular zero x ∗ . For every openneighborhood Ω of x ∗ , there is a neighborhood Ω of x ∗ such that, from everyinitial iterate x ∈ Ω , the rank- r Newton’s iteration (21) converges quadraticallyto a semiregular zero ˆ x ∈ Ω of f in the same branch as x ∗ .Proof. Let Ω ∗ ∋ x ∗ be an open convex subset of the domain for f as specified inLemma 1 so that (9) and (12) hold in Ω ∗ . From Lemma 4, we can further assume J rank- r ( x ) † f ( x ) = implies x is a semiregular zero of f in the same branch of x ∗ for every x ∈ Ω ∗ . Since f is smooth, there are constants µ, γ > (cid:13)(cid:13) f ( y ) − f ( x ) (cid:13)(cid:13) ≤ µ k y − x k (cid:13)(cid:13) f ( y ) − f ( x ) − J ( x ) ( y − x ) (cid:13)(cid:13) ≤ γ k y − x k for all x , y ∈ Ω ∗ . Denote S ε ( x ∗ ) := (cid:8) x (cid:12)(cid:12) k x − x ∗ k < ε (cid:9) for any ε >
0. Forany given open neighborhood Ω of x ∗ , there is a δ with 0 < δ < S δ ( x ∗ ) ⊂ Ω ∗ ∩ Ω with (cid:13)(cid:13) J rank- r ( x ) † (cid:13)(cid:13) (cid:0) γ k x − y k + ζ k f ( y ) k (cid:1) < h < x , y ∈ S δ ( x ∗ ) and, there is a τ with 0 < τ < δ such that (cid:13)(cid:13) J rank- r ( z ) † (cid:13)(cid:13) k f ( z ) k ≤ (1 − h ) δ < δ < z ∈ S τ ( x ∗ ). Then, for every x ∈ S τ ( x ∗ ), we have k x − x ∗ k ≤ k x − x k + k x − x ∗ k ≤ k J rank- r ( x ) † k k f ( x ) k + τ < δ. Namely x ∈ S δ ( x ∗ ). Assume x i ∈ S δ ( x ∗ ) for all i ∈ { , , . . . , k } . Since x j − x j − + J rank- r ( x j − ) † f ( x j − ) = (by (21)) J rank- r ( x j ) † J rank- r ( x j ) J rank- r ( x j ) † = J rank- r ( x j ) † (by (6)) J ( x j − ) J rank- r ( x j − ) † = J rank- r ( x j − ) J rank- r ( x j − ) † (by (7))10or all j ∈ { , , . . . , k } , we have (cid:13)(cid:13) x j +1 − x j (cid:13)(cid:13) = (cid:13)(cid:13) J rank- r ( x j ) † f ( x j ) (cid:13)(cid:13) = (cid:13)(cid:13)(cid:13) J rank- r ( x j ) † (cid:16) f ( x j ) − J ( x j − ) (cid:0) x j − x j − + J rank- r ( x j − ) † f ( x j − ) (cid:1)(cid:17)(cid:13)(cid:13)(cid:13) ≤ (cid:13)(cid:13) J rank- r ( x j ) † (cid:0) f ( x j ) − f ( x j − ) − J ( x j − ) ( x j − x j − ) (cid:1)(cid:13)(cid:13) + (cid:13)(cid:13) J rank- r ( x j ) † (cid:0) J rank- r ( x j ) J rank- r ( x j ) † − J rank- r ( x j − ) J rank- r ( x j − ) † (cid:1) f ( x j − ) (cid:13)(cid:13) ≤ k J rank- r ( x j ) † k (cid:0) γ k x j − x j − k + ζ k f ( x j − ) k (cid:1) k x j − x j − k (24) ≤ h k x j − x j − k leading to (cid:13)(cid:13) x j +1 − x j (cid:13)(cid:13) ≤ h j k x − x k ≤ h j − h δ for all j ∈ { , . . . , k } . Thus k x k +1 − x k ≤ k x k +1 − x k k + · · · + k x − x k ≤ (cid:0) h k + h k − + · · · + 1 (cid:1) k x − x k < − h − h δ = δ and x k +1 ∈ S δ ( x ∗ ) since k x k +1 − x ∗ k ≤ k x k +1 − x k + k x − x ∗ k < (cid:0) + (cid:1) δ = δ, completing the induction so all iterates (cid:8) x j (cid:9) ∞ j =0 of (21) are in S δ ( x ∗ ) from anyinitial iterate x ∈ S τ ( x ∗ ). Furthermore, the iterates (cid:8) x j (cid:9) ∞ j =0 form a Cauchysequence since, for any k, j ≥ k x k + j − x k k ≤ k x k + j − x k + j − k + · · · + k x k +1 − x k k ≤ (cid:0) h j − + · · · + h + 1 (cid:1) k x k +1 − x k k < − h k x k +1 − x k k (25) < − h h k · − h δ = δ h k (26)can be as small as needed when k is sufficient large. Consequently, the sequence (cid:8) x k (cid:9) ∞ k =0 generated by the iteration (21) converges to a certain ˆ x ∈ S δ ( x ∗ ) ⊂ Ω atwhich J rank- r (ˆ x ) † f (ˆ x ) = and thus, by Lemma 4, the limit ˆ x is a semiregularzero in the same branch of x ∗ with the identical dimension.We now have k x k − ˆ x k ≤ − h k x k − x k +1 k ≤ δ h k for all k ≥ j → ∞ in (25) and (26). Substituting k f ( x j − ) k = k f ( x j − ) − f (ˆ x ) k ≤ µ k x j − − ˆ x k ≤ µ (cid:0) k x j − − x j k + k x j − x j +1 k + k x j +1 − x j +2 k + · · · (cid:1) ≤ µ − h k x j − x j − k µ >
0, the inequality (24) yields k x j +1 − x j k < β k x j − x j − k for all j = 1 , , . . . where β = max (cid:26) , ξ (cid:18) γ + ζ µ − h (cid:19)(cid:27) Let Ω = n x ∈ S τ ( x ∗ ) (cid:12)(cid:12)(cid:12) (cid:13)(cid:13) J rank- r ( x ) † f ( x ) (cid:13)(cid:13) < hβ o Then, for every initial iterate x ∈ Ω , we have k x − x k = (cid:13)(cid:13) J rank- r ( x ) † f ( x ) (cid:13)(cid:13) < hβ ≤ h. Consequently, for all j = 1 , , . . . , we have k x j +1 − x j k < β j − k x − x k j = (cid:0) β k x − x k (cid:1) j − k x − x k < h j Thus k x k − ˆ x k = lim j →∞ k x k + j − x k k = lim j →∞ (cid:0) k x k + j − x k + j − k + · · · + k x k +1 − x k k (cid:1) ≤ − h k x k +1 − x k k ≤ − h h k (27)with h <
1. Consequently the convergence to ˆ x is at quadratic rate. ✷ Theorem 1 can serve as the universal convergence theorem of Newton’s iterationincluding the conventional version (1) as a special case. The sufficient conditionsfor Theorem 1 consist of smoothness of f , semiregularity of x ∗ and the initialiterate x being near x ∗ . These assumptions are minimal and indispensable forthe convergence of the rank- r Newton’s iteration (21). The iteration would needto be adapted if the mapping f is not smooth. Without the semiregularity of x ∗ , the limit ˆ x as a stationary point is seldom a zero from our experiment. Thenon-global convergence is an accepted imperfection of Newton’s iteration rather thana drawback. Theoretical criteria on the initial iterate beyond sufficient nearnessto a solution are of little practical meaning as trial-and-error would be easier thanverifying those conditions.The common definition of quadratic convergence of a sequence { x k } to its limit ˆ x would require that there is a constant λ such that k x k − ˆ x k ≤ λ k x k − − ˆ x k for k larger than a certain k and imply k x k − ˆ x k ≤ λ k x k − ˆ x k k − k . The inequality(27) ensures essentially the same error bound for the iterate x k toward ˆ x and canbe used as an alternative definition of quadratic convergence.12 emark 2 (Convergence near an ultrasingular zero) If ˆ x is an ultrasingularzero of f where r = rank ( J (ˆ x ) ), the iteration (21) still locally converges to acertain stationary point ˇ x at which J rank- r (ˇ x ) † f (ˇ x ) = but ˇ x is not necessarily azero of f . Such an ˇ x satisfies the necessary condition for k f ( x ) k to reach a localminimum if rank ( J (ˇ x ) ) = r . From our computational experiment, a stationarypoint to which the iteration (21) converges is rarely a zero of f when the initialiterate is near an ultrasingular zero. Remark 3 (On the projection rank r ) Application of the iteration (21) requiresidentification of the rank r of the Jacobian at a zero. There are various approachesfor determining r such as analytical methods on the application model (c.f. § § θ < (cid:13)(cid:13) J ( x ∗ ) † k − , the numerical rank rank θ ( J ( x ) )within θ is identical to r = rank ( J ( x ∗ ) ) if x is sufficiently close to x ∗ [35].Furthermore, the projection rank needs to be identified or computed only once for asolution branch and the same rank can be used repeatedly in calculating other witnesspoints in the same branch. Practical applications in scientific computing often involve equations that are giventhrough empirical data with limited accuracy. While the underlying exact equationcan have solution sets of positive dimensions, the perturbed equation may not. Suchapplications can be modeled as solving an equation f ( x , y ) = for x ∈ Ω ⊂ C m or R m (28)at a particular parameter value y ∈ Σ ⊂ C n or R n representing the data. Theequation (28) may have a nonisolated solution set only at a particular isolated pa-rameter value y = y ∗ . A typical example is given as Example 3 later in § x exist for an equation f ( x , t ) = only when t = 1exactly. The question becomes: Assuming the equation (28) has a semiregular so-lution x = x ∗ at an underlying parameter y = y ∗ but y ∗ is known only throughempirical data in ˜ y ≈ y ∗ , does the iteration x k +1 = x k − f x ( x k , ˜ y ) † rank- r f ( x k , ˜ y ) , k = 0 , , . . . (29) converge and, if so, does the limit ˜ x approximate a zero x = ˆ x of the mapping x f ( x , y ∗ ) with an accuracy k ˜ x − ˆ x k = O ( k ˜ y − y ∗ k ) in the same order of thedata? The following theorem attempts to answer that question.
Theorem 2 (Convergence Theorem on Perturbed Equations)
Let a smoothmapping ( x , y ) f ( x , y ) be defined in an open domain. Assume x ∗ is a semireg-ular zero of x f ( x , y ∗ ) at a fixed y ∗ with rank ( f x ( x ∗ , y ∗ ) ) = r > and f y ( x ∗ , y ∗ ) k > . Then there exist a neighborhood Ω ∗ × Σ ∗ of ( x ∗ , y ∗ ) , aneighborhood Ω of x ∗ and a constant h with < h < such that, at ev-ery fixed ˜ y ∈ Σ ∗ serving as empirical data for y ∗ and from any initial iterate x ∈ Ω , the iteration (29) linearly converges to a stationary point ˜ x ∈ Ω ∗ at which f x (˜ x , ˜ y ) † rank- r f (˜ x , ˜ y ) = with an error bound k ˜ x − ˆ x k ≤ − h (cid:13)(cid:13) f x ( x ∗ , y ∗ ) † (cid:13)(cid:13) (cid:13)(cid:13) f y ( x ∗ , y ∗ ) (cid:13)(cid:13) k ˜ y − y ∗ k + O (cid:0) k ˜ y − y ∗ k (cid:1) (30) toward a semiregular zero ˆ x of x f ( x , y ∗ ) in the same branch of x ∗ .Proof. Following almost the same proof of Lemma 1 using the smoothness of themapping ( x , y ) f ( x , y ), there is an open bounded convex neighborhood Ω × Σ of ( x ∗ , y ∗ ) along with constants ζ , γ, η > (cid:13)(cid:13) f x (ˆ x , ˆ y ) rank- r f x (ˆ x , ˆ y ) † rank- r − f x (ˇ x , ˆ y ) rank- r f x (ˇ x , ˆ y ) † rank- r (cid:13)(cid:13) ≤ ζ (cid:13)(cid:13) ˆ x − ˇ x (cid:13)(cid:13) (cid:13)(cid:13) f (ˆ x , ˆ y ) − f (ˇ x , ˆ y ) − f x (ˇ x , ˆ y ) (ˆ x − ˇ x ) (cid:13)(cid:13) ≤ γ (cid:13)(cid:13) ˆ x − ˇ x (cid:13)(cid:13) (cid:13)(cid:13) f x (ˆ x , ˆ y ) † rank- r − f x (ˆ x , ˇ y ) † rank- r (cid:13)(cid:13) ≤ η (cid:13)(cid:13) ˆ y − ˇ y (cid:13)(cid:13) (cid:13)(cid:13) f x (ˆ x , ˆ y ) (cid:13)(cid:13) < (cid:13)(cid:13) f x ( x ∗ , y ∗ ) (cid:13)(cid:13) and (cid:13)(cid:13) f y (ˆ x , ˆ y ) (cid:13)(cid:13) < (cid:13)(cid:13) f y ( x ∗ , y ∗ ) (cid:13)(cid:13) for all ˆ x , ˇ x ∈ Ω and ˆ y , ˇ y ∈ Σ withmax (ˆ x , ˆ y ) ∈ Ω × Σ (cid:13)(cid:13) f x (ˆ x , ˆ y ) † rank- r (cid:13)(cid:13) ≤ (cid:13)(cid:13) f x ( x ∗ , y ∗ ) † (cid:13)(cid:13) Let S ε ( x ∗ ) := (cid:8) x ∈ Ω (cid:12)(cid:12) k x − x ∗ k < ε (cid:9) for any ε > h be any fixedconstant with 0 < h <
1. There are constants δ, τ, τ ′ > τ ′ < τ < δ , S τ ( x ∗ ) ⊂ S δ ( x ∗ ) ⊂ Ω and an open neighborhood Σ ⊂ Σ of y ∗ such that (cid:13)(cid:13) f x (ˆ x , ˆ y ) † rank- r (cid:13)(cid:13) (cid:0) γ k ˆ x − ˇ x k + ζ k f (ˆ x , ˆ y ) k (cid:1) < h (cid:13)(cid:13) f x ( z , ˆ y ) † rank- r (cid:13)(cid:13) k f ( z , ˆ y ) k ≤ (1 − h ) δ < δ (cid:13)(cid:13) f x (˜ z , ˆ y ) † rank- r (cid:13)(cid:13) k f (˜ z , ˆ y ) k ≤ (1 − h ) τ < τ for all ˆ x , ˇ x ∈ S δ ( x ∗ ), z ∈ S τ ( x ∗ ), ˜ z ∈ S τ ′ ( x ∗ ) and ˆ y ∈ Σ . Using the sameargument in the proof of Theorem 1, the sequence { x k } generated by the iteration(29) starting from any x ∈ S τ ′ ( x ∗ ) satisfies k x k +1 − x k k < h k x k − x k − k forall k ≥ S τ ( x ∗ ) for every fixed ˜ y ∈ Σ and converges to an ˜ x ∈ S τ ( x ∗ ) ⊂ S τ ( x ∗ ) that depends on the choices of x and ˜ y . Resetting x = ˜ x ∈ S τ ( x ∗ ), Theorem 1 ensures the iteration (29) with y = y ∗ ∈ Σ stays in S δ ( x ∗ ) and converges to a certain semiregular zero ˆ x of themapping x f ( x , y ∗ ). We can assume (cid:13)(cid:13) ˆ y − ˇ y (cid:13)(cid:13) < − h η k f x ( x ∗ , y ∗ ) k for all ˆ y , ˇ y ∈ Σ by shrinking Σ if necessary. Subtracting both sides of˜ x = ˜ x − f x (˜ x , ˜ y ) † rank- r f (˜ x , ˜ y ) x = ˜ x − f x (˜ x , y ∗ ) † rank- r f (˜ x , y ∗ )14ields k ˜ x − x k = (cid:13)(cid:13) f x (˜ x , ˜ y ) † rank- r f (˜ x , ˜ y ) − f x (˜ x , y ∗ ) † rank- r f (˜ x , y ∗ ) (cid:13)(cid:13) ≤ (cid:13)(cid:13) f x (˜ x , ˜ y ) † rank- r − f x (˜ x , y ∗ ) † rank- r (cid:13)(cid:13) k f (˜ x , ˜ y ) − f (ˆ x , y ∗ ) k + (cid:13)(cid:13) f x (˜ x , y ∗ ) † rank- r (cid:13)(cid:13) k f (˜ x , ˜ y ) − f (˜ x , y ∗ ) k . From (cid:13)(cid:13) f x (˜ x , y ∗ ) † rank- r (cid:13)(cid:13) k f (˜ x , ˜ y ) − f (˜ x , y ∗ ) k ≤ (cid:13)(cid:13) f x ( x ∗ , y ∗ ) † (cid:13)(cid:13) (cid:13)(cid:13) f y ( x ∗ , y ∗ ) (cid:13)(cid:13) k ˜ y − y ∗ k and, since f (ˆ x , y ∗ ) = , (cid:13)(cid:13) f x (˜ x , ˜ y ) † rank- r − f x (˜ x , y ∗ ) † rank- r (cid:13)(cid:13) k f (˜ x , ˜ y ) − f (ˆ x , y ∗ ) k < η (cid:13)(cid:13) ˜ y − y ∗ (cid:13)(cid:13) (cid:0) k f x ( x ∗ , y ∗ ) k k ˜ x − ˆ x k + O ( k ˜ y − y ∗ k ) (cid:1) ≤ − h k ˆ x − ˜ x k + O (cid:0) k ˜ y − y ∗ k (cid:1) , we have k ˜ x − x k ≤ (cid:13)(cid:13) f x ( x ∗ , y ∗ ) † (cid:13)(cid:13) (cid:13)(cid:13) f y ( x ∗ , y ∗ ) (cid:13)(cid:13) k ˜ y − y ∗ k + − h k ˆ x − ˜ x k + O (cid:0) k ˜ y − y ∗ k (cid:1) . As a result, k ˆ x − ˜ x k = lim k →∞ k x k − ˜ x k ≤ lim k →∞ (cid:0) k x k − x k − k + · · · + k x − ˜ x k (cid:1) ≤ lim k →∞ ( h k + h k − + · · · + 1) k x − ˜ x k ≤ − h (cid:13)(cid:13) f x ( x ∗ , y ∗ ) † (cid:13)(cid:13) (cid:13)(cid:13) f y ( x ∗ , y ∗ ) (cid:13)(cid:13) k ˜ y − y ∗ k + k ˆ x − ˜ x k + O (cid:0) k ˜ y − y ∗ k (cid:1) , leading to (30). The theorem is proved by setting Ω ∗ × Σ ∗ = S δ ( x ∗ ) × Σ andΩ = S τ ′ ( x ∗ ). ✷ Theorem 2 leads to some intriguing implications. At the exact data parameter value y = y ∗ , solving the equation f ( x , y ∗ ) = for a nonisolated solution in x canbe an ill-posed problem as the structure of the solution can be infinitely sensitive toperturbations on the parameter y . However, Theorem 2 implies that solving thestationary equation f x ( x , y ∗ ) † rank-r f ( x , y ∗ ) = for x ∈ Ωin a neighborhood of a semiregular zero x ∗ is a well-posed problem in the sensethat the solution is Lipschitz continuous with respect to the perturbations on theparameter y from y ∗ with r = rank ( f x ( x ∗ , y ∗ ) ). When the parameter y ∗ isonly known through empirical data ˜ y , the mapping x f ( x , ˜ y ) may not havea positive-dimensional zero of its own near x ∗ and, in some cases, does not have a15ero at all. As a result, solving the equation f ( x , ˜ y ) = in exact sense for x isfutile even if we extend the machine precision. It may be a pleasant surprise thata zero ˜ x of the mapping x f x ( x , ˜ y ) † rank-r f ( x , ˜ y ) exists near x ∗ . Furthermorethat ˜ x approximates a desired zero ˆ x of the underlying mapping x f ( x , y ∗ )with an accuracy in the same order as accuracy of the data. More importantly inpractice, this numerical zero ˜ x is attainable as the rank- r Newton’s iteration (29)locally converges to it.The iteration (29) is more than an algorithm as it is, at the same time, a naturalregularization of a hypersensitive zero-finding problem. Since x ∗ is a nonisolatedzero of the mapping x f ( x , y ∗ ), the iteration does not necessarily converge to x ∗ but to some ˆ x in the same branch of the solution set.As a special case, Theorem 2 reaffirms the local convergence of the Gauss-Newton iter-ation in solving an overdetermined system f ( x ) = b for its least squares solution[11,Theorem 10.2.1, p.222][32] and extends the result with an error bound on the solutionwith respect to the perturbation on the right-hand side b . Corollary 1 (Convergence of Gauss-Newton Iteration)
Let x f ( x ) be asmooth mapping from an open domain Ω ⊂ F n to F m where F = R or C with n < m . Assume x ∗ ∈ Ω is a least squares solution of the overdetermined system f ( x ∗ ) = b with nullity ( f x ( x ∗ ) ) = 0 . Then there are open neighborhoods Ω ∗ and Σ ∗ of x ∗ and f ( x ∗ ) respectively such that the Gauss-Newton iteration x k +1 = x k − f x ( x k ) † ( f ( x k ) − b ) for k = 0 , , · · · linearly converges to x ∗ from any initial iterate x ∈ Ω ∗ provided that b ∈ Σ ∗ and the convergence is quadratic if b = f ( x ∗ ) . Furthermore, the same iteration x k +1 = x k − f x ( x k ) † ( f ( x k ) − ˜ b ) for k = 0 , , · · · on the perturbed system f ( x ) = ˜ b ∈ Σ ∗ from x ∈ Ω ∗ converges to a critical point ˜ x ∈ Ω at which f x (˜ x ) H (cid:0) f (˜ x ) − ˜ b (cid:1) = with an error bound k ˜ x − x ∗ k ≤ c (cid:13)(cid:13) f x ( x ∗ ) † (cid:13)(cid:13) (cid:13)(cid:13) ˜ b − b (cid:13)(cid:13) + O (cid:0) k ˜ b − b k (cid:1) . where c is a moderate constant depending on ˜ x and ˜ b .Proof. The proof is a straightforward verification by applying Theorem 2 on themapping ( x , y ) f ( x ) − y at ( x ∗ , y ∗ ) where y ∗ = f ( x ∗ ) with ˜ y = b andˆ x = x ∗ . ✷ Remark 4 (Condition number of a nonisolated zero)
Another implication ofTheorem 2 lies in the sensitivity measure derived from the error estimate (30), from16hich we can naturally define a condition number of a zero x ∗ of the mapping x f ( x , y ∗ ) with respect to the parameter value y = y ∗ as κ f , x ( x ∗ , y ∗ ) := (cid:26) (cid:13)(cid:13) f x ( x ∗ , y ∗ ) † (cid:13)(cid:13) k f y ( x ∗ , y ∗ ) k if x ∗ is semiregular ∞ otherwise. (31)The condition number of a semiregular zero is finite and, if it is not large, a smallperturbation in the data parameter y results in accurate approximation of the zerowith an error estimate (30). On the other hand, errors of an ultrasingular zero hasno known bound and the condition number is thus infinity. We consider a special case first: Solving a consistent linear system A x = b witha rank r matrix A ∈ C m × n using the rank- r Newton’s iteration (21). Since thesystem is consistent, namely b ∈ Range ( A ), the solution set is an ( n − r )-dimensionalaffine subspace A † b + Kernel ( A ) := (cid:8) A † b + y (cid:12)(cid:12) y ∈ Kernel ( A ) (cid:9) = (cid:8) A † b + N z (cid:12)(cid:12) z ∈ C n − r (cid:9) in which every particular solution is semiregular as defined in Definition 2 wherecolumns of N ∈ C n × ( n − r ) form an orthonormal basis for Kernel ( A ). From anyinitial iterate x ∈ C n , the nearest point in the solution set A † b + Kernel ( A ) is A † b + N z where z = z is the least squares solution to the linear system A † b + N z = x . Namely z = N H ( x − A † b ) = N H x since N H A † = O , and the nearest solution A † b + N z = A † b + N N H x = A † b + ( I − A † A ) x since ( I − A † A ) = N N H are the same orthogonal projection onto Kernel ( A ). Onthe other hand, let the mapping f : C n → C m be defined as f ( x ) = A x − b with the Jacobian J ( x ) ≡ A ≡ J rank- r ( x ). From x ∈ C n , the rank- r Newton’siteration (21) requires only one step x = x − A † ( A x − b ) = A † b + ( I − A † A ) x In other words, the rank- r Newton’s iteration converges to the nearest solution onthe ( n − r )-dimensional solution set from the initial iterate.We now consider a general nonlinear mapping f : Ω ⊂ C m → C n with the Jacobian J ( x ) at any x ∈ Ω. At an iterate x j near a semiregular ( n − r )-dimensional solutionˆ x of the equation f ( x ) = , we have f ( x ) = f ( x j ) + J ( x j ) ( x − x j ) + O ( k x − x j k ) ≈ f ( x j ) + J rank- r ( x j ) ( x − x j ) .
17n the other hand, from f ( x j ) = J (ˆ x ) ( x j − ˆ x ) + O (cid:0) k x j − ˆ x k (cid:1) we have (cid:0) I − J (ˆ x ) J (ˆ x ) † (cid:1) (cid:0) f ( x j ) − f (ˆ x ) (cid:1) = O (cid:0) k x j − ˆ x k (cid:1) and f ( x j ) = J rank- r ( x j ) J rank- r ( x j ) † f ( x j ) + (cid:0) I − J rank- r ( x j ) J rank- r ( x j ) † (cid:1) f ( x j )= J rank- r ( x j ) J rank- r ( x j ) † f ( x j ) + (cid:0) I − J (ˆ x ) J (ˆ x ) † (cid:1) (cid:0) f ( x j ) − f (ˆ x ) (cid:1) + (cid:0) J (ˆ x ) J (ˆ x ) † − J rank- r ( x j ) J rank- r ( x j ) † (cid:1) (cid:0) f ( x j ) − f (ˆ x ) (cid:1) = J rank- r ( x j ) J rank- r ( x j ) † f ( x j ) + O (cid:0) k x j − ˆ x k (cid:1) The basic principle of Newton’s iteration is to designate the next iterate x = x j +1 as a numerical solution of the linear system f ( x j ) + J ( x j ) ( x − x j ) = . (32)Since J (ˆ x ) is rank-deficient at the nonisolated zero ˆ x and x j is near ˆ x , thecoefficient matrix J ( x j ) of the system (32) is expected to be highly ill-conditionedfor being nearly rank-deficient. Consequently, solving the system (32) as it is maynot be reliable. As elaborated in [35, § J rank- r ( x j ) J rank- r ( x j ) † f ( x j ) + J rank- r ( x j ) ( x − x j ) = (33)and pick a proper vector x = x j +1 from the solution in the affine Grassman-nian. The linear system (33) is well-conditioned if the actual sensitivity measure k J (ˆ x ) k k J (ˆ x ) † k ≈ k J rank- r ( x j ) k k J rank- r ( x j ) † k is not large and (33) is an approx-imation to the system (32). The iterate x = x j +1 from (21) is the minimum normsolution of the approximate system (33) and x j +1 − x j ∈ Kernel ( J rank- r ( x j )) ⊥ ≈ Kernel ( J (ˆ x )) ⊥ = Range ( φ z (ˆ z )) ⊥ where z φ ( z ) is the parametrization of the ( n − r )-dimensional zero of f with φ (ˆ z ) = ˆ x , implying the geometric interpretation of the rank- r Newton’s iteration:From an initial iterate, the rank- r Newton’s iteration (21) aims at thenearest point on the solution branch following a normal line toward thesolution set.How accurate the iteration hitting the nearest point depends on several factors includ-ing how far the initial iterate is. Similar geometric interpretations are also observedin [6, 28] in the cases where the Jacobians are surjective.
Example 1 (Normal line direction of convergence)
Consider f : R → R as follows f ( x, y ) = (cid:18) x + x y − x + 2 x + 2 y − x y + y − y − x − y + 3 (cid:19) (34)18igure 1: Each sequence of iteration (21) asymptotically follows a normal line toward thesolution set. Illustration is plotted using actual data in Example 1 from two initial iterates. whose zeros consist of a semiregular 1-dimensional unit circle x + y = 1 anda 0-dimensional isolated point ( − , x , y ) = (1 . , . x, ˆ y ) ≈ (0 . , . x , y ) = (0 . , . . . . . , . . . . ). Both sequences of iterates asymptotically follow cor-responding normal lines of the solution set toward the particular solutions as shownin Figure 1. J rank- r ( x k ) † f ( x k ) Computing the iteration shift J rank- r ( x k ) † f ( x k ) at the iterate x k is the problem ofcalculating the minimum norm solution , or minimum norm least squares solution if b Range ( A rank- r ), of the linear system A rank- r z = b where A ∈ C m × n and r ≤ min { m, n } . (35)The most reliable method is based on the singular value decomposition in the followingAlgorithm 1. Algorithm 1: Computing A † rank- r b by SVD Input: matrix A ∈ C m × n , vector b ∈ C m , integer r ≤ min { m, n } . – By a full or partial SVD, calculate singular values and corresponding leftand right singular vectors σ j , u j , v j , for j = 1 , . . . , r . – calculate y = h u H b σ , . . . , u H r b σ r i ⊤ – calculate z = [ v , . . . , v r ] y A † rank- r b = z If the Jacobian J ( x ∗ ) at the zero x ∗ is of low rank such that r ≪ min { m, n } , apartial SVD such as the USV-plus decomposition [17] can be more efficient. For thecases of r / min { m, n } , we need the following lemma [35, Lemma 13]. Lemma 5
Let A be an m × n matrix with right singular vectors v , . . . , v n and A rank- r be its rank- r projection. Assume σ r ( A ) > σ r +1 ( A ) and N is a matrixwhose columns form an orthonormal basis for span { v r +1 , . . . , v n } . Then, for everydimension- m vector b , the following identity hold: A † rank- r b = ( I − N N H ) (cid:20) µ N H A (cid:21) † (cid:20) (cid:21) . (36)For matrices A ∈ C m × n of large sizes with r ≈ min { m, n } , the SVD can beunnecessarily expensive. There are reliable algorithms that can be substantiallymore efficient. We briefly elaborate some of those algorithms in this section.The simplest case of computing A † rank- r b is when r equals the row dimension n . This case arises when solving an underdetermined system f ( x ) = when theJacobian is surjective at the nonisolated solution. The following algorithm is a wellestablished numerical methods for computing the minimum norm solution of a fullrank underdetermined linear system [13, § Algorithm 2: Solving (35) when r = m < n Input: matrix A ∈ C m × n with m < n , vector b ∈ C m (integer r = m ). – calculate the thin QR decomposition [13, p. 248] A H = Q R – solve the lower-triangular system R H y = b for y– set z = Q y output A † rank- r b = z On the case r < m < n , the following Algorithm 3 is modified from the rank-revealing method in [20].
Algorithm 3: Solving (35) when r < m < n
Input: matrix A ∈ C m × n with m < n , vector b ∈ C m , integer r < m . – calculate the thin QR decomposition [13, p. 248] A H = Q R – set G = R H . – for k = 1 , . . . , m − r do ∗ calculate the vector u k as the terminating iterate of the iteration (see[20] for details) with τ = k A k ∞ y j +1 = y j − (cid:20) τ y H j G k − (cid:21) † (cid:20) τ y H j y j − τG k − y j (cid:21) , j = 0 , , . . . (37)from a random unit vector y ∈ C m As a by product of terminating (37), extract the thin QR decomposi-tion [2 τ u k , G H k − ] H = Q k G k end do – solve for y = y ∗ of the triangular system G H m − r y = Q H m − r (cid:20) m − r b (cid:21) – set z ∗ = Q ( I − U U H ) y ∗ where U = [ u , . . . , u m − r ]output A † rank- r b = z ∗ Notice that G is lower triangular. It is a standard technique in numerical lin-ear algebra to apply Given’s rotation [13, § τ y j , G H k − ] H = Q R where R is in lower triangular form throughout the pro-cess. The algorithm can be easily explained as follows: Let A H = [ Q , ˜ Q ] (cid:20) R O (cid:21) be a full QR decomposition of A H (while Q R is the corresponding thin version).Then Range ( ˜ Q ) is a subspace of Kernel ( A ) whose basis consists of m − r additionalvectors besides columns of ˜ Q . The vectors Q u , . . . , Q u m − r approximately forman orthonormal basis for the vector space spanned by the m − r right singular vectorsof A . The solution y ∗ of the equation G H m − r y = Q H m − r b is the least squaressolution of the linear system (cid:20) τ U H R H (cid:21) y = (cid:20) (cid:21) . Then it is a straightforward verification using Lemma 5 that z ∗ = Q ( I − U U H ) y ∗ is the least squares solution of A rank- r z = b that is orthogonal to columns of[ Q U, ˜ Q ].For the cases of solving (35) where m ≥ n , the kernel or the partial kernel of A generally can not be computed from a QR decomposition of A H like the casesfor m < n . As a result, a basis { u , . . . , u n − r } for Kernel ( A rank- r ) needs to becomputed as shown in Algrithm 4 below. Algorithm 4: Solving (35) when r ≤ n ≤ m Input: matrix A ∈ C m × n with m ≥ n , vector b ∈ C m , integer r ≤ n . – calculate the thin QR decomposition [13, p. 248] A = Q R – for k = 1 , . . . , n − r do ∗ calculate the vector u k as the terminating iterate of the iteration (see[20] for details) with τ = k A k ∞ y j +1 = y j − (cid:20) τ y H j R k − (cid:21) † (cid:20) τ y H j y j − τR k − y j (cid:21) , j = 0 , , . . . (38)from a random unit vector y ∈ C m As a by product of terminating (37), extract the thin QR decomposi-tion [2 τ u k , R k − ] H = Q k R k end do – solve for y = y ∗ of the triangular system R n − r y = Q H n − r (cid:20) n − r b (cid:21) – set z ∗ = ( I − U U H ) y ∗ where U = [ u , . . . , u n − r ]output A † rank- r b = z ∗ Figure 2:
Solution sets of the system in (39)
Models with nonisolated solu-tions arise in many applications.Unaware of the capability ofNewton’s iteration in computingsuch solutions, scientific com-puting practitioners go to greatlengths to make solutions iso-lated by various techniques suchas adding auxiliary equationsand variables. We shall elab-orate case studies in which non-isolated solutions can naturallybe modeled into well-posed com-putational problems. The di-mensions of the solution setsand their semiregularity mayalso be obtained analytically inthe modeling process so that theprojection rank r of Newton’siteration (21) becomes readilyavailable. The rank- r Newton’siteration for any integer r is implemented in our Numerical Algebraic Computationtoolbox NAClab for Matlab as a functionality Newton . We shall demonstrate theimplementation and its effectiveness in computing nonisolated solutions with no needfor auxiliary equations.
One of the main subjects of numerical algebraic geometry and its application inalgebraic kinematics is computing solutions of positive dimensions to polynomial sys- http://homepages.neiu.edu/ ∼ zzeng/naclab.html r Newton’s iteration (21) can beapplied to calculating semiregular witness points directly on the polynomial systemswithout needing extra equations and variables.
Example 2 (A polynomial system)
An illustrative example for nonisolated so-lutions of different dimensions is given in [2, p. 143] as follows: f ( x, y, z ) = ( y − x ) ( x + y + z −
1) ( x − z − x ) ( x + y + z −
1) ( y − y − x ) ( z − x ) ( x + y + z −
1) ( z − . (39)Among the solution sets, a point (1 , , { y = x , z = x } along withthree lines, and a surface { x + y + z = 1 } are of dimensions 0, 1 and 2respectively, as shown in Figure 2. Zeros are semiregular except at intersectionsof the solution sets. The iteration (21) converges at quadratic rate toward thosesolutions from proper initial iterates by setting the projection rank r = 3, 2 and1 respectively. Example 3 (Perturbed cyclic-4 system)
Cyclic- n roots are among the bench-mark problems in solving polynomial systems. Those systems arise in applicationssuch as biunimodular vectors, a notion traces back to Gauss [12]. In general, everycyclic- n system possesses positive dimensional solution sets if n is a multiple of aperfect square [1].We simulate practical computation with empirical data through the cyclic-4 systemin x = ( x , x , x , x ) ∈ C with a parameter t ∈ C : f ( x , t ) := x + x + x + x t x x + x x + x x + x x x x x + x x x + x x x + x x x x x x x − (40)With t ∗ = 1, the equation f ( x ,
1) = is the cyclic-4 system whose solutionconsists of two 1-dimensional sets { x = − x , x = − x , x x = ± , t = 1 } . (41)All zeros in (41) are semiregular except eight ultrasingular zeros in the form of( ± , ± , ± , ±
1) and ( ± i, ± i, ± i, ± i ) with proper choices of signs. When the pa-rameter t is perturbed from t ∗ = 1 to any other nearby value ˜ t , the 1-dimensionalsolution sets dissipate into 16 isolated solutions. Consequently, the parameter value t ∗ = 1 is a bifurcation point at which the solution changes structure.23or instance, suppose the parameter t ∗ is known approximately, say ˜ t = .Even though the solutions of f ( x , ˜ t ) = are all isolated, Theorem 2 ensures theiteration x k +1 = x k − f x (cid:0) x k , ˜ t (cid:1) † rank-3 f ( x k , ˜ t ) , k = 0 , , . . . converges to a numerical solution ˜ x as a zero of the mapping x f x ( x , ˜ t ) † rank-3 f ( x , ˜ t )that approximates a point in the 1-dimensional solution set of the underlying system f ( x ,
1) = with an error in the order of | ˜ t − t ∗ | = 10 − if x is sufficiently close toa semiregular point in the solution set (41), as demonstrated in the following callingsequence applying the NAClab module
Newton : >> P = {’x1+x2+x3+x4’,’0.9999*x1*x2+x2*x3+x3*x4+x4*x1’,...; % enter the perturbed cyclic-4 system as’x1*x2*x3+x2*x3*x4+x3*x4*x1+x4*x1*x2’, ’x1*x2*x3*x4-1’}; % a cell array of character strings>> v = {’x1’,’x2’,’x3’,’x4’}; % enter cell array of the variable names>> J = PolynomialJacobian(P,v); % Jacobian of P w.r.t. the variable names in v>> f = @(x,P,J,v) PolynomialEvaluate(P,v,x); % the function handle for evaluate the system P at x>> fjac = @(x,x0,P,J,v) PolynomialEvaluate(J,v,x0)*x; % function for evaluating J at x>> domain = ones(4,1); parameter = {P,J,v}; % domain (space of 4x1 vectors) and parameters>> z0 = [0.8;1.2;-0.8;-1.2]; % initial z0>> [z,res,fcond] = Newton({f,domain,parameter},{fjac,3},z0,1); % call rank-3 Newton iteration on f% projection from z0 using display type 1Step 0: residual = 7.8e-02Step 1: residual = 2.4e-03 shift = 2.4e-02Step 2: residual = 1.0e-04 shift = 6.8e-04Step 3: residual = 1.0e-04 shift = 5.8e-07Step 4: residual = 1.0e-04 shift = 4.3e-13Step 5: residual = 1.0e-04 shift = 3.6e-16Step 6: residual = 1.0e-04 shift = 1.5e-16 The iteration does not converges to a solution to f ( x , . as the residual k f ( x j , . k can only be reduced to 10 − but the shifts k x j +1 − x j k = k f x ( x j , . † rank-3 f ( x j , . k −→ . × − approaches to the unit roundoff, indicating the iteration converges to a stationarypoint ˜ x at which f x (˜ x , . † rank-3 f (˜ x , . . The module Newton termi-nates at ˜ x that approximates the nearest solution ˆ x of the underlying equation f ( x ,
1) = where˜ x = ( )ˆ x = ( )with a forward error 2 . × − much smaller than the data error 10 − = | ˜ t − t ∗ | as asserted in Theorem 2. As a confirmation of the geometric interpretation of theiteration elaborated in §
6, the point ˆ x is also only about 6 . × − away from thenearest point in the solution set to the initial iterate x = ( ). Example 4 (A bifurcation model)
The parameterized cyclic-4 system (40) leads24o a bifurcation model in the form of the equation f ( x , t ) = for ( x , t ) ∈ R × R depending on the parameter t . The solution of the equation f ( x , t ) = for( x , t ) ∈ R × R consists of 16 branches x = ψ j ( t ) for j = 1 , . . . ,
16 and t ∈ R .On the plane t = 1, eight pairs of those 16 branches intersect at 8 bifurcation pointsthat are embeded in two additional branches of solution curves (41). All solutionsare semiregular except at the intersections. The parameter value t ∗ = 1 is thebifurcation point at which the solution changes the structure.After finding ˜ x in Example 2 as an approximate zero of x f ( x , t ∗ ) from theempirical data ˜ t = 0 . x , ˜ t ) is very close to a branch in (41).The point (˜ x , ˜ t ) is likely to be closer to a solution branch in (41) than to other 17branches. If so, the iteration( x j +1 , t j +1 ) = ( x j , t j ) − f x t ( x j , t j ) † rank-4 f ( x j , t j ) , j = 0 , , . . . starting from the initial iterate ( x , t ) = (˜ x , ˜ t ) converges roughly to the nearestpoint in the branch in (41) on the plane t = 1 in R × R by the geometric inter-pretation elaborated in §
6. This is significant because the sequence { t j } convergesto the important bifurcation point t ∗ = 1. This observation can easily be confirmedby NAClab
Newton and obtains a solution ( x , x , x , x , t ) as( )that is accurate with at least 15 digits, particularly the bifurcation point at t = 1 . (cid:13)(cid:13) f ( x j , t j ) (cid:13)(cid:13) reduces from 10 − to 4 . × − , indicating convergenceto a solution of f ( x , t ) = . Also notice the x component is identical to nearestpoint ˆ x for at least 15 digits, confirming the geometric interpretation in § Remark 5 (Comparison with existing methods)
Regularization has been re-quired for computing nonisolated solutions. An established strategy for regularizingsolutions of positive dimensions is called linear slicing [24, § f ( x ) = canthen be squared by randomization in the form of Λ f ( x ) = where Λ is a randommatrix of proper size [24, § • The singular solutions may disappear when the system is given with empiricaldata as we simulated in Example 3. The augmented slicing equations may notintersect a solution. • Does the randomized square system with perturbed data has zeros approximatethe solutions of the underlying system with exact data?25ssuming the linear slicing isolates regular solutions, Theorem 2 actually implies thatdirect applying the Gauss-Newton iteration on the regularized system without ran-domization results in local linear convergence to a stationary point approximating anunderlying solution even if the solution set is annihilated by the data perturbation.The randomization likely achieves a similar result. From our experiment, linearslicing combined with randomization produces end results of near identical accuracywith a moderate increase of condition numbers. The advantage of the rank- r New-ton’s iteration is to eliminate the need of the solution isolation and randomizationprocesses and to enable direct convergence to the solution manifolds.
For a polynomial pair p and q of degrees m and n respectively with a greatestcommon divisor (GCD) u of degree k along with cofactors v and w , the GCDproblem can be modeled as a zero-finding problem f ( u, v, w ) = (0 , ∈ P m × P n (42)for ( u, v, w ) ∈ P k × P m − k × P n − k with the holomorphic mapping f : P k × P m − k × P n − k −→ P m × P n ( u, v, w ) ( u v − p, u w − q ) (43)where P m is the vector space of polynomials with degrees up to m , etc. Throughstandard bases of monomials, the mapping f in (43) can be represented by a mappingfrom C m + n − k +3 to C m + n +2 . Let ( u ∗ , v ∗ , w ∗ ) denote a particular solution of (42).The general solution of (42) is a 1-dimensional set (cid:8) (cid:0) t u ∗ , t v ∗ , t w ∗ (cid:1) (cid:12)(cid:12) t ∈ C \ { } (cid:9) (44)on which the Jacobian is rank-deficient. Adding one auxiliary equation, however,the Jacobian of the modified mapping becomes injective [33], implying the Jacobianof f in (43) is of rank deficient by one columnwise. As a result, the set (44)consists of semiregular zeros of the mapping f in (43) and the iteration (21) islocally quadratically convergent by setting the projection rank r = ( k + 1) + ( m − k + 1) + ( n − k + 1) − m + n − k + 2 . Extra equations are not needed. Although there are infinitely many solutions forminga 1-dimensional set, there is practically no difference in finding anyone over the other.
Example 5 (Numerical greatest common divisor)
For a simple demo of theiteration (21) and its
NAClab implementation
Newton , let p = − . − . x − x − . x − . x − x , and q = − . + x + x x (45)26hat are considered perturbed data of an underlying polynomial pair with a GCD1 + x + x . The process of identifying the GCD degree and computing the initialapproximations of the GCD along with the cofactors can be found in [33]. With NAClab installed, the following sequence of Matlab statements carries out the rank-8Newton’s iteration( u j +1 , v j +1 , w j +1 ) = ( u j , v j , w j ) − J rank-8 ( u j , v j , w j ) † f ( u j , v j , w j )for j = 0 , , . . . . >> p = ’-1.3333-2.3333*x-4*x^2-3.6667*x^3-2.6667*x^4-x^5’; % enter p as a character string>> q = ’-1.9999+x+x^2+3*x^3’; % enter q similarly>> f = ... % enter the mapping as function handle for (u,v,w) -> (u*v-p, u*w-q)@(u,v,w,p,q) {pminus(ptimes(u,v),p),pminus(ptimes(u,w),q)};>> J = ... % enter the Jacobian J at (u0,v0,w0) as the mapping (u,v,w) -> (u0*v+u*v0, u0*w+u*w0)@(u,v,w,u0,v0,w0,p,q){pplus(ptimes(u0,v),ptimes(u,v0)),pplus(ptimes(u0,w),ptimes(u,w0))};>> domain = {’1+x+x^2’,’1+x+x^2+x^3’, ’1+x’}; % representation of the domain for the mapping f>> parameter = {p,q}; % parameters for the mapping f>> u0 = ’x^2+1.4*x+1.6’; v0 = ’-1.5-x-1.6*x^2-x^3’; w0 = ’-2+2.8*x’; % initial (u0,v0,w0)>> [z,res,fcond] = Newton({f,domain,parameter},{J,8},{u0,v0,w0},1); % Newton on f with J in% rank-8 projection from the initial iterate (u0,v0,w0) with display setting 2Step 0: residual = 1.5e+00 1.6+14.*x+x^2Step 1: residual = 1.1e-01 shift = 4.9e-01 1.114771189108 + 1.114523702576*x + 1.088690420621*x^2Step 2: residual = 1.2e-03 shift = 5.9e-02 1.089839864820 + 1.090023690710*x + 1.089788605779*x^2Step 3: residual = 8.4e-06 shift = 1.0e-03 1.089756319215 + 1.089767203330*x + 1.089783432095*x^2Step 4: residual = 8.3e-06 shift = 1.4e-07 1.089756333892 + 1.089767171466*x + 1.089783428226*x^2Step 5: residual = 8.3e-06 shift = 5.1e-13 1.089756333892 + 1.089767171466*x + 1.089783428226*x^2Step 6: residual = 8.3e-06 shift = 7.0e-15 1.089756333892 + 1.089767171469*x + 1.089783428226*x^2 The computed GCD . + . x + . x is of the scaling independent distance1 . × − that is in the same order of the data error 2 . × − .Similar to Example 4, the system with inexact data for ( p, q ) does not have asolution as the residual k f ( u j , v j , w j ) k can not be reduced below 8 . × − butthe shift k ( u j +1 , v j +1 , w j +1 ) − ( u j , v j , w j ) k reduces to near unit roundoff, implying J rank-8 ( u j , v j , w j ) † f ( u j , v j , w j ) approaches zero. Remark 6 (On existing methods)
The 1-dimensional solution (44) can be iso-lated by linear slicing (c.f. Remark 5) such as requiring the GCD to be monic or tosatisfy a linear equation [33]. The rank- r Newton’s iteration achieves the results ofthe same or better accuracy in practical computation without needing an isolationprocess.
Accurate computation of defective eigenvalues requires regularization due to hyper-sensitivities to data perturbations [34]. Let ˆ λ ∈ C be an eigenvalue of A ∈ C n × n with what we call a multiplicity support m × k . Namely ˆ λ is of the geometric27ultiplicity m with the smallest Jordan block of size k × k . The problem ofcomputing ˆ λ can be naturally modeled as solving the equation g ( λ, X, A ) = O ∈ C n × k (46)for ( λ, X ) ∈ C × C n × k where g : C × C n × k × C n × n −→ C n × k ( λ, X, G ) ( G − λ I ) X − X S (47)with a constant nilpotent upper triangular matrix parameter S ∈ C k × k of rank k −
1. It can be verified with a standard linear algebra that the solution set of (46)is of dimension ℓ = m k in the form of ( λ ∗ , α X + · · · + α ℓ X ℓ ) with a constanteigenvalue component λ ∗ . By [34, Lemma 2(ii)], the Jacobian g λX (ˆ λ, ˆ X, A ) is ofnullity m k , implying the solution is semiregular. Similar to the case study in § n k − m k + 1) Newton’s iteration (21). Example 6 (Defective eigenvalue computation)
We demonstrate the effective-ness of the rank- r Newton’s iteration in modeling and computation of defective eigen-values with a simple example. Let the matrix A and a perturbed matrix E be asfollows. A = − − − − − − − − − − − − and E = 10 − . − . − . − . . . − . . − . − . . . . − . − . − . . . − . . . . − . . − . − . − . − . − . − . − . − . . − . − . − . Matrix A ∈ C × possesses an exact eigenvalue λ ∗ = 3 with the multiplicitysupport 2 × λ = 2 . NAClab module
LinearSolve to solve the homogeneous linearsystem g ( λ , X, A ) = O with the rank-8 projection for a solution X = X andobtain the initial iterate ( λ , X ). With A and E above entered in Matlab, weapply the rank- r Newton’s iteration (21) with r = ( n − m ) k + 1 = 9 by executingthe NAClab module Newton in the following calling sequence. >> g = @(lambda,X,G,S) G*X-lambda*X-X*S; % enter the mapping g as a function handle>> J = @(lambda,X,lambda0,X0,G,S) G*X-lambda*X0-lambda0*X-X*S; % enter the Jacobian>> domain = {1,ones(n,k)}; % representation of the domain for the mapping g>> parameter = {A,S}; % parameters of the mapping g>> [z,res,fcond] = ... % Call Newton on the mapping g with J in rank-9 projection fromNewton({g,domain,parameter},{J,9},{lambda0,X0},2); % (lambda0,X0) using display type 2Step 0: residual = 2.5e-02 2.900000000000000Step 1: residual = 4.1e-03 shift = 1.0e-01 3.000493218848026Step 2: residual = 3.0e-06 shift = 6.5e-03 2.999999313752690Step 3: residual = 2.6e-12 shift = 4.1e-06 3.000000000001909 tep 4: residual = 4.4e-16 shift = 3.2e-12 3.000000000000000Step 5: residual = 2.2e-16 shift = 5.1e-16 3.000000000000000 The λ components of the iteration accurately converges to the defective eigenvalue λ ∗ = 3 at roughly the quadratic rate with an accuracy at the order of the unitroundoff.We simulate practical computation with imperfect empirical data by using the per-turbed matrix ˜ A = A + E in the iteration (cid:0) ˜ λ j +1 , ˜ X j +1 (cid:1) = (cid:0) ˜ λ j , ˜ X j (cid:1) − g λX (cid:0) ˜ λ j , ˜ X j , A + E (cid:1) † rank-9 g (cid:0) ˜ λ j , ˜ X j , A + E (cid:1) . The iteration reaches the numerical solution (˜ λ, ˜ X ) where ˜ λ = 3 . . × − about the same as the data perturbation k E k ≈ . × − . Example 7 (Nearest matrix with the defective eigenvalue)
Furthermore, theiteration (21) can again be applied to refine the eigenvalue computation by solvingthe equation g ( λ, X, G ) = O for ( λ, X, G ) ∈ C × C n × k × C n × n (48)in the specific iteration (cid:0) ˆ λ j +1 , ˆ X j +1 , ˆ A j +1 (cid:1) = (cid:0) ˆ λ j , ˆ X j , ˆ A j (cid:1) − g λXG (cid:0) ˆ λ j , ˆ X j , ˆ A j (cid:1) † rank-12 g (cid:0) ˆ λ j , ˆ X j , ˆ A j (cid:1) (49)for j = 0 , , . . . , starting from (˜ λ , ˜ X , ˆ A (cid:1) = (˜ λ, ˜ X, A + E ). The projection rank r = n k = 12 because, by [34, Lemma 2(ii)], the (full) Jacobian of g is surjective,implying the solution is semiregular with dimension (1 + n k + n ) − n k = n + 1.In this example, the refinement (49) needs only 1 step to reduce the residual of theequation (48) from 2 . × − to 1 . × − . The iteration terminates at (cid:0) ˆ λ, ˆ X, ˆ A (cid:1) with ˆ λ ≈ . A is the matrix withan exact defective eigenvalue ˆ λ of multiplicity support m × k , implying the backwarderror is k ˜ A − ˆ A k ≈ . × − . From the geometric interpretation elaborated in §
6, the iteration (49) converges to (ˆ λ, ˆ X, ˆ A ) that is roughly the nearest point in thesolution set of (48) from the initial iterate (˜ λ, ˜ X, A + E ). Remark 7 (On the existing method)
The method proposed in [34] appears to bethe only published method for accurate computation of defective eigenvalues. Themethod employs m k auxiliary equations in the eigenvalue model to isolate a solutionpoint in the solution set of (46). Although the computational results are practicallyequivalent for Example 6, those arbitrary extra equations are unnatural and are nowunnecessary. Linear slicing for finding the nearest matrix with defective eigenvaluesin Example 7 requires quadruple the number of equations (37 extra equations on top29f the 12 existing ones) unnecessarily. The rank- r Newton’s iteration substantiallyreduces the computing cost in problems such as Example 7.
As formulated in §
3, ultrasingularity occurs if the nullity of the Jacobian is higher at azero than the dimension of the zero. Difficulties arise in computing ultrasingular zerosincluding slow convergence rate of iterative methods (c.f. [8]) and, more importantly,barriers of low attainable accuracy [23, 31]. As pointed out in [15], “ direct attempts atsuch computation may easily fail or give inconclusive results”. Quadratic convergencerate of Newton’s iteration at semiregular zeros is not expected at ultrasingular zeros.At an isolated ultrasingular zero x ∗ of a smooth mapping f : Ω ⊂ F m → F n where F = R or C , the Jacobian J ( x ∗ ) is of rank r < m . At k -th step of conventionalNewton’s iteration (1) or the Gauss-Newton iteration when n > m , the Jacobian J ( x k ) is usually highly ill-conditioned so that the computation of the iterate x k +1 isgenerally inaccurate, substantially limiting the attainable accuracy of the computedzero even if the iteration converges.A depth-deflation strategy [7] deflates the isolated ultrasingularity and transforms thezero x ∗ into a component of a possibly regular zero ( x ∗ , y ∗ ) of an expanded mapping g : ( x , y ) (cid:0) f ( x ) , J ( x ) y , R y − e (cid:1) (50)where R is a random ( m − r ) × m matrix and e = . If ( x ∗ , y ∗ ) is stillan ultrasingular zero of ( x , y ) g ( x , y ), the deflation process can be continuedrecursively. It is proved in [7] that the number of deflation steps is bounded by the so-called depth of an ultrasingular isolated zero. In practice, however, one deflation stepis likely to be enough except for the cases where the breadth nullity ( J ( x ∗ ) ) = 1.At the terminating step of depth deflation, the ultrasingular zero x ∗ of f is acomponent of the semiregular zero of the final expanded mapping. As a result,the Gauss-Newton iteration locally converges at quadratic rate. More importantly,the zero x ∗ can be computed with an accuracy proportional to the data precisionor unit round-off, circumventing the barrier of the perceived attainable accuracy atultrasingular zeros. An earlier deflation strategy in [19] is proven to terminate withthe number of steps bounded by the multiplicity. A so-called strong deflation methodin symbolic-numerical computation is proposed in [14] and also proved to terminatein finitely many steps.The deflation strategy applies to systems at ultrasingular nonisolated solutions aswell. We illustrate the deflation process in the following examples. Example 8 (Isolated ultrasingularity embedded in a nonisolated zero set)
In the cyclic-4 system in Example 3, the mapping x f ( x ,
1) as in (40) possesses8 ultrasingular zeros embedded in the two solution curves (41). Those 8 points are30onisolated zeros with isolated ultrasingularity since each point is a unique ultrasin-gular zero in a small open neighborhood. For instance, the point x ∗ = (1 , − , − , rank ( f x ( x ∗ ,
1) ) = 2. Interestingly, for almostall R ∈ R × , there is a unique y ∗ ∈ R such that ( x ∗ , y ∗ ) is a semiregu-lar isolated zero of the expanded system ( f ( x , , f x ( x , y , R y − (1 , x ∗ , y ∗ ). Evenif the system is given with imperfect empirical data, the Gauss-Newton iteration onthe expanded system still converges linearly to a stationary point that approximates( x ∗ , y ∗ ) with an accuracy at the same order of the data. The depth deflation methodbeing applicable to nonisolated solutions of isolated ultrasingularities has apparentlynot been observed before. The theoretical termination and the bound on the numberof deflation steps are still unknown. Example 9 (An ultrasingular branch of zeros)
Nonisolated zeros can be ultra-singular on an entire branch. Let x = ( x , x , x , x ) and the mapping f ( x ) = x + x + x x − x + x + x x − x + x + x x − with a 1-dimensional solution curve S = { x = (0 , , s, /s ) | s = 0 } on whichthe Jacobian J ( x ) is of rank 1. All the solutions in S are ultrasingular due to nullity ( J ( x ) ) = 3 = 1 = dim f ( x ) on S . As a result, the equation f ( x ) = isunderdetermined since, for every x ∈ S and almost all matrix R ∈ R × , there is aunique y ∈ R such that ( x , y ) is a 1-dimensional zero of the expanded mapping( x , y ) g ( x , y ) in (50) where e can be any nonzero vector, say (1 , , x k +1 , y k +1 ) = ( x k , y k ) − g xy ( x k , y k ) † rank-7 g ( x k , y k ) , k = 0 , , . . . can be carried out using the following NAClab calling sequence. >> P = {’x1^3+x2^2+x3^2*x4^2-1’; ’x1^2+x2^3+x3^2*x4^2-1’; ’x1^2+x2^2+x3^3*x4^3-1’}; % enter system>> x = {’x1’;’x2’;’x3’;’x4’}; J = pjac(P,x); % variable name array x, Jacobian of P w.r.t. x>> y = [’y1’;’y2’;’y3’;’y4’]; R = Srand(3,4); % expanded variable array y, random 3x4 matrix>> F = [P; ptimes(J,y); pminus(ptimes(R,y),[1;0;0])]; K = pjac(F,[x;y]); % new system and Jacobian>> g = @(u,v,F,K,x,y) peval(F,[x;y],[u;v]); % function handle for expanded mapping g(x,y)>> gjac = @(u,v,u0,v0,F,K,x,y) peval(K,[x;y],[u0;v0])*[u;v]; % function handle for expanded Jacobian>> u0 = [0.001; 0.003; 0.499; 2.002]; % an initial estimate of the solution>> [~,~,V] = svd(peval(J,v,u0)); v0 = V(:,2:4)*((R*V(:,2:4))\[1;0;0]); % new initial iterate>> [Z,res,fcd] = Newton({g,{ones(4,1),ones(4,1)},{F,K,x,y}}, {gjac,7},{u0,v0}, 1) % rank-7 NewtonStep 0: residual = 6.2e-03Step 1: residual = 1.5e-05 shift = 3.0e-03Step 2: residual = 1.2e-09 shift = 2.4e-05Step 3: residual = 2.2e-16 shift = 1.6e-09Step 4: residual = 4.4e-16 shift = 8.2e-16 shift and residual show both k g xy ( x k , y k ) † rank-7 g ( x k , y k ) k and k g ( x k , y k ) k approach zero, implying the limit (ˆ x , ˆ y ) is both a stationary point anda zero of g . The first component of the output 1 × Z is the computedzero ˆ x = ( , , , )that is accurate to the 16th digits with residual near the unit roundoff. The Jacobian g xy (ˆ x , ˆ y ) is of nullity 1 that is identical to the dimension of the zero (ˆ x , ˆ y ) to theexpanded mapping g . Namely, the ultrasingularity of the mapping f is deflatedby expanding it to g , resulting in a semiregular zero curve in R × R whose firstcomponent is the zero curve of f .At the current stage, theories of the depth-deflation approach (50) for ultrasingularnonisolated zeros are still in development. It is proved in [14] that the Hauenstein-Wampler strong-deflation method terminates in finite many steps. It is highly de-sirable in practical numerical computation for the depth-deflation approach (50) toterminate with the number of deflation steps bounded. The theories in [7, 14] leadto a conjecture: Assume x ∗ is a k -dimensional ultrasingular zero of an analyticmapping f whose Jacobian J ( x ) maintains a constant nullity n > k for all x ∈ Ω ∩ f − ( ) where Ω is an open neighborhood of x ∗ . Then the recursive depth-deflation process (50) terminates in finitely many steps so that x ∗ is a componentof a semiregular k -dimensional zero of the final expanded system.
10 Conclusion
Newton’s iteration is not limited to solving for isolated regular solutions of nonlinearsystems. Nonisolated solutions of f ( x ) = can be solved by the rank- r Newton’siteration x k +1 = x k − J rank- r ( x k ) † f ( x k ) for k = 0 , , . . . that locally quadratically converges to a semiregular zero of positive dimension where J rank- r ( x k ) † is the Moore-Penrose inverse of the rank- r approximation to the Jacobian J ( x k ). When the system is given with empirical data, the rank- r Newton’s iterationstill converges locally to a stationary point that accurately approximates a nonisolatedsolution of the underlying exact system with an accuracy in the same order of thedata, even if the solution of the given system is substantially altered by the dataperturbation.
References [1] J. Backelin. Square multiples n give infinitely many cyclic n-roots. Reports,Matematiska Institutionen 8, Stockholms universitet, 1989.322] D. J. Bates, A. J. Sommese, J. D. Hauenstein, and C. W. Wampler.
NumericalSolving Polynomial Systems with Bertini . SIAM, Philadelphia, 2013.[3] A. Ben-Israel. A Newton-Raphson method for the solution of systems of equa-tions.
J. Math. Anal. Appl. , 15:243–252, 1966.[4] P. T. Boggs. The convergence of the Ben-Israel iteration for nonlinear leastsquares problems.
Math. Comp. , 30:512–522, 1976.[5] X. Chen, M. Z. Nashed, and L. Qi. Convergence of Newton’s method for sin-gular smooth and nonsmooth equations using adaptive outer inverses.
SIAM J.Optim. , 7:445–462, 1997.[6] M. T. Chu. On a numerical treatment for the curve-tracing of the homotopymethod.
Numer. Math. , 42:323–329, 1983.[7] B. Dayton, T.-Y. Li, and Z. Zeng. Multiple zeros of nonlinear systems.
Math-ematics of Computation , 80:2143–2168, 2011. DOI. 10.1090/S0025-5718-2011-02462-2.[8] D. W. Decker, H. B. Keller, and C. T. Kelley. Convergence rate for Newton’smethod at singular points.
SIAM J. Numer. Anal , pages 296–314, 1983. DOI.10.1137/0720020.[9] J.-P. Dedieu and M.-H. Kim. Newton’s method for analytic systems of equationswith constant rank derivatives.
J. Complexity , 18:187–209, 2002.[10] P. Deuflhard. A short history of Newton’s method.
Documenta Mathematica ,Extra Volume:25–30, 2012.[11] J. E. Dennis and R. B. Schnabel.
Numerical Methods for Unconstrained Op-timization and Nonlinear Equations . Prentice-Hall Series in ComputationalMathematics, Prentice-Hall, Englewood Cliffs, New Jersey, 1983.[12] H. F¨uhr and Z. Rzeszotnik. On biunimodular vectors for unitary matrices.
Linear Algebra and its Applications , 484:86–129, 2015.[13] G. H. Golub and C. F. Van Loan.
Matrix Computations . The John HopkinsUniversity Press, Baltimore and London, 4th edition, 2013.[14] J. D. Hauenstein and C. W. Wampler. Isosingular sets and deflation.
Found.of Comput. Math. , 13:371–403, 2013. DOI: 10.1007/s10208-013-9147-y.[15] H. B. Keller. Geometrically isolated nonisolated solutions and their approxima-tion.
SIAM J. Numer. Anal. , 18:822–838, 1981. DOI. 10.1137/0718056.[16] N. Kollerstrom. Thomas Simpson and ‘Newton’s method of approximation’: anenduring myth.
Brit. J. Hist. Sci. , 25:347–354, 2012.3317] T.-L. Lee, T.-Y. Li, and Z. Zeng. A rank-revealing method with updating,downdating and applications, Part II.
SIAM J. Matrix Anal. Appl. , 31:503–525,2009. DOI. 10.1137/07068179X.[18] Y. Levin and A. Ben-Israel. A Newton’s method for systems of m equations in n variables. Nonlinear Anal. , 47:1961–1971, 2001.[19] A. Leykin, J. Verschelde, and A. Zhao. Newton’s method with deflation forisolated singularities of polynomial systems.
Theoretical Computer Science ,pages 111–122, 2006.[20] T.-Y. Li and Z. Zeng. A rank-revealing method with updating, downdat-ing and applications.
SIAM J. Matrix Anal. Appl. , 26:918–946, 2005. DOI.10.1137/S0895479803435282.[21] M. Z. Nashed and X. Chen. Convergence of Newton-like methods for singularoperator equations using outer inverses.
Numer. Math. , 66:235–257, 1993.[22] J. M. Ortega and W. C. Rheinboldt.
Iterative Solution of Nonlinear Equations inSeveral Variables . SIAM Publ., Philadelphia, 2000. First published by AcademicPress, New York and London, 1977.[23] V. Y. Pan. Solving polynomial equations: Some history and recent progress.
SIAM Review , 39:187–220, 1997.[24] A. J. Sommese and C. W. Wampler.
The Numerical Solution of Systems ofPolynomials . World Scientific Pub., Hackensack, NJ, 2005.[25] G. W. Stewart. On the perturbation of pseudo-inverses, projections, and linearleast squares problems.
SIAM Review , 19:634–662, 1977.[26] G. W. Stewart.
Matrix Algorithms, Volume I, Basic Decompositions . SIAM,Philadelphia, 1998.[27] G. W. Stewart and J. Sun.
Matrix Perturbation Theory . Academic Press, Inc,Boston, San Diego, New York, London, Sydney, Tokyo, Toronto, 1990.[28] K. Tenabe. Continuous Newton-Raphson method for solving an underdeter-mined system of nonlinear equations.
Nonlinear Analysis, Theory, Methods andApplications , 3:495–503, 1979.[29] C. W. Wampler and A. J. Sommese. Numerical algebraic geometryand algebraic kinematics.
Acta Numerica , 20:469–567, 2011. DOI.10.1017/S0962492911000067.[30] P.-˚A. Wedin. Perturbation bounds in connection with singular value decompo-sition.
BIT , 12:99–111, 1972.[31] T. J. Ypma. Finding a multiple zero by transformations and Newton-like meth-ods.
SIAM Review , 25:365–378, 1983.3432] Z. Zeng. The approximate irreducible factorization of a univariate polynomial.Revisited. Proc. of ISSAC ’09, ACM Press, pp. 367–374, 2009.[33] Z. Zeng. The numerical greatest common divisor of univariate polynomials. InJ. R. L. Gurvits, P. P´ebay and D. Thompson, editors,
Contemporary Mathemat-ics Vol. 556, Amer. Math. Society,
Randomization, Relaxation and Complexityin Polynomial Equation Solving, pages 187–217, Providence, RI, 2011.[34] Z. Zeng. Sensitivity and computation of a defective eigenvalue.
SIAM J. MatrixAnalysis and Applications , 37(2):798–817, 2016. DOI. 10.1137/15M1016266.[35] Z. Zeng. On the sensitivity of singular and ill-conditioned linear systems.