D Burago
email: burago@math.psu.edu

D Grigoriev

A Slissenko
email: slissenko@univ-paris12.fr

Approximating shortest path for the skew lines problem in time doubly logarithmic in 1/epsilon

We consider two 3-dimensional situations when a polytime algorithm for approximating a shortest path can be constructed. The main part of the paper treats a well known problem of constructing a shortest path touching lines in R 3 : given a list of straight lines L = (L 1 , . . . , L n ) in R 3 and two points s and t , find a shortest path that, starting from s , touches the lines L i in the given order and ends at t . We remark that such a shortest path is unique. We show that it can be length-position ε-approximated (i. e. both its length and its position can be found approximately) in time Rn d α 16

+ O n 2 log log 1 ε , where d is the minimal distance between consecutive lines of L, α is the minimum of sines of angles between consecutive lines, and R is the radius of a ball where the initial approximation can be placed (such a radius can be easily computed from the initial data).

As computational model we take real RAM extended by square and cubic roots extraction. This problem of constructing a shortest path touching lines is known for quite some time to be a challenging problem. The existing methods for approximating shortest paths based on adding Steiner points which form a grid and subsequently applying Dijkstra's algorithm for finding a shortest path in the grid, provide a complexity bound which depends polynomially on 1 ε , while our algorithm for the problem under consideration has complexity linear in log log 1 ε . Our algorithm is motivated by the observation that the shortest path in question is a geodesic in a certain length space of nonpositive curvature (in the sense of A. D. Alexandrov), and it relies on the (elementary) theory of CAT(0)-spaces.

In the second part of the paper we analyze very simple grid approximations. We assume that a parameter a > 0 describing separability of obstacles is given and Preprint submitted to Elsevier Science the part of a grid with mesh size a outside the obstacles is built (for semi-algebraic obstacles all these pre-calculations are polytime). We show that there is an algorithm of time complexity O1 a 6 which, given a-separated obstacles in a unit cube, finds a path (between given vertices s and t of the grid) whose length is bounded from above by (84π * + 96a), where π * is the length of a shortest path. On the other hand, as we show by an example, one cannot approximate the length of a shortest path better than 7π * if one uses only grid polygons (constructed only from grid edges). For semi-algebraic obstacles our computational model is bitwise. For a general type of obstacles the model is bitwise modulo constructing the part of the grid admissible for our paths. Observe that the existing methods for approximating shortest paths are not directly applicable for semi-algebraic obstacles since they usually place the Steiner points forming a grid on the edges of polyhedral obstacles.

Introduction

We consider two 3-dimensional situations when a polytime algorithm for approximating shortest path can be constructed. The main part of the paper concerns a known problem of constructing a shortest path touching lines in R3 in a specified order : given a list of straight lines L = (L 1 , . . . , L n ) in 3dimensional space and two points s and t , find a shortest path that starts from s , touches the lines L i in the given order and ends at t . We will call this problem the skew lines problem (cf. [START_REF] Papadimitriou | An algorithm for shortest-path motion in three dimensions[END_REF]). The second situation is, in a way, more simple: we look for a length approximation of a shortest path (i. e. for a path such that its length, but not necessarily its position, is close to the length of shortest paths) amidst obstacles that are separated.

Related Results

The general problem of constructing a shortest path or even a 'sufficiently good approximation' for such a path is well known [START_REF] Mitchell | Geometric shortest paths and network optimization[END_REF]. It was shown to be NPhard even for convex polyhedral obstacles [START_REF] Canny | New lower bound technique for robot motion planning problemms[END_REF] (by "polyhedral obstacles" we mean unions of polyhedra). The skew lines problem in R 3 is also well known. This problem is mentioned in [START_REF] Papadimitriou | An algorithm for shortest-path motion in three dimensions[END_REF] as, presumably, representing the basic difficulties in constructing shortest paths in 3-dimensional Euclidean space. The paper [START_REF] Papadimitriou | An algorithm for shortest-path motion in three dimensions[END_REF] states "For example, there is no efficient algorithm known for finding the shortest path touching several such lines in a specified order", where "such lines" means "skew lines in three-dimensional space". The same problem of finding a shortest path touching straight lines was mentioned in a survey [START_REF] Schwartz | Algorithmic motion planning in robotics[END_REF] as presumably difficult, though a possibility of using numerical methods was mentioned without any elaboration.

In the problem of approximation of shortest paths we distinguish two types of approximations: length approximation and length-position approximation.

Length approximation means that, given an ε, we look for a path whose length is ε-close to the length of shortest path (ε-close may mean either additive or multiplicative error). Position approximation presumes that a path we wish to construct is in a given neighborhood of a shortest path. An ε-neighborhood of a path Π is a union of balls of radius ε centered at points of Π. To construct a length-position approximation means that, given an ε, we look for a path that is in an ε-neighborhood of a shortest path and whose length is ε-close to the length of shortest path.

In [START_REF] Canny | New lower bound technique for robot motion planning problemms[END_REF] one can find the following result that can be related to the complexity of position approximation: for the case of polyhedral obstacles, finding the sequence of edges touched by a shortest path is NP-hard. In the same article [START_REF] Canny | New lower bound technique for robot motion planning problemms[END_REF] it is proven that for polyhedral obstacles with the size of description N , determining O( √ N ) bits of the length of a shortest path is also NP-hard.

Concerning weaker approximations, for the case of polyhedral obstacles, paper [START_REF] Papadimitriou | An algorithm for shortest-path motion in three dimensions[END_REF] describes an algorithm that finds a path which is of length at most (1 + ε) times of the length of shortest path, and whose running time is polynomial in the size of the input data and 1 ε . Considerable gaps in this paper have been fixed in [START_REF] Choi | Approximate Euclidean shortest path in 3-space[END_REF] and the complexity of the algorithm [START_REF] Choi | Approximate Euclidean shortest path in 3-space[END_REF] is roughly O( n 4 N ε ), where n is the total number of edges in the polyhedra. Further, another method was exhibited in [C87] with complexity bound depending on n and on ε as O( n 2 ε 4 ). An algorithm with complexity bound O( n 4 ε 6 ) is designed in [H99] which in addition, constructs for a fixed source s a preprocessed data structure allowing one to compute subsequently an approximation for any target t with the complexity bound O(log( n ε )).

For the easier problem of constructing a shortest path on a given polyhedral surface one can find a shortest path exactly. First, an algorithm with complexity O(n 3 log n) was exhibited in [START_REF] Sharir | On shortest paths in polyhedral spaces[END_REF] in case of a convex surface, thereupon it was generalized in [START_REF] Mitchell | The discrete geodesic problem[END_REF] to arbitrary polyhedral surfaces, and its complexity is O(n 2 log n). Finally the latter algorithm was improved in [START_REF] Chen | Shortest paths on a polyhedron[END_REF] providing a complexity O(n 2 ). To improve the quadratic complexity, several algorithms which give ε-approximations of a shortest path were produced: in [START_REF] Agarwal | Approximate shortest paths on a convex polytope in three dimensions[END_REF] with complexity bound O(n + 1 ε 3 ) in case of a convex surface, then in [START_REF] Alexandrov | An ǫ-approximation algorithm for weighted shortest paths on polyhedral surfaces[END_REF] with complexity roughly O( n ε ) in case of a simple polyhedral surface. In [H99] for a fixed s , a preprocessed data structure is constructed with complexity O( n ε 3 ) in case of a convex surface and with complexity O(n 2 + n ε ) in case of an arbitrary polyhedral surface, respectively, which allows one for any target t to output an approximation within complexity O(log( n ε )).

All the mentioned approximation results use an approach that places Steiner points on the edges of the polyhedral obstacles to form a grid and after that to find a shortest path in this grid applying the algorithm of Dijkstra [START_REF] Aho | The Design and Analysis of Computer Algorithms[END_REF].

Since the number of the placed Steiner points depends polynomially on n and on 1 ε , the complexity of this approach should depend on n and on 1 ε as well.

The approach of [START_REF] Mitchell | The weighted region problem: finding shortest paths through a weighted planar subdivision[END_REF] to solve the weighted region problem gives an algorithm whose running time depends polynomially on log 1 ε . Combining the binary search with the local optimization as in [START_REF] Mitchell | The weighted region problem: finding shortest paths through a weighted planar subdivision[END_REF] one may probably find an approximation for the skew lines problem; however, in any case, our algorithm is exponentially faster in 1 ε (and our text is shorter than that of [START_REF] Mitchell | The weighted region problem: finding shortest paths through a weighted planar subdivision[END_REF]).

On Computational Models

In the present paper we also study the problem of approximating shortest paths in R 3 trying to improve the computational complexity in some aspect. Even in simple 2-dimensional situations, exact constructions of shortest paths involve substantial extensions of the computational model by operations like square root extraction or more powerful ones, see e. g. [START_REF] Heintz | Search for shortest path around semialgebraic obstacles in the plane[END_REF][START_REF] Grigoriev | Computing minimum-link path in a homotopy class amidst semi-algebraic obstacles in the plane[END_REF][START_REF] Grigoriev | Polytime algorithm for the shortest path in a homotopy class amidst semi-algebraic obstacles in the plane[END_REF]. Note that a shortest path in R 3 between two rational points, with even one cylindrical obstacle, may have transcendental length and necessitate using transcendental points in its construction [START_REF] Heintz | Une borne inférieure pour la construction de chemins polygonaux dans r n[END_REF]. The algorithm we propose for length-position approximation of the shortest path touching lines is based on methods of steepest descent. Each iteration of such a method involves measuring distances between points and other algebraic but not rational operations. If one uses a bitwise computational model, the error analysis becomes difficult and demands an extra considerable work. So to give a solution to the skew lines problem we exploit here an "extended" computational model over reals, which we make precise below.

As it was mentioned above, in this paper we present two results: the first one treats the skew lines problem, where we give a length-position approximation. The second result is an analysis of the complexity of a length approximation of a simple grid method to find a path amidst separated semi-algebraic obstacles. Though the obtained estimation shows that the quality of approximation is rather low, the simplicity of the method makes it worthy of analysis.

For our main result, concerning a length-position approximation of the shortest path touching lines in R3 our computational model is real RAMs [START_REF] Blum | On a theory of computation and complexity over real numbers: NP-completeness, recursive functions and machines[END_REF] with square and cubic roots extraction. This model extends (by allowing also extracting cubic roots) the common model used (implicitly) in [C87,SS86,MMP87,CH90,AHSV97,ALMS98,H99] which admits rational operations and extracting square roots. We mention also that in [Pap85,CSY94] a bitwise model was used which takes into account the bit complexity.

Our second result uses the bitwise model modulo constructing a grid; the latter can be done in polytime for semi-algebraic obstacles.

Our Results and Methods We Use

We start with a remark that the shortest path Π * that we seek, is unique. We show that this path Π * can be length-position ε-approximated in time Rn d α

16

+O(n 2 log log 1 ε ), where d is the minimal distance between consecutive lines of L, α is the minimum of sines of angles between consecutive lines, and R is the radius of a ball where the initial approximation can be placed (we have this formula under the condition that d ≤ 1; it does not diminish the generality; without this condition one must take maximum of ( 35) and (39)). Such a satisfactory radius R can be easily computed from the initial data: R ∈ [

√ n|Π * |, 2n
The algorithm is based on the following geometric idea, which shows that there is a unique minimum for the problem in question, and it is a strict minimum (the latter is quantified by means of a convexity estimate). Given a sequence of lines, one can consider a sequence of copies of Euclidean space, and glue them into a "chain" by attaching "neighboring" spaces along the corresponding line. The resulting (singular) space happens to be of nonpositive curvature (see [START_REF] Ballmann | Lectures on Spaces of Nonpositive Curvature[END_REF][START_REF] Berestovsij | Multidimensional generalized Riemann spaces[END_REF] and the subsection 2.1 below). Now a shortest path we want to construct is a geodesic in this new space, and this immediately implies its uniqueness (by the Cartan-Hadamard Theorem [START_REF] Ballmann | Lectures on Spaces of Nonpositive Curvature[END_REF]). Furthermore, convexity comparisons for the distance functions in non-negatively curved spaces allow us to estimate the rate of convergence in gradient descent methods. This approach is somewhat similar to applications of Alexandrov geometry to certain problems concerning hard ball gas models and other semi-dispersing billiard systems, see [START_REF] Burago | Uniform estimates on the number of collisions in semi-dispersing billiards[END_REF], [START_REF] Burago | Hard Balls Gas and Alexandrov Spaces of Curvature Bounded Above[END_REF]. To avoid excessive use of Alexandrov geometry, we formulate the mentioned convexity property in terms of the Hessian of the distance function and prove it by direct computation.

The exposed geometric idea helps us to achieve the principal feature of our algorithm, namely that its complexity depends linearly on log log( 1 ε ) rather than polynomially on 1 ε as in the existing algorithms from [Pap85,CSY94,C87,AHSV97,ALMS98,H99] based on the grid method (at the expense of a worse dependency on n and on other geometric parameters). It would be interesting to describe more situations when a better dependency of the complexity on 1 ε (logarithmic, cf. [START_REF] Mitchell | The weighted region problem: finding shortest paths through a weighted planar subdivision[END_REF], or even double logarithmic) would be possible simultaneously with better dependency on geometric parameters.

However, our algorithm does not allow to solve the skew lines problem for bitwise models of computations. The latter would require, in particular, estimating the bitsize of the output data. For bitwise models one can consider the following two settings.

(SkewLines Bitwise exact) Given a list of n lines and two points s and t in R 3 , find n algebraic numbers on the consecutive lines such that the path going through these points is the shortest path between s and t touching the lines.

(SkewLines Bitwise approximate) Given a list of n lines and two points s and t in R 3 , and given an ε, find an ε-approximation to the shortest path between s and t touching the lines.

We believe that our result is a step towards solving (SkewLines Bitwise Approximate) problem for which it would be interesting to design an algorithm with complexity polynomial in log(1/ǫ).

Very likely the method we use can be generalized to find a shortest path touching cylinders in R 3 in a prescribed order.

It is straightforward (a rather general case of 3-dimensional TSP with Neighborhoods [START_REF] Mitchell | Geometric shortest paths and network optimization[END_REF]) that if the order of touching lines is not prescribed the problem of finding a shortest path becomes N P -hard. Indeed, take the lines to be parallel and place the points s and t at a plane orthogonal to these lines. Then the problem of finding a shortest path touching all the lines is equivalent to the Euclidean traveling salesman problem in which it is necessary to find a shortest path between s and t and passing all intersecting points of the lines with the plane and which is known to be NP-complete [START_REF] Papadimitriou | Euclidean TSP is NP-complete[END_REF]. We can slightly incline lines if we wish not to have parallel ones.

Our second result is an analysis of very simple grid approximations. The question concerning possibilities of 'generalized' grid approximations is mentioned in the conclusion of [START_REF] Choi | Approximate Euclidean shortest path in 3-space[END_REF]. The matter is that the result of [START_REF] Papadimitriou | An algorithm for shortest-path motion in three dimensions[END_REF][START_REF] Choi | Approximate Euclidean shortest path in 3-space[END_REF] is based on such a method. The initial idea used in [START_REF] Papadimitriou | An algorithm for shortest-path motion in three dimensions[END_REF][START_REF] Choi | Approximate Euclidean shortest path in 3-space[END_REF] is straightforward: partition the edges of polyhedra that constitute the obstacles into sufficiently small segments and take them as vertices of a graph of visibility. Then determine whether two such segments are mutually visible, in which case connect them by edges. It is not a simple task to implement correctly this idea for the bitwise computational model because of the necessity to 'approximate approximations', and this was the main source of gaps in [START_REF] Papadimitriou | An algorithm for shortest-path motion in three dimensions[END_REF]. The primitive grid method we use avoids such difficulties even for semi-algebraic obstacles. But one cannot get approximations that are, in the worst case, better than 7π * , where π * is the length of the shortest path. More precisely, our result is as follows.

When considering grid approximations we assume that a parameter a > 0 describing a separability of obstacles and the part of a grid admissible for the paths is given, and s and t are vertices of the grid lying in the space admissible for the paths. We show that there is an algorithm that in the bitwise machine model has the running time O 1 a 6 and that for given aseparated obstacles in a unit cube finds a path (between given points s and t ) whose length is bounded from above as (84π * +96a), where π * is the length of shortest paths. On the other hand, we show by an example that one cannot approximate the length of shortest paths better than 7π * . We conjecture that 7 is the exact bound for the 3-dimensional case for the method we consider. We do not discuss in detail how to construct the parameter a and the part of the grid admissible for the paths. It depends on the type of obstacles. For semi-algebraic obstacles in R 3 it can be done (precisely, in terms of algebraic numbers) in polytime by a bitwise machine. After that, the construction of a shortest path approximation deals only with natural numbers of a modest size.

As we mentioned above the methods of placing Steiner points and forming a grid developed in [Pap85,CSY94,C87,AHSV97,ALMS98,H99] cannot be directly applied to semi-algebraic obstacles since they place Steiner points just on the edges of polyhedral obstacles. On the other hand, in our simple grid method the dependency on the complexity of obstacles is reduced to the construction of the grid; after that the complexity of the algorithm depends only on the parameter a of the obstacles; this is in accordance with the proposal expressed in [START_REF] Alexandrov | An ǫ-approximation algorithm for weighted shortest paths on polyhedral surfaces[END_REF]: "...while studying the performance of geometric algorithms, geometric parameters (e.g. fatness, density, aspect ratio, longest, closest) should not be ignored...".

Structure of the Paper

The paper is divided into two parts. Its structure is as follows.

The first and the main part of the paper deals with shortest paths touching lines. In section 2 we give basic properties of the spaces under consideration, in particular, we explain what geometric properties ensure the uniqueness of the shortest path. Our main theorem is formulated in section 2.5. In fact, this part is not needed for the proof as we reduce the problem to a problem of computing a minimum of a strictly convex function of many variables on a compact space. However, these geometric considerations, that are easy to apply, show the road to take, and they are very instructive indeed for the problem of shortest path that has rather various technical contexts and rarely has a unique solution. Section 3 contains main technical estimates. The central point is to bound from below the eigenvalues of the Hessian of the path length function (section 3.2, Proposition 2). This bound is crucial for the estimation of complexity of the gradient descent. The latter provides an initial approximation for the application of Newton's method that follows the gradient descent. Section 4 gives our approximation algorithm. It starts with a construction of an initial approximation, then applies a gradient descent and after that Newton's method. In the last section 5 of Part 1 we estimate the complexity of the algorithm.

Part 2 presents our grid algorithm. The definition of separability is given in section 6. The same section (Proposition 7) gives a lower bound for the probability to get well separated obstacles when randomly choosing balls as obstacles. Then in section 7 we present the algorithm that constructs a length approximation.

Part 1: Shortest Paths Touching Lines 2 Basic Properties of Shortest Paths Touching Lines

Our construction is based on the following initial observation: we are seeking a shortest path in a simply connected metric space of nonpositive curvature. This observation guarantees the uniqueness of the shortest path; moreover, it suggests that the distance functions in the space enjoys certain convexity properties, which will actually permit us to use the method of steepest descent to reach an approximate position of the shortest path.

As we mentioned in the Introduction, formally speaking we can omit these considerations and go directly to the estimations of convexity of the length function within our parametrization of paths. However, these geometric arguments are simple and direct, and can be used in other problems of construction of shortest paths. We remark that usually the uniqueness of the shortest path does not take place though the metric function is quite good (e. g. a sphere as an obstacle in R 3 ).

First, we recall some known facts about spaces of nonpositive curvature.

On Spaces of Nonpositive Curvature

Let M be a simply connected metric space. Denote by |XY | M , or by |XY | if M is clear from the context, the distance between points X and Y in M. A path in M between two points A and B is a continuous mapping ϕ of [0, 1] into M such that ϕ(0) = A and ϕ(1) = B. The length of a path ϕ can be defined as the limit of sums 0≤i≤(n-1) |ϕ(i

1 n )ϕ((i + 1) 1 n )| as n → ∞.
A shortest path between two points is a path connecting these points and having the smallest length. We suppose that for any two points of M there exists a shortest path between them.

By the angle between two intervals going out of the same point on the plane R 2 we mean the (non directed) angle in [0, π]. Given 3 points A, B and C the angle between intervals AB and AC will be denoted by ∠BAC or ∠CAB, the angle between 2 intervals ξ and ζ will be denoted by ∠ξζ, and the angle between 2 vectors will be treated as the angle between two intervals obtained from the vectors by translating them to the same origin. If different spaces are considered in the same context the angle will be denoted by ∠ M . Take any three points A, B and C in M (Figure 1). Let Π AB and Π AC be shortest paths respectively between A and B, and A and C. On these paths we take two arbitrary points Π AB (t) and Π AC (τ ), where t and τ are the values of parameters determining these points.

A B C Π ΑC (τ) B' C' M (BAC) Π ΑΒ (t) B' t C' τ
Then we choose any point A ′ in R 2 and take in R 2 two points B ′ t and C ′ τ that satisfy the equalities:

|A ′ B ′ t | R 2 = |Π AB ([0, t])| M , |A ′ C ′ τ | R 2 = |Π AC ([0, τ ])| M , |B ′ t C ′ τ | R 2 = |Π AB (t)Π AC (τ )| M . We are interested in the triangle △A ′ B ′ t C ′ τ in the plane R 2 .
Denote by θ(t, τ ) the angle between the intervals A ′ B ′ t and A ′ C ′ τ in R 2 . The angle between Π AB and Π AC at A is defined as lim t,τ →0 θ(t, τ ), if the latter exists.

Define the angle at A to B and C, or shorter ∠ M (BAC) as the supremum of angles at A between a shortest path from A to B and a shortest path from A to C.

Consider again three points A, B and C in M. Draw in R 2 a triangle A ′ B ′ C ′ determined by the following conditions:

|A ′ B ′ | R 2 = |AB| M , |A ′ C ′ | R 2 = |AC| M and ∠ M (BAC) = ∠ R 2 (B ′ A ′ C ′ ).
By definition, the space M is of nonpositive curvature if the angle ∠ M (BAC) exists and

|B ′ C ′ | R 2 ≤ |BC| M for any three points A, B and C of M.
It is known that (see, for instance, [START_REF] Ballmann | Lectures on Spaces of Nonpositive Curvature[END_REF]) (Uniqueness of the Shortest Path.) In any simply connected space of nonpositive curvature, there is a unique shortest path between any two points.

The Configuration Space

Let L = (L 1 , . . . , L n ) be a list of straight lines in R 3 , and s and t two points. Without loss of generality, we assume that these points are not on the mentioned lines. Clearly, the shortest paths touching the lines of L in the given order are polygonal chains (broken lines).

The space where we are looking for such a shortest broken line can be obtained by the following operation. For each pair (L i , L i+1 ) take a copy R i of R 3 . Glue consecutively R i and R i+1 along the line L i+1 . Denote the obtained space by R L . Each space R i is a space of nonpositive curvature, and they are glued along isometric convex sets. The Reshetnyak's theorem (see [START_REF] Berestovsij | Multidimensional generalized Riemann spaces[END_REF]) implies that the resulting space R L is of nonpositive curvature.

The next two paragraphs are not necessary to describe the main algorithm, but rather to clarify some geometric background behind it. Similarly to non-positively curved spaces, one can consider spaces with other curvature bounds (say, spaces with curvature bounded above by -1). They feature even stronger convexity of distance functions, and algorithms for constructing shortest paths converge in such spaces even faster. Shortest paths in such spaces have remarkable trapping properties, as it is asserted by the classic Morse Lemma. The spaces constructed in this paper by gluing several copies of Euclidean space along a number of lines are non-positively curved. Of course, they are not spaces of curvature bounded above by any negative number, for they have flat parts (regions in the Euclidean spaces used for this construction). However, if one considers two points A, B in different copies of Euclidean space and such that the shortest path between these points meets al least two lines used to glue the spaces together, and the gluing lines are not parallel, then a distance function |A • | restricted to a neighborhood of B has the same convexity properties as if the space was of curvature bounded above by a negative constant. This allows us to treat such shortest paths as if they were in a space with negative upper curvature bound.

Note that it is easy to see geometrically where this negative curvature comes from if the gluing lines are skew lines, and what happens if the gluing lines are parallel. First assume that the gluing lines are parallel. Consider the following plane lying in our space and passing through A and B. It is made of three parts: a half-plane containing A and bounded by the first gluing line (this part lies in the same copy of Euclidean spaces as A), then a strip between the two gluing lines (it lies in the copy of Euclidean space which the shortest path connecting A and B passes through), and finally a half-plane containing B and bounded by the second gluing line. These three parts indeed form a plane (a surface isometric to a plane) lying in our space totally geodesically. It contains A and B, and hence the convexity properties of the distance function |A • | restricted to a neighborhood of B are no better than in a flat space. To see where the negative curvature comes from in case of skew lines, consider a short Euclidean segment Γ containing B and parallel to the second gluing line. Connect all points of Γ with A. If the gluing lines were parallel, we would get just a part of the plane discussed above (more precisely, a triangular region). However, in case of skew lines the family of segments connecting A with Γ form a non-flat surface. This surface has two flat parts (a triangle with a vertex at A, and a trapezoid (one of the bases of the trapezoid is Γ, and the other one lies in the gluing line). However, the third part of this surface is a ruled surface connecting two segments on the gluing lines, and this is indeed a negatively curved surface.

A similar construction can be used to obtain a space of nonpositive curvature to seek shortest paths consecutively touching convex cylinders Z = (ζ 1 , . . . , ζ n ) in R N . In this case the gluings are made along the corresponding cylinders.

Consider the shortest path problem in R L , where two points s and t are fixed.

From the uniqueness of shortest path between two points in a simply connected nonpositively curved space we immediately conclude that:

Claim 1 There is a unique shortest path between any two given points in R L ; hence a shortest path touching the lines in a prescribed order and connecting two given points is unique.

Furthermore, a standard distance comparison argument with Euclidean development of a broken line representing an approximate solution implies that length-approximation guarantees position-approximation (this is proven in detail in Lemma 4, section 4.1):

Claim 2 If the length of the actual shortest path is L, every path whose length is less than L+ǫ belongs to the Lǫ/2-neighborhood of the actual shortest path.

To simplify our constructions we will assume that (see Assumption 2.5 in section 2.5) the consecutive lines of L are pairwise disjoint. The case with intersecting consecutive lines needs a particular attention when our approximation to the shortest path is close to a point of intersection of lines. By a path in R L we mean a polygon P consisting of n + 1 links (straightline intervals), denoted by P i , 1 ≤ i ≤ n, connecting s with t via consecutive lines, i. e. the link P 1 connects s with L 1 , the link P i+1 connects L i with L i+1 (1 ≤ i ≤ n -1), and the link P n+1 connects L n with t . When speaking about the order of points on a path we mean the order corresponding to going along the path starting from s . Thus any link of a path is directed from s to t .

Notice that for a shortest path P the angle between P i (incident ray) and L i must be equal to the angle between L i and P i+1 (emergent ray). Clearly, the incident angle being fixed, the emergent rays constitute a cone (which may have up to two intersection points with L i+1 -this means that a geodesic in our configuration space can branch, having two different continuations after passing through a line L i ).

Technical Notations

Throughout the text we use some notations that we summarize here -the reader can consult the list when coming across a new notation. The list is divided into 4 parts: notations concerning vectors, lengths, angles, paths. Some more local notation will appear later.

Notations 1 (Vectors, Lines, Points):

• U, V is the scalar product of vectors U and V .

• |U | or U is the L 2 -norm of a vector U (the length of U ).

• A , where A is a matrix, is the spectral norm of A. For positive definite symmetric matrices (our main case) A = sup{ AU, U : U = 1}.

• L = (L 1 , . . . , L n ) is the list of straight lines in R 3 we consider, and s and t are respectively the first and the last point to connect by a polygon touching the lines of L in the order of the list.

• Λ 0 i is a fixed vector determining a point on L i . We use it as the origin of coordinates on L i .

• Λ i is a unit vector defining the direction on L i .

• T = (t 1 , . . . , t n ) is a list of n reals that will serve as coordinates respectively on L 1 , . . . , L n . Our applications of gradient and Newton's methods take place in the space R n of such points T .

• To make the notations more uniform we assume that Λ 0 0 = s , Λ 0 n+1 = t , and that always t 0 = t n+1 = 0.

• W i (t i )= df (Λ 0 i + t i Λ i )
; clearly, it is a point on L i determined by a single real parameter t i . Note that W 0 (t 0 ) = s and W n+1 (t n+1 ) = t .

• V i = V i (t i , t i )= df W i+1 (t i+1 ) -W i (t i ) is a vector connecting W i with W i+1 (in this order), 0 ≤ i ≤ n.
• D 0 is a vector from s to L 1 perpendicular to L 1 , D n is a vector from L n to t perpendicular to L n , and D i is a vector connecting L i and L i+1 and perpendicular to both of them.

• For a directed interval σ in R 3 we denote by σ -and σ + its beginning and its end respectively.

Notations 2 (Lengths):

• ρ will be used to denote various L 2 -distances, e. g. ρ(s , L i ) will denote the distance between s and L i .

• v i = df |V i | is the length of V i , 0 ≤ i ≤ n.
• d i is the square of the distance between L i and L i+1 for 1 ≤ i ≤ n -1, i. e.

d i = ρ(L i , L i+1 ) 2 , d 0 = df ρ(s , L 1 ) 2 , d n = df ρ(L n , t ) 2 . • d= df min{d i : 0 ≤ i ≤ n}. • d= df + √ d.
To simplify formulas we assume that d ≤ 1; we do not loose the generality as we can change the coordinates in an appropriate way, and within our model of computation this change does not affect the complexity. In the proofs we will also give formulas without this assumption.

Notations 3 (Angles:)

• ∠ABC, where A, B and C are points, is the angle at B in the triangle determined by these 3 points. In fact, we will consider angles in R 3 , though formally speaking we are in configuration space R L defined in section 2.2.

• ∠V V ′ , where V and V ′ are vectors, is the angle in [0, π] between these vectors.

• α i = df (sin ∠Λ i Λ i+1 ) 2 , 1 ≤ i ≤ n -1, α= df min{α i : 1 ≤ i ≤ n -1}.
• α= df + √ α.

• β i = df (cos ∠Λ i Λ i+1 ) 2 , β i = df cos ∠Λ i Λ i+1 , 1 ≤ i ≤ n -1.
Notations 4 (Paths):

• A path Π(T ) is determined by a list of reals T that gives the consecutive points W (t i ) of Π(T ) on the lines of L i connected by this path. Thus, W 0 (t 0 ) = s , W i (t i ) = Π(T ) ∩ L i for 1 ≤ i ≤ n and W n+1 (t n+1 ) = t , where W i , t i and L i are defined in Notations 1.

• By P i for a path P we denote the ith link of P . We consider P i as an interval directed from L i to L i+1 , 0 ≤ i ≤ n. Regarded as a vector, it coincides with V i .

• Π * is the shortest path between s and t .

• T * = (t * 1 , . . . , t * n ) is the point in R n defining Π * , i. e. Π * = Π(T * ).

• Π 0 is the initial approximation of Π * defined below in section 2.4.

• T 0 is the point in R n defining Π 0 , i. e. Π 0 = Π(T 0 ).

• r= df |Π 0 | is the length of the initial approximation. It will be clear from its construction that d < r. Moreover, without loss of generality, just to simplify some formulas that will appear in the estimations of complexity, we suppose that r > 1.

• R= df r √ n is the radius of a ball in R n that will contain all paths in question.

• B is the closed ball in R n of radius R centered at T 0 ; it contains all paths under consideration.

Initial Approximation for the Shortest Path

Denote by P 0 i the base of the perpendicular from s to L i .

Take as initial approximation to the shortest path the polygon Π 0 that starts at s , then goes to P 0 1 , then to P 0 2 and so on to P 0 n and finally to t . Clearly Π 0 ⊆ B and

|P i P i+1 | ≤ (|s P i | + |s P i+1 |) for 1 ≤ i ≤ n -1.
(1)

The length of Π 0 can be bounded in terms of the length of the shortest path Π * as

|Π * | ≤ |Π 0 | = r ≤ 2 n-1 i=1 ρ(s , L i ) + |s t | ≤ 2n|Π * |. (2)
The first inequality follows from the fact that the length of any path connecting s and t and touching lines from L is greater than or equal to |Π * |, the second inequality is implied by (1), and the third one is due to the fact that the distance from s to any line is not greater than |Π * |.

For R the estimations (2) give

√ n|Π * | ≤ R ≤ 2n 3 2 |Π * |. (3)
We remark that if L i = {W 0 i + tΛ i : t ∈ R}, then P i = W 0 i + t 0 i Λ i is determined by P i , Λ i = 0, and thus,

P i = W 0 i -W 0 i , Λ i Λ i . (4)
The parameter R = |Π 0 | √ n intervenes in our estimation of complexity, and one may ask what is the dependence of R on the length of the input. Suppose that the input lines are represented by W 0 i and Λ i , and the maximum of the lengths of involved numbers is ν. Then the length of the input is O(nν). Clear that in the worst case the value R is exponential in this input length. But it can be as well linear, logarithmic and so on. The precise bound on R is an open question.

Main Theorem

Hereafter we suppose that Assumption: α > 0 and d > 0.

The first inequality says that there are no consecutive parallel lines among L i , and the second one says that the consecutive lines do not intersect and the points s and t are not on the lines.

Now we can formulate our main result:

Main Theorem.

There is an algorithm that finds a length-position ε-approximation of the shortest path touching lines in time

Rn d α 16 + O n 2 log log 1 ε ,
where R, d and α are defined above in Notations 4, 2, 3.

Remark 1. The complexity bound of the Main Theorem depends actually on ε if the latter is sufficiently small (see the end of section 5.2), namely,

ε < O d 4 α 2 r 3 √ n .
Otherwise, the second summand of the estimation which contains ε is dominated by the first summand in the complexity bound. As we remarked in the Introduction, in order to simplify the first summand we supposed, without loss of generality, that d ≤ 1; without this assumption one must take maximum of (35) and (39) which have a similar form.

Bounds on Eigenvalues of the Hessian of the Path Length

This section contains technical bounds on the norm of derivatives of the path function, in particular the main estimation, namely, a lower bound on the eigenvalues of the Hessian ('second derivative') of the path length. The fact that the length function is at "least as convex as in Euclidean space" follows from the curvature bound. We will make direct computation though, due to the fact that we use a particular coordinate system. Consider a path Π(T ) = Π(t 1 , . . . , t n ) represented by points W i , 0 ≤ i ≤ n + 1 (Notations 4 from section 2.3).

Notations 5 :

• f (T )= df |Π| = n i=0 v i is the length of Π(T )
, where v i is defined in Notations 2.

• g(T )= df f ′ (T )= df gradf (T ) = ( ∂f ∂t 1 (T ), . . . , ∂f ∂tn (T )) is its gradient (first derivative).

Let us begin with the formula for g, that is the first variation of length. For 1 ≤ i ≤ n we set

π i = ∂f ∂t i = ∂ ∂t i v i-1 + v i = ∂v i-1 ∂t i + ∂v i ∂t i = V i-1 , Λ i v i-1 - V i , Λ i v i , (5)
and

g = (π 1 (t 1 , t 2 ), π 2 (t 1 , t 2 , t 3 ), . . . , π n-1 (t n-2 , t n-1 , t n ), π n (t n-1 , t n )) (6)
(which stresses the variables on which each component of g depends).

3.1 The Second Variation Formula for the Path Length

The Hessian of |Π|, which we will denote ∆ = ∆(T ) (Jacobian of |Π| ′ ), looks as shown in Figure 2. The matrix of ∆ is symmetric 3-diagonal with positive diagonal entries. Here are the formulas for derivatives involved in ∆ in arbitrary coordinates (for further references).

                                    ∂ 2 v 0 ∂t 2 1 + ∂ 2 v 1 ∂t 2 1 ∂ 2 v 1 ∂t 1 ∂t 2 0 . . . 0 0 ∂ 2 v 1 ∂t 1 ∂t 2 ∂ 2 v 1 ∂t 2 2 + ∂ 2 v 2 ∂t 2 2 ∂ 2 v 2 ∂t 2 ∂t 3 . . . 0 0 0 ∂ 2 v 2 ∂t 2 ∂t 3 ∂ 2 v 2 ∂t 2 3 + ∂ 2 v 3 ∂t 2 3 . . . 0
∂ 2 v n-1 ∂t n-1 ∂tn ∂ 2 v n-1 ∂t 2 n + ∂ 2 vn ∂t 2 n                                    
∂ 2 v i ∂t i ∂t i+1 = ∂ ∂t i V i , Λ i+1 v i = ∂V i ∂t i , Λ i+1 v i - V i , Λ i+1 V i , ∂V i ∂t i v 3 i = - Λ i , Λ i+1 v i + V i , Λ i+1 V i , Λ i v 3 i = (7) = - 1 v i cos ∠Λ i Λ i+1 -cos ∠V i Λ i+1 • cos ∠V i Λ i (8) ∂ 2 v i ∂t 2 i = 1 v i 1 - V i , Λ i 2 v 2 i = (sin ∠V i , Λ i ) 2 v i . ( 9
) ∂ 2 v i ∂t 2 i+1 = 1 v i 1 - V i , Λ i+1 2 v 2 i = (sin ∠V i , Λ i+1 ) 2 v i . ( 10
)

Lower and Upper Bound on the Eigenvalues of the Hessian of the Path Length

Decompose ∆ into a sum of matrices with 2 × 2-blocks corresponding to links of Π plus two 1 × 1-blocks for the first and the last links.

Notations 6

• ∆ 0 = df ∂ 2 v 0 ∂t 2 1 • ∆ i = df        ∂ 2 v i ∂t 2 i ∂ 2 v i ∂t i ∂t i+1 ∂ 2 v i ∂t i ∂t i+1 ∂ 2 v i ∂t 2 i+1        f or 1 ≤ i ≤ n -1 • ∆ n = df ∂ 2 v n ∂t 2 n
• ∆ i is the n × n-matrix consisting of ∆ i situated at its proper place in ∆ and with zeros at all other places. More precisely, (1, 1)-element of ∆ i is placed at (i, i) for 1 ≤ i ≤ n and at (1, 1) for i = 0. Within these notations:

∆ = 0≤i≤n ∆ i . (11)
Lower Bound on Eigenvalues of ∆.

Strict convexity of metric implies that the second derivative ∆ is positive definite. However, we need constructive upper and lower bounds on eigenvalues of ∆. So we do not use explicitly the just mentioned property of the metricit will follow from the estimations below.

We start with a more difficult question of obtaining a lower bound on the least eigenvalue of ∆, i. e. on inf{ ∆U, U : U = 1}. To do it we reduce this problem to the same problem for matrices ∆ i .

Lemma 1 The only eigenvalue of ∆ 0 and that of ∆ n is greater than d 2r 3 , and both eigenvalues of ∆ i for 1 ≤ i ≤ n -1 are greater than dα 2r 3 . Thus dα 2r 3 is a lower bound on the norm of ∆ i for all i (and for any t i ).

Proof. As supposed in the beginning of section 2.5 there are no consecutive parallel lines among L i . Thus vectors Λ i , Λ i+1 and D i (see notations in section 2.3) constitute a basis in R 3 , and we can represent the link V i as

V i = t i • Λ i + t i+1 • Λ i+1 + D i (
here t i and t i+1 are different from those used in notations in section 2.3).

To simplify computations rewrite this formula as W = xA + yB + D, where For i = 0 we have that the only eigenvalue of ∆ 0 is its element

W = V i , x = t i , A = Λ i , y = t i+1 , B = Λ i+1
∂ 2 v 0 ∂t 2 1 = 1 v 0 1 -W,B 2 v 2 0 .
Here

W, B = yB + D, B = y, W, W = v 2 0 = yB + D, yB + D = y 2 + d 0 and thus ∆ 0 = 1 v 0 1 - y 2 v 2 0 = 1 v 3 0 y 2 + d 0 -y 2 = d 0 v 3 0 > d r 3 . Similarly for i = n we have ∆ n ≥ dn v 3 n > d r 3 .
Consider i ∈ {1, . . . , n -1}. Notice that in our notations

v 2 = xA + yB + D, xA + yB + D = x 2 + 2xy β i + y 2 + d i and ∆ i =        ∂ 2 v ∂x 2 ∂ 2 v ∂x∂y ∂ 2 v ∂x∂y ∂ 2 v ∂y 2       
. Now we compute the second derivative of v:

∂ 2 v ∂x 2 = 1 v 1 - xA + yB + D, A 2 v 2 = 1 v 3 v 2 -(x + y β i ) 2 = 1 v 3 y 2 + d i -(y β i ) 2 = 1 v 3 y 2 α i + d i . ∂ 2 v ∂y 2 = 1 v 1 - xA + yB + D, B 2 v 2 = 1 v 3 v 2 -(x β i + y) 2 = 1 v 3 x 2 α i + d i . ∂ 2 v ∂x∂y = - 1 v i β i - xA + yB + D, B xA + yB + D, A v 2 = - 1 v 3 β i (x 2 + 2xy β i + y 2 + d i ) -(x β i + y)(x + y β i ) = - 1 v 3 (x 2 β i + 2xy β 2 i + y 2 β i + d i β i ) -(x 2 β i + xy + xy β 2 i + y 2 β i ) = - 1 v 3 xy β 2 i + d i β i -xy = - 1 v 3 d i β i -xyα i .
Let us estimate the determinant

det∆ i = 1 v 6 (y 2 α i + d i )(x 2 α i + d i ) -(d i β i -xyα i ) 2 = 1 v 6 x 2 y 2 α 2 i + d i α i (x 2 + y 2 ) + d 2 i -d 2 i β 2 i -x 2 y 2 α 2 i + 2xyα i d i β i = 1 v 6 d i α i (x 2 + y 2 ) + d 2 i α i + 2xyα i d i β i = α i v 6 d i (x 2 + y 2 ) + d 2 i + 2xyd i β i ≥ α i v 6 d i (x 2 + y 2 ) + d 2 i -2|x||y|d i ≥ d 2 i α i v 6 . 20 Note that Tr∆ i = 1 v 3 (α i (x 2 + y 2 ) + 2d i ). The smallest of two eigenvalues of ∆ i is not less than det∆ i Tr∆ i ≥ v 3 d 2 i α i v 6 (α i (x 2 + y 2 ) + 2d i ) ≥ d 2 i α i v 3 2d i = d i α i 2v 3 > dα 2r 3 .
Proposition 2 All eigenvalues of ∆(T ) are greater than dα 2r 3 . Thus ∆(T ) is positive definite, the metric function is strictly convex and ∆(T ) ≥ dα r 3 , and

∆ -1 (T ) ≤ r 3 dα for any T .
Proof. Now we estimate ∆U, U for U = 1. Formula (11) gives

∆U, U = i ∆ i U, U = i ∆ i U i , U i ,
where U i is the appropriate 2-vector. From Lemma 1 we have

∆U, U ≥ i dα 2r 3 U i , U i = dα 2r 3 i U i , U i = dα 2r 3 • 2 = dα r 3 .
This inequality gives a lower bound on eigenvalues of ∆. Hence, ∆ ≥ dα r 3 , and ∆ -1 ≤ r 3 dα .

Upper Bound on the Norm of ∆.

The matrix ∆ is symmetric, therefore, using formulas (8), ( 9), (10), Notations 2 and standard inequalities for the norms one gets

∆ = ∆(T ) ≤ √ n • ∆ 1 = √ n • max i j |∆ i,j | ≤ 4 √ n max i 1 v i ≤ 4 √ n d . (12)
for all T .

Newton's Method for the Search of Zero of the Gradient of Path Length

Our algorithm consists of 3 phases described below in section 4. The third phase is an application of Newton's method to approximate the zero of the Jacobian f ′ of the length function. The main difficulty, however, is to find an appropriate approximation for Newton iterations. That will be done by a gradient descent described in section 4.3. We will use the bounds we just have obtained to estimate the rate of convergence of Newton's method.

First let us recall standard facts about Newton's method. Consider a mapping H = (h 1 , . . . , h n ) from R n into R n . Letter X without or with superscripts and subscripts will denote vector-columns from R n . Newton's iterations are defined by the formula

X k+1 = X k -(H ′ (X k )) -1 H(X k ), ( 13
)
where H ′ is the Jacobian of H.

Proposition 3 gives a sufficient condition to ensure a fast convergence of Newton's iterations (it can be found in textbooks on numerical methods).

Proposition 3 Suppose that X 0 is a zero of H, i. e. H(X 0 ) = 0, where 0 = (0, . . . , 0) ∈ R n .

Let a, a 1 , a 2 be reals such that 0 < a and 0 ≤ a 1 , a 2 < ∞, and denote

Ω a = {X : X -X 0 < a}, c = a 1 a 2 and b = min{a, 1 c }. If X 0 ∈ Ω b and (A) (H ′ (X)) -1 ≤ a 1 for X ∈ Ω a , (B) H(X 1 ) -H(X 2 ) -H ′ (X 2 )(X 1 -X 2 ) ≤ a 2 X 1 -X 2 2 for X 1 , X 2 ∈ Ω a then X k -X 0 ≤ 1 c c X 0 -X 0 2 k . ( 14
)
In Proposition 3 there are two constants a 1 and a 2 : the norm of the inverse of H ′ has to be at most a 1 (condition (A)); and a 2 is, in fact, un upper bound on the norm of the second derivative of H ′′ (condition (B)). Both bounds must hold in a-neighborhood of the zero X 0 of H. In our case a, that appears only in (B), will be 'big' (more precisely, we will take a = R), so we can ignore it for now. Then the convergence is determined by two things: a parameter

b = 1 a 1 a 2
, and the choice of initial approximation X 0 that must be in an open b-neighborhood of the zero X 0 . We are interested in making c X 0 -X 0 smaller than 1 with a known upper bound less than 1. We will construct X 0 such that c X 0 -X 0 ≤ 1 2 .

We will take the bound from Proposition 2 for a 1 . As for a 2 , it is not hard to estimate it using formulas for elements of ∆. We will choose X 0 in section 4.3.

Choosing parameters a 1 and a 2 to satisfy (A) and (B) of Proposition 3.

We start with calculating a 2 . To satisfy Condition (B) it is sufficient to show that the (partial) second derivatives of f are bounded in some neighborhood of X 0 . We use Taylor's formula with the third derivatives of f . Take any T 1 , T 2 ∈ B. In our case

H(X 1 ) -H(X 2 ) -H ′ (X 2 )(X 1 -X 2 ) = f ′ (T 1 ) -f ′ (T 2 ) -∆(T 2 )(T 1 -T 2 ). ( 15
)
We assume that the vectors involved in (15) are represented as columns. To bound from above the norm of this vector, firstly estimate its components. Using notations T j = (t

(j)
1 , . . . , t (j) n ), j = 1, 2, we can write the ith component of (15) as

π i (t (1) i-1 , t (1) i , t (1) i+1 ) -π i (t (2) i-1 , t (2) i , t (2) i+1 ) + j ∂π i (t (2) j ) ∂t j (t (1) j -t (2) j ). ( 16
)
Taylor's formula says that ( 16) is equal to

1 2 j=i-1,i,i+1 ∂ 2 π i (ξ) ∂t j ∂( t ′ j -t j ) 2 + j,k=i-1,i,i+1 ∂ 2 π i (ξ) ∂t j ∂t k (t ′ j -t j )(t ′ k -t k ) (17)
for some vector ξ = (ξ j ) j whose jth component is between t

(1) j and t

(2) j .

The equalities (7), ( 8), ( 9), (10) show that the second derivatives of π i (notation from (6)) involved in (17) are of the form

1 v 2 i (sin O(1) ϕ, cos O(1) ψ, (O(1))) + 1 v 2 i-1 (sin O(1) ϕ, cos O(1) ψ, (O(1))),
where (sin O(1) ϕ, cos O(1) ψ, (O(1)) is a sum of expressions in parenthesis ("arguments of Σ"). Thus the absolute value of each second derivative of π i is bounded by O( 1 d ). Hence, the absolute value of each component of vector ( 15) is bounded by O( 1 d ). This implies that in our case (see ( 15))

H(X 1 ) -H(X 2 ) -H ′ (X 2 )(X 1 -X 2 ) = f ′ (T 1 ) -f ′ (T 2 ) -∆(T 2 )(T 1 -T 2 ) ≤ C 1 √ n d T 1 -T 2 2 (18)
for some constant C 1 > 0 and for T 1 , T 2 ∈ B.

Thus we set

a 2 = C 1 • √ n d . ( 19
)
As in our case H ′ = ∆, Proposition 2 permits to take as a 1 the lower bound for ∆ -1 from this proposition.

Now we can define all the parameters for Newton's method:

• a = R (justified by ( 18)),

•

a 1 = r 3 dα (see Proposition 2), a 2 = C 1 n d (see (18)),
• c = a 1 a 2 = C 1 nr 3 αd 2 for the constant C 1 from (18).

Our assumptions about r (see Notations 4) imply that a = r ≥ 1 c .

Hence we take

b = 1 c = αd 2 C 1 nr 3 = α 2 d 4 C 1 nr 3 . (20)
For further references we rewrite (20 as (we use Notations 2, 3)

• b 2 = c 2 α 2 d 4 nr 3 = c 2 αd 2 nr 3 , (21)
where

c 2 = df 1 2C 1 .
To ensure fast convergence we will need a 'good' initial approximation X 0 for Newton's method, namely such that

X 0 ∈ Ω b 2 . ( 22
)
For such X 0 we have

c X 0 -X 0 ≤ 1 2 , (23)
and the rate of convergence guaranteed by Proposition 2 will be

1 2 2 k
, where k is the number of iterations (see Proposition 3). It remains to notice that we use each

|AA * | = |t i -t * i | twice for 1 ≤ i ≤ n.
Thus, to obtain a length-approximation with precision ε, we will construct a position ε √ n -approximation.

Preliminary computations

Phase 1: Compute the path Π 0 and the values d, d, α, α r, R and save also all intermediate values that appear within these calculations (all that will be used in the estimations that govern the two other phases).

Our construction of the initial approximation Π 0 was described in section 2.4. The values d, d are defined in Notations 2, the values α, α are defined in Notations 3 and the values r, R are defined in Notations 4.

Initial Gradient Descent to Initial Approximation of Newton's Method

Recall that f (t 1 , . . . , t n ) denotes the length of a path Π represented by points on the lines of L, i. e. f (t 1 , . . . , t n ) = f (T ) = |Π(T )| (see Notations 5). This function f is convex (see Proposition 2), and the point at which it attains its minimum is denoted by T * (section 2.3, Notations 4). Proposition 2 gives a positive lower bound for the norm of the Hessian ∆ of f ; we denote this lower bound by δ. Let δ be the upper bound on the norm of ∆ from (12).

Consider g= df gradf = ∂f ∂t 1 , . . . , ∂f ∂tn . Notice that g(T * ) = 0. We have g(T ) = T T * ∆, where one can integrate ∆ along any path from T * to T .

Denote by η the distance from T * to T . Since ∆ is positive definite and its norm is at least δ, we have ∆(V ), V / V ≥ δ V for all vectors V . Hence the projection of the integral of ∆ along the segment from T * to T to this segment is at least δη, and in particular the norm of the integral itself is at least δη.

Recalling that this integral is equal to the gradient of f at T , we conclude that the absolute value ∂f ∂t i (T ) of one of the coordinates of the gradient g(T ) is at least δ

•η √ n . Choose i such that ∂f ∂t i (T ) ≥ δ • η √ n . (24)
Let us estimate how one can decrease the length of the path by changing t i (which geometrically represents moving the point where the path meets L i ).

To study how f depends on t i , introduce a function in one variable τ given by ϕ(τ ) = f (t 1 , . . . , t i-1 , τ, t i+1 , . . . , t n ). Function ϕ is the restriction of f to a straight line and hence convex. Let ϕ attain its minimum at τ 0 .

We summarize the just introduced notations for further references:

Notations 7 :

• δ = 4 √ n d
is an upper bound for ∆ (formula (12)).

• δ = dα r 3 is a lower bound for ∆ (Proposition 2).

• g= df gradf = ∂f ∂t 1 , . . . , ∂f ∂t n .

• η is the distance from T * to T , i. e. η= df T * -T , where T * is a point where the length of our path Π(T * ) is minimal and T is the current point in the space of parameters determining our path.

• ϕ(τ ) = f (t 1 , . . . , t i-1 , τ, t i+1 , . . . , t n ) (f is defined in Notations 5).

• i is chosen such that (24).

We wish to estimate |ϕ(t i ) -ϕ(τ 0 )| from below. Notice that our bounds on the Hessian of f imply that δ ≤ ϕ ′′ ≤ δ. We need the following lemma.

Lemma 5

Let h be a strictly convex function, i. e. h ′′ > 0 on a segment σ. Suppose that h attains its minimum at a point ζ 0 , and assume that h ′′ ≤ δ.

Then

h(ζ) -h(ζ 0 ) ≥ (h ′ (ζ)) 2 4 δ (25)
for all ζ ∈ σ.

Proof.

Suppose without loss of generality that ζ 0 < ζ, and let a point

ζ 1 be such that h ′ (ζ 1 ) = h ′ (ζ)
2 . This defines ζ 1 uniquely since the derivative h ′ of h is strictly increasing.

Then

h(ζ) -h(ζ 1 ) ≥ (ζ -ζ 1 ) • h ′ (ζ) 2 . ( 26
)
On the other hand, integrating

h ′′ from ζ 1 to ζ yields h ′ (ζ) 2 ≤ δ • (ζ -ζ 1 ). (27)
Multiplying ( 26) and ( 27) we obtain (25).

Lemma 5 gives

Lemma 6 In the notations introduced above (Notations 7)

|ϕ(t i ) -ϕ(τ 0 )| ≥ δ 2 η 2 4n δ .
This Lemma 6 describes what one gains (in terms of shortening the path in question) by setting t i = τ 0 . To use this lemma we introduce Notations 8 :

• gain= df |ϕ(t i ) -ϕ(τ 0 )| is the gain in the shortening of the path Π(T ) after setting t i = τ 0 , where τ 0 is the point where ϕ attain its minimum.

• U pBnd= df (gain • 4n δ) 1 2
δ is an upper bound on η (the latter denotes the distance between T * and T -see Notations 7).

Phase 2: Gradient Method Procedure.

Phase 2.1: Find i such that ∂f ∂t i (T ) = max j ∂f ∂t j (T ) . (See formulas (5) for

∂f ∂t i = ∂|Π| ∂t i )
. This is done by using formulas (5), which involve only arithmetic operations and square root extractions.

Claim 3 The complexity of Phase 2.1 is linear in n.

Phase 2.2: Find the point τ 0 where ϕ(τ ) attains its minimum. Since ϕ is strictly convex, this happens at the point where its derivative vanishes. Using our notations (see Notations 5 for ϕ and Notations 2 of section 2.3 for v i ) we have

ϕ(τ ) = v 0 (t 1 ) + v 1 (t 1 , t 2 ) + • • • + v i-1 (t i-1 , τ ) + v i (τ, t i+1 ) + • • • + v n (t n ).
Hence, using (5), we have

dϕ dτ (τ ) = dv i-1 (t i-1 , τ ) dτ + dv i (τ, t i+1 ) dτ = V i-1 (t i-1 , τ ), Λ i v i-1 - V i (τ, t i ), Λ i v i . (28)
Equating the expression (28) to zero, then squaring both sides and expressing the square of the length of a vector v j by the scalar product of the vector by itself, we get

V i-1 (t i-1 , τ ), Λ i 2 V i-1 (t i-1 , τ ), V i-1 (t i-1 , τ ) = V i (τ, t i ), Λ i 2 V i (τ, t i ), V i (τ, t i ) . ( 29
)
Notice that both V i-1 (t i-1 , τ ) and V i (τ, t i ) depend in τ linearly (see Notations 1). Thus, the equation ( 29) is of the 4th degree in τ . The classical method of L. Ferrari reduces this equation to equations of the 2nd and of the 3rd degree. This method uses a fixed number of arithmetical operations. The solutions of these equations of the 2nd and 3rd degree (formulas of G. Cardano for the latter) can be found using explicit formulas over arithmetic operations and square and cubic root extraction -all these operations are admissible in our computational model (see Introduction).

Hence,

Claim 4 The complexity of Phase 2.2 is O(1).

Note that, for a standard convexity reason, the modification of the path described in this phase cannot make it leave a ball where the path lied. Hence the path stays in the ball B (Notations 4).

Phase 2.3: Calculate (we use Notations 8): gain, the point T ′ = (t 1 , . . . , t i-1 , τ 0 , t i , . . . , t n ) and U pBnd, where τ 0 is from Phase 2.2 just above. 21) and ( 22)) then go to Phase 3 with T ′ as the initial approximation for Newton's method.

If ε √ n < U pBnd < b 2 (cf. (
Otherwise, repeat Phase 2 with T = T ′ . Thus, the recalculation of gain, T ′ and U pBnd is iterated while

U pBnd ≥ max{ ε √ n , b 2 }.
As one can see from the formulas defining the values of gain, T ′ and U pBnd (Notations 8, 5 and the respective notations from section 2.3),

Claim 5 The complexity of computing gain, T ′ and U pBnd is linear in n

Application of Newton's Method

Phase 3: Newton's method:

Apply Newton's method with the value of T , found by Phase 2.3, as X 0 . Iterate (13) until obtaining a sufficiently close approximation; see section 5 below, formula (46). The rate of convergence is estimated in Proposition 3, formula (14).

Complexity of the Algorithm

Now we summarize the comments on the complexity made in the previous sections to count the number of steps used by the algorithm to find a lengthposition ε-approximation of the shortest path touching the lines of L. By C i , c j we will denote various positive constants (using capital C for upper bounds, and lower case c for lower bounds).

Complexity of Preliminary computations

Phase 1 finds perpendiculars from s to lines L i . These perpendiculars fall at points P i ∈ L i . To find a point P i we solve a linear equation with one unknown t i : W i (t i )s , Λ i = 0. This can be done using only arithmetic operations.

The complexity of solving one equation is constant, thus the total complexity of finding all P i 's is linear in n. The points P i represent the polygon Π 0 , and computing its length |Π 0 | also involves square root extraction. However, the complexity is still linear in n.

Thus, the parameter R has been computed with linear complexity.

To compute sines α i = sin ∠Λ i Λ i+1 and cosines β i = cos ∠Λ i Λ i+1 we compute the scalar products Λ i , Λ i+1 , and use the standard formulas that involve arithmetic operations and square root extractions only. The complexity is again linear in n.

The next parameter to find is the minimal distance between consecutive lines. For two consecutive lines L i and L i+1 , the distance between them is the length of the segment connecting the lines and perpendicular to both of them; this is again a standard "analytic geometry" computation, which yields d and d in linear time. The calculation of δ and δ add only to the complexity (see formulas in Notations 5). Hence, Claim 6 The complexity of Phase 1 is linear in n.

Complexity of Gradient Descent

The complexity of one iteration of the gradient descent of Phase 2 was estimated from above at the end of section 4.3. This complexity is O(n).

So we are to estimate a sufficient number of iterations.

Recall that T -T * = η ≤ U pBnd (Notations 5 and Lemma 6). Phase 2.3 iterates the calculation of U pBnd until the latter becomes less than max{ ε

√ n , b 2 }.
It is evident that initially we have gain ≤ r (r is the length of the initial approximation to the shortest path -see Notations 4). Hence, initially for the initial gain (we use the expression for U pBnd from Notations 8)

U pBnd ≤ (gain • 4n δ) 1 2 δ ≤ c 3 (rn • n 1 2 ) 1 2 r 3 d 1 2 dα = c 3 • r 7 2 • n 3 4 d 5 2 • α = df B 0 . (30)
(We replaced d by d 2 and δ, δ by their expressions using respectively Notations 2 and Notations 7.)

So we can estimate the number of iterations in the Phase 2.3 as the ratio of the right-hand side value in formula (30) over a lower bound for gain that is valid while U pBnd goes down to the demanded value. Such a lower bound can be found from the condition U pBnd ≥ max ε √ n , b 2 that controls the continuation of the iterations:

U pBnd = (gain • 4n δ) 1 2 δ ≥ max ε √ n , b 2 (31)
This condition (31) gives

gain ≥ δ 2 4n δ • max ε 2 n , b 2 4 = (dα) 2 • d r 6 • 4n • 4 √ n • max ε 2 n , b 2 4 = (32) 5 • α 2 r 6 • 16n 3 2 • max ε 2 n , b 2 4 =                          d 5 α 2 r 6 • 16n 3 2 • b 2 4 = df g 1 , if ε 2 n < b 2 4 d 5 α 2 r 6 • 16n 3 2 • ε 2 n = df g 2 , if ε 2 n ≥ b 2 4 (33)
To estimate the number of iterations of recalculations of U pBnd in Phase 2 we estimate B 0 g 1 and B 0 g 2 . Denote the number of iterations by N bIter, and consider 2 cases corresponding to the cases in (33).

g 1 = d 5 α 2 r 6 • 16n 3 2 • c 2 αd 2 nr 3 = c ′ 2 • d 9 • α 3 r 9 • n 5 2 (34)
Divide B 0 from (30) by g 1 from (34)

N bIter ≤ O   r 7 2 • n 3 4 • r 9 • n 5 2 d 5 2 • α • d 9 • α 3   = O r 12.5 • n 3.25 d 11.5 • α 4 . ( 35
)
We take into account that d ≤ 1 < r (Notations 2, 4), α ≤ 1 and α = α 2 (Notations 3) and n ≥ 1, and rewrite (35) as

N bIter ≤ O r d 12.5 • n α 4 ≤ O rn d α 13 . (36)
Case 2:

ε 2 n ≥ b 2 4
. Divide B 0 from (30), this time, by g 2 from (33):

N bIter ≤ B 0 g 2 = O   r 7 2 • n 3 4 • r 6 • n 3 2 • n d 5 2 • α • d 5 • α 2 • ε 2   = O   r 19 2 • n 13 4 d 15 2 • α 3 • ε 2   (37)
In Case 2 we can bound 1 ε 2 using Case 2 condition:

1 ε 2 ≤ 4 b 2 • n = O r 6 • n d 4 • α 2 = O r 6 • n d 8 • α 2 , ( 38
)
where we used also (21)) to replace b and Notations 2 to replace d by d 2 . From (38) and (37) we get

N bIter ≤ O   r 19 2 • n 13 4 • r 6 • n d 15 2 • α 3 • d 8 • α 2   = O   r 31 2 • n 17 4 d 31 2 • α 5   = O r 15.5 • n 4.25 d 15.5 • α 5 ≤ (39) O r d 15.5 • n α 5 ≤ O rn d α 16 . (40)
The bound in (40) majorizes the bound in (36), so we can take the bound from (40) for further references. Hence, in any case .

N bIter ≤ O rn d α 16 . (41)

Complexity of Newton's Method

Once ε √ n ≤ b 2 , and this is the case we are mainly interested in, Phase 3 of our algorithm applies Newton's method. Formula (14) for the convergence of Combining complexity estimates for the three phases stated in Claims 6, 8 and 9 and observing that the term O n 2 (log r + log 1 α ) is much smaller than the bound on the complexity of the gradient procedure, we finally obtain the estimation of the Main Theorem.

Thus, the Main Theorem is proved.

Part 2: A Simple Grid Approximation of Shortest Paths Amidst Separated Obstacles 6 Separated Obstacles

We consider obstacles in R 3 though all the notions and constructions can be easily generalized to an arbitrary dimension. The dimension under consideration influences only the constants that appear in estimations of the complexity of algorithms.

The problem we study here is, as in the previous part, to construct a path that connects two given points s and t , does not intersect given obstacles W and whose length is close to that of a shortest path.

Obstacles and Paths

Let W be an arbitrary set in R 3 , and s and t be two points in its compliment.

Denote by cl(S), int(S) and bnd(S) respectively the closure, the interior and the boundary of a set S.

Denote by B v (a), where a ∈ R ≥0 and v ∈ R 3 , the ball of radius a centered at v.

The boundary bnd(W ) = cl(W ) \ int(W ) may contain 'degenerate' pieces. For example, an isolated point or a point with a neighborhood homeomorphic to a segment. Such pieces can hardly be considered as obstacles. So we assume that every point of bnd(W ) has a 2-dimensional neighborhood. For technical reasons it is convenient to assume a neighborhood of each point of bnd(W ) intersects int(W ). To achieve this we can 'slightly inflate' W . It is known how to do it efficiently for semi-algebraic obstacles, see e. g. [START_REF] Grigoriev | Polytime algorithm for the shortest path in a homotopy class amidst semi-algebraic obstacles in the plane[END_REF].

A path is a continuous piecewise smooth image of a closed segment. A simple path or a quasi-segment is a path without self-intersections (which is, clearly, homeomorphic to a segment).

The set R 3 \ W will be called the free space, and its closure will be called the space admissible for trajectories.

We consider only paths lying in the admissible space and not intersecting the interior of W .

Separability and Random Separated Balls

We say that obstacles W are a-separated if for any v ∈ R 3 the set W ∩ B v (a) is connected.

For example, if each connected component of W is convex and the distance between each two connected components is greater than 2a, then W is aseparated. However, the convexity of the obstacles is not assumed in the general case. What is imposed by a-separation is a certain smoothness of concave (from the point of view of an observer outside the obstacles) pieces of the boundary -a ball of radius 2a that goes inside 'fjords' of the an obstacle cannot touch or intersect two pieces of the obstacles that are 'remote' if to follow the boundary.

Our main goal is to describe an algorithm that constructs a length approximation to a shortest path under the condition of separability of obstacles. But now we will make a digression estimating the expectation of separability of obstacles constituted by n randomly chosen balls. The centers of the balls are independently chosen in the unit cube [0, 1] 3 , and their radii are independently chosen from an interval [0, r] under the uniform distribution.

Proposition 7 There exist constants c 1 , c 2 > 0 such that for any r < c 1 n -2 3 the union of n randomly chosen balls in [0, 1] 3 is (c 2 n -2 3 )-separated with probability greater than 2 3 .

Proof.

We claim that n randomly chosen points in the unit cube are c 3 n 2 3

-separated with probability greater than 2 3 for some constant c 3 > 0.

Indeed, a choice of n points in [0, 1] 3 is equivalent to a choice of a point p = (x 1 , y 1 , z 1 , . . . , x n , y n , z n ) ∈ [0, 1] 3n . For any pair of natural numbers i, j such that 1 ≤ i < j ≤ n, the measure of of those points p for which

|x i -x j |, |y i -y j |, |z i -z j | ≤ c 3 n -2 3 is not greater than √ 2c 3 n -2 3
3 . This is true for any constant c 3 . Hence the measure of those p for which there exist a pair 1 ≤ i < j ≤ n such that this condition holds is not greater than n(n-1) 2 √ 2c 3 n -2 3 3 ≤ 1 3 for an appropriate constant c 3 .

For finish the proof of the lemma we choose c 1 and c 2 so that 2c 1 + c 2 < c 3 .

A Grid Algorithm for a Shortest Path Approximation

In this section we assume that the obstacles W ⊆ [0, 1] 3 are a-separated and that s , t ∈ [0, 1] 3 .

Approximation Algorithm and its Complexity

Denote by L ⊆ [0, 1] 3 a cubic grid with mesh (edge of the basic cube) a.

Without loss of generality we assume that 1 a is an integer and that s and t are nodes of the grid.

Consider a graph G whose vertices are those nodes of the grid that do not belong to W and whose edges are those edges of the grid that do not intersect W .

Assume that the length of every edge in G is a.

The length of a path P will be denoted by |P |.

Theorem on Simple Grid Method. There is an algorithm of time complexity O( 1 a 6 ) that, using the graph G defined above, constructs a path Π G connecting s and t and whose length satisfies the inequality |Π G | ≤ 392|Π * | + 448a, where Π * is a shortest path amidst the obstacles W connecting s and t .

Remark 2. If a is not known but it is possible to construct the part of the grid admissible for paths, an algorithm can consecutively try a = 1 2 , 1 4 , 1 8 , . . . until it obtains a path with the bound from Theorem on Simple Grid Method.

From Proposition 7 one concludes

Corollary 8 When the set of obstacles W is a union of n random balls whose radii are independently chosen from [0, c 1 n -2 3 ], then a ≤ c 2 n -2 3 with probability greater than 2 3 . Hence the algorithm will compute an approximation (as described in Theorem on Simple Grid Method) to a shortest path with probability greater than 2 3 ) within time O(n 4 ).

We say that a path is a grid polygon if it goes only along the edges of the grid.

To prove the Theorem, it suffices to show that the length of a shortest grid polygon Π G connecting s and t and satisfies

|Π G | ≤ 84|Π * | + 96a. ( 49
)
Then this grid polygon can be constructed in quadratic time as a shortest path Π G in the graph G by using any polytime algorithm for a shortest path in a graph (e. g. see [START_REF] Aho | The Design and Analysis of Computer Algorithms[END_REF]).

Estimate (49), and hence Theorem on Simple Grid Method will follow from Lemmas 9 and 10.

We say that a cube of the grid is visited by a path P if the closure of this cube without vertices intersects with P . Notice that P does not determine uniquely the order of visited cubes as P may visit two cubes simultaneously by going along their common edge.

Lemma 9 If two nodes v 1 and v 2 of the grid are connected by a path P that visits s distinct cubes of the grid, then one can connect v 1 and v 2 by a grid polygon whose length is not greater than ≤ 12sa.

Proof.

For each cube K, we can consider a maximal segment of P contained in K.

Thus P is subdivided into intervals, each contained in one cube of the mesh and such that its continuation in either direction leaves the cube. Consider such a segment (a subpath) [wu] for a cube K, that is P enters K via a point w and leaves it via a point u (of course, P can make several visits like that to a cube). Let u lie in a face with vertices u 1 , u 2 , u 3 , u 4 . Then one of the four intervals uu i , 1 ≤ i ≤ 4, does not intersect the obstacles W .

Indeed, assume the contrary. Then, because of a-separability of W , the intersection W with any cube of the grid and, thus, with any of its faces is a convex set. The fact that a segment lying in the face intersects W means that there are points of the segment in the interior of W . But u ∈ int(W ). So if u is different from any u i , then our assumption implies that u ∈ int(W ), for u is the convex hull of the intersection of the face with int(W ). If u = u i for some i then uu i = u have the desired property.

Mark the points when P visits for the first time the 1st, the 9th, the 16th, ..., the (7l + 2)th, ... cube (we count only the visits to new cubes). The pieces of P between two consecutive marked points will be called intervals. Clearly, the number of intervals is at least s-2 7 .

It remains to show that the length of each interval is at least a. Reasoning by contradiction, note that the projection of an interval whose length is less than a to each coordinate axis is a segment of a length less than a. Hence it is contained in the interior of the union of two adjacent intervals of the form [ka, (k + 1)a]. Now the product of these intervals, which is the interior of the union of eight grid cubes with a common vertex, contains the interval. This contradicts the assumption that the interval visited at least 9 cubes.

Recall that s denotes the number of different cubes of the grid visited by a shortest path Π * . Lemma 10 gives |Π * | ≥ s-2 7 a ≥ s-2 7 -6 7 a. On the other hand, for the length of a shortest grid polygon Π G Lemma 9 gives |Π G | ≤ 12sa. Estimate (49) immediately follows from these two inequalities, and this concludes the proof of Theorem on Simple Grid Method.

Remark 3. The estimation given by Theorem on Simple Grid Method is, obviously, far from the exact one. We conjecture that |Π G | ≤ 7|Π * |. The latter bound cannot be essentially improved as one can see from the following example. Denote by v 1 , v 2 , v 3 , v 4 the vertices of the bottom face of the unit cube (ordered counter clockwise), and by u 1 , u 2 , u 3 , u 4 the respective vertices of the top face, see Figure 4. Let W 0 be the polyhedron with the following 5 vertices: 3 vertices at the middles of edges v 1 v 2 , v 3 v 4 , u 3 u 4 , and 2 vertices on edges u 1 v 1 and u 2 v 2 close to points v 1 and v 2 respectively. To obtain W , apply to W 0 a homothety (scaling) with coefficient (1 + ε), for a small enough ε > 0 centered at the center of the cube.

Fig. 1 .

 1 Fig. 1. Comparing triangles in M and on the plane.

Fig. 2 .

 2 Fig. 2. Hessian (of path length) ∆.

  and D = D i . For i = 0 we put x = 0, and for i = n we put y = 0. And let v = W, W . Within these notations d i = D, D , D ⊥ A and D ⊥ B.

Fig. 3 .

 3 Fig. 3. Comparing the lengths of two paths. and thus, |AB| -|A * B * | ≤ |AA * | + |BB * |. Similarly, we get |A * B * | -|AB| ≤ |AA * | + |BB * |. Hence, ||AB| -|A * B * || ≤ |AA * | + |BB * |.

Fig. 4 .

 4 Fig. 4. Some faces of W 0 .

Partially supported by an Alferd P. Sloan Fellowship and NSF Grant DMS-9803129.

2 Member of St-Petersburg Steklov Mathematical Institute, Russian Academy of Sciences, St-Petersburg, Russia

.3 Member of St-Petersburg Institute for Informatics, Russian Academy of Sciences, St-Petersburg, Russia.

2 |Π * |],where |Π * | is the length of Π * (estimate (3) in section 2). Observe that when α = 0 or d = 0 it could happen that the gradient descent and Newton's method which are invoked in our algorithm, do not converge. (There are considerations how to treat these particular cases but they need another, geometrically more 'invariant' approach.)

Acknowledgements. We are thankful to Misha Gromov for sharing his geometric insight, and to Yuri Burago for his useful comments and his help with the references. We are also grateful to anonymous referee for many useful remarks that helped to improve the presentation.

Algorithm for the Shortest Path Touching Lines

The algorithm takes as input a list L of lines, points s and t and ε, 1 ≥ ε > 0.

The algorithm outputs a path, represented as the list of points where it meets the lines of L. The output path is in the ε-neighborhood of the shortest path and its length is ε-close to the length of the shortest path.

We assume that the lines of L are represented by vectors Λ 0 i (point on the line L i ) and Λ i (unit vector directed along the line L i ) introduced in section 2.3. The complexity of transforming another usual representation of lines and points into this form is of linear complexity for our computation model.

The algorithm consists of three phases:

(Phase 1). Preliminary computations.

(Phase 2). Application of a gradient method to find an approximation for Newton's method.

(Phase 3). Application of Newton's method.

Below in this section we will describe the Phases. Their descriptions are rather short within the technique developed before, however some new notions will be needed. Together with this description we will make some estimations of complexity that concerns 'local' computations. The global estimation of the complexity will be done in the next section 5.

Length Approximation from a Position Approximation

The algorithm seeks an approximation of the position of the shortest path. An appropriate approximation of the length in our case is 'automatic' as follows from Lemma 4 below.

Proof.

The last inequality is standard. Let us prove the first one. We compare Π = Π(T ) and Π * = Π * (T * ) linkwise, see Figure 3, where

Newton's method estimates the error after k iterations by

which we want to be smaller than ε √ n . For our choice of initial approximation for Newton's method, see (23), and for our value of 1 c , see (20), the value (42) takes the form

Thus we need 1 2

Inequality ( 44) is equivalent to

Where all log's are with base 2; recall also that (without loss of generality) we assumed n ≤ r. From (45) we can bound k as

It remains to estimate the complexity of computing the kth iterate X k+1 using formula (13). In our case F (T ) = gradf (T ) and F ′ (T ) -1 = ∆(T ) -1 . The complexity of calculating F (T ) is shown to be linear in n. Computing the inverse of a 3-diagonal n × n-matrix ∆(T ) takes O(n 2 ) operations. Note that the complexity of computing one element of ∆(T ) is constant. Thus the total complexity of Newton's method is

Taking into consideration that log(x + y) ≤ (log x + log y) for x, y ≥ 2 , we can simplify this estimation (47), obtaining the following bound for the complexity of Newton's method Claim 9 The complexity of Phase 3 is

Now replace P by a path P ′ obtained by a sequence of the following modifications. Take any cube K visited by P , together with a maximum segment [wu] of P in K. Let i, 1 ≤ i ≤ 4, be such that the segment uu i does not intersect W . Insert the intervals uu i and u i u into P to force P to visit a vertex of the cube by going there and back. Repeating the same procedure for the entry point w (which is in its turn an exit point for some other cube), we modified our segment so that it ends and begins at a vertex of K. We will say that the visit of P to K terminates at u i . Perform this operation for all cubes visited by P . This gives P ′ , a new path which visits the same cubes, and for which every maximum subpath contained in a cube begins at and leaves the cube via a vertex.

We remark that though the number of cubes visited by P ′ is the same as for P , i. e. s, the length of P ′ might have increased.

Delete from P ′ all loops connecting the same vertex, thus obtaining a new path P ′′ . Notice that the number of vertices visited by P ′′ is at most 8s, for any vertex of any of the visited cubes is visited at most once. Moreover, P ′′ , still connects v 1 and v 2 . Now consider a maximum subpath of P ′′ contained in a grid cube K and connecting two vertices of K. Denote by W 1 (convex) intersection W ∩ K.

Consider the following homotopy of this subpath in K: pull each of its points v in the direction form the point of W 1 nearest to v. (The nearest point is unique since W 1 is convex.) Let the magnitude of the velocity of v (with respect to the parameter of the homotopy) be equal to the distance from v to the boundary ∂K of K. This way we "push" the maximum subpath of P ′′ in question away from W 1 until we transform it into a path lying in the boundary ∂K of K, connecting the same vertices of K and not intersecting W . Of course, its length still might have increased. Repeat the same procedure for each maximal subpath contained in one cube, having constructed a path going along faces.

Having repeated again the same operation for each maximum subpath contained in one boundary square of a grid cube (and now using a homotopy pushing this subpath to the boundary of this square), we end up with a polygonal path P ′′′ that consists of edges of the mesh, connects the same vertices v 1 and v 2 , and is contained in at most s cubes. After removing all loops from this path, we obtain a path P ′′′′ that traverses each edge at most once. Since s cubes together have at most 12s edges, we have constructed a path with required properties and whose length is at most 12as.

Lemma 10 If a path P connecting nodes v 1 and v 2 of the grid has visited s different cubes of the grid, then |P | ≥ s-2 7 a.

Proof.

For this obstacle the only grid polygon connecting v 1 and v 2 visits consecutively v 1 , v 4 , u 4 , u 1 , u 2 , u 3 , v 3 , v 2 . On the other hand, one can see that v 1 and v 2 can be connected by a path (avoiding W and contained in the cube) of length close to 1.