Approximating shortest path

for the skew lines problem

in time doubly logarithmic in 1/epsilon

D. Burago^{1}^{♯}^{∗}

Department of Mathematics, Penn State University, University Park, USA

D. Grigoriev^{2}^{♭}

Institut Mathématique de Rennes, CNRS Rennes, France

A. Slissenko^{3}^{∗}

Laboratoire d’Algorithmique, Complexité et Logique, University Paris-12, France

Version of September 10, 2003

Abstract

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 ℝ^{3}: given a list of straight lines L = (L_{1}, …, L_{n}) in ℝ^{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 lengthposition ε-approximated (i. e. both its length and its position can be found approximately) in time

( \frac{Rn}{d̃ ; α̃} )^{16} + O (n^{2} log log \frac{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 \frac{1}{ε}, while our algorithm for the problem under consideration has complexity linear in log log \frac{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 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 O (( \frac{1}{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.

1Address: Dept. of Comput. Sci., Penn State University, University Park, PA 16802, USA.

E-mail: burago@math.psu.edu

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

2Address: IMR, Université Rennes 1, Beaulieu 35042, Rennes, France.

E-mail: dima@math.univ-rennes1.fr

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

3Address: LACL, Dept. of Informatics, University Paris 12, 61 Av. du Gén. de Gaulle, 94010, Créteil, France.

E-mail: slissenko@univ-paris12.fr

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

1 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 ℝ^{3} in a specified order: given a list of straight lines L = (L_{1}, …, L_{n}) in 3-dimensional 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. [Pap85]). 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.

1.1 Related Results

The general problem of constructing a shortest path or even a ‘sufficiently good approximation’ for such a path is well known [Mit00]. It was shown to be NP-hard even for convex polyhedral obstacles [CR87] (by “polyhedral obstacles” we mean unions of polyhedra). The skew lines problem in ℝ^{3} is also well known. This problem is mentioned in [Pap85] as, presumably, representing the basic difficulties in constructing shortest paths in 3-dimensional Euclidean space. The paper [Pap85] 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 [SS90] 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 [CR87] 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 [CR87] it is proven that for polyhedral obstacles with the size of description N, determining O(\sqrt{N}) bits of the length of a shortest path is also NP-hard.

Concerning weaker approximations, for the case of polyhedral obstacles, paper [Pap85] 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 \frac{1}{ε}. Considerable gaps in this paper have been fixed in [CSY94] and the complexity of the algorithm [CSY94] is roughly O(\frac{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(\frac{n^{2}}{ε^{4}} ). An algorithm with complexity bound O(\frac{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(\frac{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 [SS86] in case of a convex surface, thereupon it was generalized in [MMP87] to arbitrary polyhedral surfaces, and its complexity is O(n^{2} log n). Finally the latter algorithm was improved in [CH90] providing a complexity O(n^{2}). To improve the quadratic complexity, several algorithms which give ε-approximations of a shortest path were produced: in [AHSV97] with complexity bound O(n + \frac{1}{ε^{3}} ) in case of a convex surface, then in [ALMS98] with complexity roughly O(\frac{n}{ε}) in case of a simple polyhedral surface. In [H99] for a fixed s , a preprocessed data structure is constructed with complexity O(\frac{n}{ε^{3}} ) in case of a convex surface and with complexity O(n^{2} + \frac{n}{ε}) in case of an arbitrary polyhedral surface, respectively, which allows one for any target t to output an approximation within complexity O(log(\frac{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 [AHU76]. Since the number of the placed Steiner points depends polynomially on n and on \frac{1}{ε}, the complexity of this approach should depend on n and on \frac{1}{ε} as well.

The approach of [MP91] to solve the weighted region problem gives an algorithm whose running time depends polynomially on log \frac{1}{ε}. Combining the binary search with the local optimization as in [MP91] one may probably find an approximation for the skew lines problem; however, in any case, our algorithm is exponentially faster in \frac{1}{ε} (and our text is shorter than that of [MP91]).

1.2 On Computational Models

In the present paper we also study the problem of approximating shortest paths in ℝ^{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. [HKSS94, GS99, GS98]. Note that a shortest path in ℝ^{3} between two rational points, with even one cylindrical obstacle, may have transcendental length and necessitate using transcendental points in its construction [HKSS93]. 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 ℝ^{3} our computational model is real RAMs [BSS89] 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.

1.3 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 (\frac{Rn}{d̃ ; α̃})^{16} + O(n^{2} log log \frac{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 ∈ [\sqrt{n}|Π^{∗}|, 2n^{\frac{3}{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.)

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 [Bal95, BN93] 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 [Bal95]). 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 [BFK98], [Bur98]. 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(\frac{1}{ε}) rather than polynomially on \frac{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 \frac{1}{ε} (logarithmic, cf. [MP91], 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 ℝ^{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 ℝ^{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 [Mit00]) that if the order of touching lines is not prescribed the problem of finding a shortest path becomes NP -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 [Pap77]. 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 [CSY94]. The matter is that the result of [Pap85, CSY94] is based on such a method. The initial idea used in [Pap85, CSY94] 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 [Pap85]. 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 ((\frac{1}{a})^{6}) and that for given a-separated 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 ℝ^{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 [ALMS98]: ”...while studying the performance of geometric algorithms, geometric parameters (e.g. fatness, density, aspect ratio, longest, closest) should not be ignored...”.

1.4 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 1). 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 3) 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 ℝ^{3}).

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

2.1 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\frac{1}{n})φ((i + 1)\frac{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 ℝ^{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.

Then we choose any point A′ in ℝ^{2} and take in ℝ^{2} two points B′_{t} and C′_{τ} that satisfy the equalities:

|A′B′_{t}|_{ℝ^{2}} = |Π_{AB}([0, t])|_{M}, |A′C′_{τ} |_{ℝ^{2}} = |Π_{AC}([0, τ])|_{M},

|B′_{t}C′_{τ} |_{ℝ^{2}} = |Π_{AB}(t)Π_{AC}(τ)|_{M}.

We are interested in the triangle △A′B′_{t}C′_{τ} in the plane ℝ^{2}.

Figure 1: Comparing triangles in M and on the plane.

Denote by θ(t, τ) the angle between the intervals A′B′_{t} and A′C′_{τ} in ℝ^{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 ℝ^{2} a triangle A′B′C′ determined by the following conditions:

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

|BC|_{M} for any three points A, B and C of M.

It is known that (see, for instance, [Bal95])

(Uniqueness of the Shortest Path.) In any simply connected space of nonpositive curvature, there is a unique shortest path between any two points.

2.2 The Configuration Space

Let L = (L_{1}, …, L_{n}) be a list of straight lines in ℝ^{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 ℝ^{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 [BN93]) 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 ℝ^{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 2, 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 \sqrt{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}).

2.3 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 𝕃_{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 ℝ^{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.

• Λ_{i}^{0} 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 ℝ^{n} of such points T.

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

• W_{i}(t_{i})=_{df} (Λ_{i}^{0} + 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 ℝ^{3} we denote by σ^{−} and σ^{+} its beginning and its end respectively.

Notations 2 (Lengths):

• ρ will be used to denote various 𝕃_{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} + \sqrt{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 ℝ^{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} + \sqrt{α}.

• β_{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 ℝ^{n} defining Π^{∗}, i. e. Π^{∗} = Π(T^{∗}).

• Π^{0} is the initial approximation of Π^{∗} defined below in section 2.4.

• T^{0} is the point in ℝ^{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\sqrt{n} is the radius of a ball in ℝ^{n} that will contain all paths in question.

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

2.4 Initial Approximation for the Shortest Path

Denote by P_{i}^{0} 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_{1}^{0}, then to P_{2}^{0} and so on to P_{n}^{0} 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 ∑_{i=1}^{n−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

\sqrt{n}|Π^{∗}| ≤ R ≤ 2n^{\frac{3}{2}} |Π^{∗}|. (3)

We remark that if L_{i} = {W_{i}^{0} + tΛ_{i} : t ∈ ℝ}, then P_{i} = W_{i}^{0} + t_{i}^{0}Λ_{i} is determined by ⟨P_{i}, Λ_{i}⟩ = 0, and thus,

P_{i} = W_{i}^{0} − ⟨W_{i}^{0}, Λ_{i}⟩ Λ_{i}. (4)

The parameter R = |Π^{0}|\sqrt{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_{i}^{0} 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.

2.5 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

(\frac{Rn}{d̃ ; α̃})^{16} + O (n^{2} log log \frac{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 (\frac{d̃ ^{4} α ^{2}}{r^{3} \sqrt{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.

3 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} |Π| = ∑_{i=0}^{n} v_{i} is the length of Π(T), where v_{i} is defined in Notations 2.

• g(T)=_{df} f′(T)=_{df} gradf(T) = (\frac{∂f}{∂t_{1}} (T), …, \frac{∂f}{∂t_{n}} (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} = \frac{∂f}{∂t_{i}} = \frac{∂}{∂t_{i}} ( v_{i−1} + v_{i} ) = \frac{∂v_{i−1}}{∂t_{i}} + \frac{∂v_{i}}{∂t_{i}} = \frac{⟨V_{i−1}, Λ_{i}⟩}{v_{i−1}} − \frac{⟨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.

[ \begin{matrix}\frac{∂^{2}v_{0}}{∂t_{1}^{2}} + \frac{∂^{2}v_{1}}{∂t_{1}^{2}} & \frac{∂^{2}v_{1}}{∂t_{1}∂t_{2}} & 0 & … & 0 & 0 \\ \frac{∂^{2}v_{1}}{∂t_{1}∂t_{2}} & \frac{∂^{2}v_{1}}{∂t_{2}^{2}} + \frac{∂^{2}v_{2}}{∂t_{2}^{2}} & \frac{∂^{2}v_{2}}{∂t_{2}∂t_{3}} & … & 0 & 0 \\ 0 & \frac{∂^{2}v_{2}}{∂t_{2}∂t_{3}} & \frac{∂^{2}v_{2}}{∂t_{3}^{2}} + \frac{∂^{2}v_{3}}{∂t_{3}^{2}} & … & 0 & 0 \\ . . . . . . & . . . . . . & . . . . . . & … & . . . . . . & . . . . . . \\ 0 & 0 & 0 & … & \frac{∂^{2}v_{n−1}}{∂t_{n−1}∂t_{n}} & \frac{∂^{2}v_{n−1}}{∂t_{n}^{2}} + \frac{∂^{2}v_{n}}{∂t_{n}^{2}} \end{matrix} ]

Figure 2: Hessian (of path length) Δ.

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

\frac{∂^{2}v_{i}}{∂t_{i}∂t_{i+1}} = \frac{∂}{∂t_{i}} (\frac{⟨V_{i}, Λ_{i+1}⟩}{v_{i}}) = \frac{⟨\frac{∂V_{i}}{∂t_{i}} , Λ_{i+1}⟩}{v_{i}} − \frac{⟨V_{i}, Λ_{i+1}⟩ ⟨V_{i}, \frac{∂V_{i}}{∂t_{i}} ⟩}{v_{i}^{3}}

= −\frac{⟨Λ_{i}, Λ_{i+1}⟩}{v_{i}} + \frac{⟨V_{i}, Λ_{i+1}⟩ ⟨V_{i}, Λ_{i}⟩}{v_{i}^{3}} = (7)

= −\frac{1}{v_{i}} ( cos ∠Λ_{i}Λ_{i+1} − cos ∠V_{i}Λ_{i+1} · cos ∠V_{i}Λ_{i} ) (8)

\frac{∂^{2}v_{i}}{∂t_{i}^{2}} = \frac{1}{v_{i}} (1 − \frac{⟨V_{i}, Λ_{i}⟩^{2}}{v_{i}^{2}} ) = \frac{(sin ∠V_{i}, Λ_{i})^{2}}{v_{i}} . (9)

\frac{∂^{2}v_{i}}{∂t_{i+1}^{2}} = \frac{1}{v_{i}} (1 − \frac{⟨V_{i}, Λ_{i+1}⟩^{2}}{v_{i}^{2}} ) = \frac{(sin ∠V_{i}, Λ_{i+1})^{2}}{v_{i}} . (10)

3.2 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} [\frac{∂^{2}v_{0}}{∂t_{1}^{2}} ]

• Δ_{i}=_{df} [ \begin{matrix}\frac{∂^{2}v_{i}}{∂t_{i}^{2}} & \frac{∂^{2}v_{i}}{∂t_{i}∂t_{i+1}} \\ \frac{∂^{2}v_{i}}{∂t_{i}∂t_{i+1}} & \frac{∂^{2}v_{i}}{∂t_{i+1}^{2}} \end{matrix} ] for 1 ≤ i ≤ n − 1

• Δ_{n}=_{df} [\frac{∂^{2}v_{n}}{∂t_{n}^{2}} ]

• Δ̃_{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 metric – it 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 \frac{d}{2r^{3}} , and both eigenvalues of Δ_{i} for 1 ≤ i ≤ n − 1 are greater than \frac{dα}{2r^{3}} . Thus \frac{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 ℝ^{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 W = V_{i}, x = t_{i}, A = Λ_{i}, y = t_{i+1}, B = Λ_{i+1} and D = D_{i}. For i = 0 we put x = 0, and for i = n we put y = 0. And let v = \sqrt{⟨W, W ⟩}.

Within these notations d_{i} = ⟨D, D⟩, D ⊥ A and D ⊥ B.

For i = 0 we have that the only eigenvalue of Δ_{0} is its element

\frac{∂^{2}v_{0}}{∂t_{1}^{2}} = \frac{1}{v_{0}} (1 − \frac{⟨W,B⟩^{2}}{v_{0}^{2}} ). Here ⟨W, B⟩ = ⟨yB + D, B⟩ = y, ⟨W, W ⟩ = v_{0}^{2} = ⟨yB + D, yB + D⟩ = y^{2} + d_{0} and thus

∥Δ_{0}∥ = |\frac{1}{v_{0}} (1 − \frac{y^{2}}{v_{0}^{2}} )| = \frac{1}{v_{0}^{3}} (y^{2} + d_{0} − y^{2}) = \frac{d_{0}}{v_{0}^{3}} > \frac{d}{r^{3}} .

Similarly for i = n we have ∥Δ_{n}∥ ≥ \frac{d_{n}}{v_{n}^{3}} > \frac{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} = [ \begin{matrix}\frac{∂^{2}v}{∂x^{2}} & \frac{∂^{2}v}{∂x∂y} \\ \frac{∂^{2}v}{∂x∂y} & \frac{∂^{2}v}{∂y^{2}} \end{matrix} ] .

Now we compute the second derivative of v:

\frac{∂^{2}v}{∂x^{2}} = \frac{1}{v} (1 − \frac{⟨xA + yB + D, A⟩^{2}}{v^{2}} ) = \frac{1}{v^{3}} (v^{2} − (x + yβ̃_{i})^{2})

= \frac{1}{v^{3}} (y^{2} + d_{i} − (yβ̃_{i})^{2}) = \frac{1}{v^{3}} (y^{2}α_{i} + d_{i}) .

\frac{∂^{2}v}{∂y^{2}} = \frac{1}{v} (1 − \frac{⟨xA + yB + D, B⟩^{2}}{v^{2}} ) = \frac{1}{v^{3}} (v^{2} − (xβ̃_{i} + y)^{2})

= \frac{1}{v^{3}} (x^{2}α_{i} + d_{i}) .

\frac{∂^{2}v}{∂x∂y} = −\frac{1}{v_{i}} (β̃_{i} − \frac{⟨xA + yB + D, B⟩ ⟨xA + yB + D, A⟩}{v^{2}})

= −\frac{1}{v^{3}} (β̃_{i}(x^{2} + 2xyβ̃_{i} + y^{2} + d_{i}) − (xβ̃_{i} + y)(x + yβ̃_{i}))

= −\frac{1}{v^{3}} ((x^{2}β̃_{i} + 2xyβ̃_{i}^{2} + y^{2}β̃_{i} + d_{i}β̃_{i}) − (x^{2}β̃_{i} + xy + xyβ̃_{i}^{2} + y^{2}β̃_{i}))

= −\frac{1}{v^{3}} (xyβ̃_{i}^{2} + d_{i}β̃_{i} − xy) = −\frac{1}{v^{3}} (d_{i}β̃_{i} − xyα_{i}) .

Let us estimate the determinant

detΔ_{i} = \frac{1}{v^{6}} ((y^{2}α_{i} + d_{i})(x^{2}α_{i} + d_{i}) − (d_{i}β̃_{i} − xyα_{i})^{2})

= \frac{1}{v^{6}} (x^{2}y^{2}α_{i}^{2} + d_{i}α_{i}(x^{2} + y^{2}) + d_{i}^{2} − d_{i}^{2}β̃_{i}^{2} − x^{2}y^{2}α_{i}^{2} + 2xyα_{i}d_{i}β̃_{i})

= \frac{1}{v^{6}} (d_{i}α_{i}(x^{2} + y^{2}) + d_{i}^{2}α_{i} + 2xyα_{i}d_{i}β̃_{i})

= \frac{α_{i}}{v^{6}} (d_{i}(x^{2} + y^{2}) + d_{i}^{2} + 2xyd_{i}β̃_{i}) ≥

\frac{α_{i}}{v^{6}} (d_{i}(x^{2} + y^{2}) + d_{i}^{2} − 2|x||y|d_{i}) ≥ \frac{d_{i}^{2}α_{i}}{v^{6}} .

Note that TrΔ_{i} = \frac{1}{v^{3}} (α_{i}(x^{2} + y^{2}) + 2d_{i}). The smallest of two eigenvalues of Δ_{i} is not less than

\frac{detΔ_{i}}{TrΔ_{i}} ≥ \frac{v^{3}d_{i}^{2}α_{i}}{v^{6}(α_{i}(x^{2} + y^{2}) + 2d_{i})} ≥ \frac{d_{i}^{2}α_{i}}{v^{3}2d_{i}} = \frac{d_{i}α_{i}}{2v^{3}} > \frac{dα}{2r^{3}} .

■

Proposition 1 All eigenvalues of Δ(T) are greater than \frac{dα}{2r^{3}} . Thus Δ(T) is positive definite, the metric function is strictly convex and ∥Δ(T)∥ ≥ \frac{dα}{r^{3}}, and ∥Δ^{−1}(T)∥ ≤ \frac{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} \frac{dα}{2r^{3}} ⟨U_{i}, U_{i}⟩ = \frac{dα}{2r^{3}} ∑_{i} ⟨U_{i}, U_{i}⟩ = \frac{dα}{2r^{3}} · 2 = \frac{dα}{r^{3}}.

This inequality gives a lower bound on eigenvalues of Δ. Hence, ∥Δ∥ ≥ \frac{dα}{r^{3}}, and ∥Δ^{−1}∥ ≤ \frac{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)∥ ≤ \sqrt{n} · ∥Δ∥_{1} = \sqrt{n} · max_{i} ∑_{j} |Δ_{i,j}|

≤ 4\sqrt{n} max_{i} \{\frac{1}{v_{i}} \} ≤ \frac{4\sqrt{n}}{d̃}. (12)

for all T.

3.3 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 ℝ^{n} into ℝ^{n}. Letter X without or with superscripts and subscripts will denote vector-columns from ℝ^{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 2 gives a sufficient condition to ensure a fast convergence of Newton’s iterations (it can be found in textbooks on numerical methods).

Proposition 2 Suppose that X_{0} is a zero of H, i. e. H(X_{0}) = \overline{0}, where \overline{0} = (0, …, 0) ∈ ℝ^{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, \frac{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}∥ ≤ \frac{1}{c} ( c∥X^{0} − X_{0}∥ )^{2^{k}} . (14)

In Proposition 2 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 = \frac{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}∥ ≤ \frac{1}{2}.

We will take the bound from Proposition 1 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 2. 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_{1}^{(j)}, …, t_{n}^{(j)}), j = 1, 2, we can write the ith component of (15) as

π_{i}(t_{i−1}^{(1)}, t_{i}^{(1)}, t_{i+1}^{(1)}) − π_{i}(t_{i−1}^{(2)}, t_{i}^{(2)}, t_{i+1}^{(2)}) + ∑_{j} \frac{∂π_{i}(t_{j}^{(2)})}{∂t_{j}}(t_{j}^{(1)} − t_{j}^{(2)}). (16)

Taylor’s formula says that (16) is equal to

\frac{1}{2} ( ∑_{j=i−1,i,i+1} \frac{∂^{2}π_{i}(ξ)}{∂t_{j}∂(}t′_{j} − t_{j})^{2} + ∑_{j,k=i−1,i,i+1} \frac{∂^{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_{j}^{(1)} and t_{j}^{(2)}.

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

\frac{1}{v_{i}^{2}} ∑(sin^{O(1)} φ, cos^{O(1)} ψ, (O(1))) + \frac{1}{v_{i−1}^{2}} ∑(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(\frac{1}{d}). Hence, the absolute value of each component of vector (15) is bounded by O(\frac{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} (\frac{\sqrt{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} · \frac{\sqrt{n}}{d}. (19)

As in our case H′ = Δ, Proposition 1 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} = \frac{r^{3}}{dα} (see Proposition 1), a_{2} = C_{1} (\frac{n}{d}) (see (18)),

• c = a_{1}a_{2} = C_{1} (\frac{nr^{3}}{αd^{2}} ) for the constant C_{1} from (18).

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

Hence we take

b = \frac{1}{c} = \frac{αd^{2}}{C_{1}nr^{3}} = \frac{α^{2}d̃^{4}}{C_{1}nr^{3}} . (20)

For further references we rewrite (20 as (we use Notations 2, 3)

• \frac{b}{2} = c_{2} \frac{α^{2}d̃^{4}}{nr^{3}} = c_{2} \frac{αd^{2}}{nr^{3}} , (21)

where c_{2}=_{df} \frac{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} ∈ Ω_{\frac{b}{2}} . (22)

For such X^{0} we have

c∥X^{0} − X_{0}∥ ≤ \frac{1}{2}, (23)

and the rate of convergence guaranteed by Proposition 2 will be (\frac{1}{2})^{2^{k}}, where k is the number of iterations (see Proposition 2).

4 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 Λ_{i}^{0} (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.

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

Lemma 2 |Π(T)| − |Π(T^{∗})| ≤ 2∥T − T^{∗}∥_{1} ≤ \sqrt{n}∥T − T^{∗}∥.

Proof. The last inequality is standard. Let us prove the first one. We compare Π = Π(T) and Π^{∗} = Π^{∗}(T^{∗}) linkwise, see Figure 3, where Π_{i} = AB, Π_{i}^{∗} = A^{∗}B^{∗}, |AA^{∗}| = |t_{i} − t_{i}^{∗}|, |BB^{∗}| = |t_{i+1} − t_{i+1}^{∗}|. The triangle inequality immediately implies that

Figure 3: Comparing the lengths of two paths.

|AB| − |A^{∗}B| ≤ |AA^{∗}|, |A^{∗}B| − |A^{∗}B^{∗}| ≤ |BB^{∗}| and thus, |AB| − |A^{∗}B^{∗}| ≤ |AA^{∗}| + |BB^{∗}|.

Similarly, we get |A^{∗}B^{∗}| − |AB| ≤ |AA^{∗}| + |BB^{∗}|. Hence, ||AB| − |A^{∗}B^{∗}|| ≤ |AA^{∗}| + |BB^{∗}|.

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

\frac{ε}{\sqrt{n}}-approximation.

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

4.3 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 1), and the point at which it attains its minimum is denoted by T^{∗} (section 2.3, Notations 4). Proposition 1 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 = (\frac{∂f}{∂t_{1}} , …, \frac{∂f}{∂t_{n}} ). 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 |\frac{∂f}{∂t_{i}} (T)| of one of the coordinates of the gradient g(T) is at least \frac{δ·η}{\sqrt{n}}.

Choose i such that

|\frac{∂f}{∂t_{i}} (T)| ≥ \frac{δ · η}{\sqrt{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 :

• δ̃ = \frac{4\sqrt{n}}{d̃} is an upper bound for Δ (formula (12)).

• δ = \frac{dα}{r^{3}} is a lower bound for Δ (Proposition 1).

• g=_{df} gradf = (\frac{∂f}{∂t_{1}} , …, \frac{∂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 3

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}) ≥ \frac{(h′(ζ))^{2}}{4δ̃} (25)

for all ζ ∈ σ.

Proof. Suppose without loss of generality that ζ_{0} < ζ, and let a point ζ_{1} be such that h′(ζ_{1}) = \frac{h′(ζ)}{2}. This defines ζ_{1} uniquely since the derivative h′ of h is strictly increasing.

Then

h(ζ) − h(ζ_{1}) ≥ (ζ − ζ_{1}) · \frac{h′(ζ)}{2}. (26)

On the other hand, integrating h′′ from ζ_{1} to ζ yields

\frac{h′(ζ)}{2} ≤ δ̃ · (ζ − ζ_{1}). (27)

Multiplying (26) and (27) we obtain (25). ■ Lemma 3 gives Lemma 4 In the notations introduced above (Notations 7)

|φ(t_{i}) − φ(τ_{0})| ≥ \frac{δ^{2}η^{2}}{4nδ̃} .

This Lemma 4 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.

• UpBnd=_{df} \frac{(gain · 4nδ̃)^{\frac{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 |\frac{∂f}{∂t_{i}} (T)| = max_{j} \{|\frac{∂f}{∂t_{j}} (T)|\}. (See formulas (5) for \frac{∂f}{∂t_{i}} = \frac{∂|Π|}{∂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

\frac{dφ}{dτ}(τ) = \frac{dv_{i−1}(t_{i−1}, τ)}{dτ} + \frac{dv_{i}(τ, t_{i+1})}{dτ} =

\frac{⟨V_{i−1}(t_{i−1}, τ), Λ_{i}⟩}{v_{i−1}} − \frac{⟨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

\begin{matrix}⟨V_{i−1}(t_{i−1}, τ), Λ_{i}⟩ & 2 \\ ⟨ & V_{i−1}(t_{i−1}, τ), V_{i−1}(t_{i−1}, τ) & ⟩ \end{matrix} = \begin{matrix}⟨V_{i}(τ, t_{i}), Λ_{i}⟩ & 2 \\ ⟨ & V_{i}(τ, t_{i}), V_{i}(τ, t_{i}) & ⟩ \end{matrix}. (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 UpBnd, where τ_{0} is from Phase 2.2 just above.

If \frac{ε}{\sqrt{n}} < UpBnd < \frac{b}{2} (cf. (21) and (22)) then go to Phase 3 with T′ as the initial approximation for Newton’s method.

Otherwise, repeat Phase 2 with T = T′. Thus, the recalculation of gain, T′ and UpBnd is iterated while UpBnd ≥ max\{\frac{ε}{\sqrt{n}}, \frac{b}{2}\}.

As one can see from the formulas defining the values of gain, T′ and UpBnd (Notations 8, 5 and the respective notations from section 2.3),

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

4.4 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 2, formula (14).

5 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 length-position ε-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).

5.1 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 O(1) to the complexity (see formulas in Notations 5). Hence,

Claim 6 The complexity of Phase 1 is linear in n.

5.2 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^{∗}∥ = η ≤ UpBnd (Notations 5 and Lemma 4). Phase 2.3 iterates the calculation of UpBnd until the latter becomes less than max\{\frac{ε}{\sqrt{n}}, \frac{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 UpBnd from Notations 8)

UpBnd ≤ \frac{(gain · 4nδ̃)^{\frac{1}{2}}}{δ} ≤ c_{3} \frac{(rn · n^{\frac{1}{2}} )^{\frac{1}{2}} r^{3}}{d̃^{\frac{1}{2}} dα} = c_{3} · \frac{r^{\frac{7}{2}} · n^{\frac{3}{4}}}{d̃^{\frac{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 UpBnd goes down to the demanded value. Such a lower bound can be found from the condition UpBnd ≥ max \{\frac{ε}{\sqrt{n}}, \frac{b}{2}\} that controls the continuation of the iterations:

UpBnd = \frac{(gain · 4nδ̃)^{\frac{1}{2}}}{δ} ≥ max \{\frac{ε}{\sqrt{n}}, \frac{b}{2}\} (31)

This condition (31) gives

gain ≥ \frac{δ^{2}}{4nδ̃} · max \{\frac{ε^{2}}{n} , \frac{b^{2}}{4} \} = \frac{(dα)^{2} · d̃}{r^{6} · 4n · 4\sqrt{n}} · max \{\frac{ε^{2}}{n} , \frac{b^{2}}{4} \} = (32)

\frac{d̃^{5} · α^{2}}{r^{6} · 16n^{\frac{3}{2}}} · max \{\frac{ε^{2}}{n} , \frac{b^{2}}{4} \} = \{ \begin{matrix}\frac{d̃^{5}α^{2}}{r^{6} · 16n^{\frac{3}{2}}} · \frac{b^{2}}{4} =_{df} g_{1}, & if \frac{ε^{2}}{n} < \frac{b^{2}}{4} \\ \frac{d̃^{5}α^{2}}{r^{6} · 16n^{\frac{3}{2}}} · \frac{ε^{2}}{n} =_{df} g_{2}, & if \frac{ε^{2}}{n} ≥ \frac{b^{2}}{4} \end{matrix} (33)

To estimate the number of iterations of recalculations of UpBnd in Phase 2 we estimate \frac{B_{0}}{g_{1}} and \frac{B_{0}}{g_{2}} . Denote the number of iterations by NbIter, and consider 2 cases corresponding to the cases in (33).

Case 1: \frac{ε^{2}}{n} < \frac{b^{2}}{4} . Replace \frac{b^{2}}{4} by its value from (21)

g_{1} = \frac{d̃^{5}α^{2}}{r^{6} · 16n^{\frac{3}{2}}} · c_{2} \frac{αd^{2}}{nr^{3}} = c′_{2} · \frac{d̃^{9} · α^{3}}{r^{9} · n^{\frac{5}{2}}} (34)

Divide B_{0} from (30) by g_{1} from (34)

NbIter ≤ O (\frac{r^{\frac{7}{2}} · n^{\frac{3}{4}} · r^{9} · n^{\frac{5}{2}}}{d̃^{\frac{5}{2}} · α · d̃^{9} · α^{3}} ) = O (\frac{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

NbIter ≤ O ((\frac{r}{d̃})^{12.5} · (\frac{n}{α})^{4}) ≤ O (\frac{rn}{d̃ ; α̃})^{13} . (36)

Case 2: \frac{ε^{2}}{n} ≥ \frac{b^{2}}{4} . Divide B_{0} from (30), this time, by g_{2} from (33):

NbIter ≤ \frac{B_{0}}{g_{2}} = O (\frac{r^{\frac{7}{2}} · n^{\frac{3}{4}} · r^{6} · n^{\frac{3}{2}} · n}{d̃^{\frac{5}{2}} · α · d̃^{5} · α^{2} · ε^{2}}) = O (\frac{r^{\frac{19}{2}} · n^{\frac{13}{4}}}{d̃^{\frac{15}{2}} · α^{3} · ε^{2}} ) (37)

In Case 2 we can bound \frac{1}{ε^{2}} using Case 2 condition:

\frac{1}{ε^{2}} ≤ \frac{4}{b^{2} · n} = O (\frac{r^{6} · n}{d^{4} · α^{2}} ) = O (\frac{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

NbIter ≤ O (\frac{r^{\frac{19}{2}} · n^{\frac{13}{4}} · r^{6} · n}{d̃^{\frac{15}{2}} · α^{3} · d̃^{8} · α^{2}} ) = O (\frac{r^{\frac{31}{2}} · n^{\frac{17}{4}}}{d̃^{\frac{31}{2}} · α^{5}} ) = O (\frac{r^{15.5} · n^{4.25}}{d̃^{15.5} · α^{5}} ) ≤ (39)

O ((\frac{r}{d̃})^{15.5} · (\frac{n}{α})^{5}) ≤ O (\frac{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

NbIter ≤ O (\frac{rn}{d̃ ; α̃})^{16} . (41)

From this bound (41) and Claim 5 that says that the complexity of calculating of gain, T′ and UpBnd in Phase 2.3 is linear we get

Claim 7 The complexity of Phase 2.3 is O (\frac{rn}{d̃ ; α̃})^{16} that is majorized by O (\frac{Rn}{d̃ ; α̃})^{16}.

From Claims 3, 4 and 7 we deduce

Claim 8 The complexity of Phase 2 is O (\frac{rn}{d̃ ; α̃})^{16} or, in other terms,

O (\frac{Rn}{d̃ ; α̃})^{16}.

5.3 Complexity of Newton’s Method

Once \frac{ε}{\sqrt{n}} ≤ \frac{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 Newton’s method estimates the error after k iterations by

\frac{1}{c} ( c∥X^{0} − X_{0}∥ )^{2^{k}} , (42)

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

\frac{1}{c} (\frac{1}{2})^{2^{k}} = b (\frac{1}{2})^{2^{k}} = \frac{αd^{2}}{C_{1}nr^{3}} (\frac{1}{2})^{2^{k}} . (43)

Thus we need

(\frac{1}{2})^{2^{k}} < \frac{C_{1}\sqrt{n}r^{3}ε}{αd^{2}}. (44)

Inequality (44) is equivalent to 2^{k} < log \frac{1}{ε} + \frac{1}{2} log n + 3 log r + log \frac{1}{α} − 2 log d + O(1) ≤

log \frac{1}{ε} + O (log r + log \frac{1}{α}) . (45)

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

k < log (log \frac{1}{ε} + O(log r + log \frac{1}{α})) . (46)

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

O (n^{2} log(log \frac{1}{ε} + O(log r + log \frac{1}{α}))) . (47)

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

O (n^{2} log log \frac{1}{ε}) + O (n^{2} (log r + log \frac{1}{α})) . (48)

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 \frac{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 ℝ^{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.

6.1 Obstacles and Paths

Let W be an arbitrary set in ℝ^{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 ∈ ℝ_{≥0} and v ∈ ℝ^{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. [GS98].

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 ℝ^{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.

6.2 Separability and Random Separated Balls

We say that obstacles W are a-separated if for any v ∈ ℝ^{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 a-separated. 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 3 There exist constants c_{1}, c_{2} > 0 such that for any r < c_{1}n^{− \frac{2}{3}} the union of n randomly chosen balls in [0, 1]^{3} is (c_{2}n^{− \frac{2}{3}} )-separated with probability greater than \frac{2}{3}.

Proof. We claim that n randomly chosen points in the unit cube are \frac{c_{3}}{n^{\frac{2}{3}}} -separated with probability greater than \frac{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^{− \frac{2}{3}} is not greater than (\sqrt{2}c_{3}n^{− \frac{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 \frac{n(n−1)}{2} (\sqrt{2}c_{3}n^{− \frac{2}{3}} )^{3} ≤ \frac{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}. ■

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

7.1 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 \frac{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((\frac{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 = \frac{1}{2}, \frac{1}{4}, \frac{1}{8}, … until it obtains a path with the bound from Theorem on Simple Grid Method.

From Proposition 3 one concludes

Corollary 1 When the set of obstacles W is a union of n random balls whose radii are independently chosen from [0, c_{1}n^{− \frac{2}{3}} ], then a ≤ c_{2}n^{− \frac{2}{3}} with probability greater than \frac{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 \frac{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 [AHU76]).

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

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

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} the (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 6 If a path P connecting nodes v_{1} and v_{2} of the grid has visited s different cubes of the grid, then |P | ≥ ⌊\frac{s−2}{7}⌋ a.

Proof. 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 ⌊\frac{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 6 gives |Π^{∗}| ≥ ⌊\frac{s−2}{7}⌋ a ≥ (\frac{s−2}{7} − \frac{6}{7}) a. On the other hand, for the length of a shortest grid polygon Π_{G} Lemma 5 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

Figure 4: Some faces of W_{0}.

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.

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.

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.

References

[AHSV97] P.K. Agarwal, S. Har-Peled, M. Sharir, and K.R. Varadarajan. Approximate shortest paths on a convex polytope in three dimensions. J.Assoc. Comput. Mach., 44:567–584, 1997.

[AHU76] A. Aho, J. Hopcroft, and J. Ullman. The Design and Analysis of Computer Algorithms. Addison-Wesley, 1976.

[ALMS98] L. Alexandrov, M. Lanthier, A. Maheshwari, and J.-R.Sack. An ϵ-approximation algorithm for weighted shortest paths on polyhedral surfaces. In Proc. 6th Scand. Workshop Algorithm Theory, vol. 1432 of Lecture Notes Comput. Sci., pp.11–22, Springer-Verlag, 1998.

[Bal95] W. Ballmann. Lectures on Spaces of Nonpositive Curvature. DMV Seminar, Band 25. Birkhäuser, Basel, Boston, Berlin, 1995.

[BN93] V. Berestovsij and I. Nikolaev. Multidimensional generalized Riemann spaces. In Yu. Reshetnyak, editor, Geometry 4. Non-regular Riemann Geometry, Encycloædia of Math. Sciences, vol. 70. Springer-Verlag, 1993.

[BSS89] L. Blum, M. Shub, and S. Smale. On a theory of computation and complexity over real numbers: NP-completeness, recursive functions and universal machines. Bull. Amer. Math. Soc., 1:1–46, 1989.

[Bur98] D.Burago. Hard Balls Gas and Alexandrov Spaces of Curvature Bounded Above. Proceedings of the ICM, Invited Lecture, vol 2, pages 289-298, 1998.

[BFK98] D. Burago, S. Ferleger, A. Kononenko Uniform estimates on the number of collisions in semidispersing billiards. Annals of Mathematics, 147, pages 695-708, 1998.

[CR87] J. Canny and J. Reif. New lower bound technique for robot motion planning problemms. In Proc. 28th Annu. IEEE Symp. on Foundations of Comput. Sci., pages 49–60, 1987.

[CH90] J. Chen, and Y. Han. Shortest paths on a polyhedron. In Proc. 6th Annu. ACM Symp. on Computational Geometry, (Berkeley, June 6-8), New York, pp.360–369, 1990.

[CSY94] J. Choi, J. Sellen, and C.-K. Yap. Approximate Euclidean shortest path in 3-space. In Proc. 10th ACM Symp. on Computational Geometry, pages 41–48, 1994.

[C87] K. L. Clarkson. Approximate algorithms for shortest path motion planning. In Proc. 19th Annu. ACM Symp. Th. Comput., pp. 56–65, 1987.

[GS98] D. Grigoriev and A. Slissenko. Polytime algorithm for the shortest path in a homotopy class amidst semi-algebraic obstacles in the plane. In Proc. of the 1998 Int. Symp. on Symbolic and Algebraic Computations (ISSAC’98), pages 17–24. ACM Press, 1998.

[GS99] D. Grigoriev and A. Slissenko. Computing minimum-link path in a homotopy class amidst semi-algebraic obstacles in the plane. St. Petersburg Math. J., 10(2):315–332, 1999.

[H99] S. Har-Peled Constructing approximate shortest path maps in three dimensions SIAM J. Comput., 28(4):1182–1197, 1999.

[HKSS93] J. Heintz, T. Krick, A. Slissenko, and P. Solernó. Une borne inférieure pour la construction de chemins polygonaux dans r^{n}. In Publications du département de mathématiques de l’Université de Limoges, pages 94–100. l’Université de Limoges, France, 1993.

[HKSS94] J. Heintz, T. Krick, A. Slissenko, and P. Solernó. Search for shortest path around semialgebraic obstacles in the plane. J. Math. Sciences, 70(4):1944–1949, 1994. Translation into English of the paper published in Zapiski Nauchn. Semin. LOMI, vol. 192(1991), p. 163–173.

[Mit00] J. S. B. Mitchell. Geometric shortest paths and network optimization. In Jörg-Rüdiger Sack and Jorge Urrutia, editors, Handbook of Computational Geometry, pages 633–701. Elsevier Science Publishers B.V. North-Holland, Amsterdam, 2000.

[MMP87] J.S.B. Mitchell, D. Mount, and C. H. Papadimitriou. The discrete geodesic problem. SIAM J. Comput., 16(4):647–668, 1987.

[MP91] J. Mitchell and C. Papadimitriou. The weighted region problem: finding shortest paths through a weighted planar subdivision. J. ACM, 38:18–73, 1991.

[Pap77] C. H. Papadimitriou. Euclidean TSP is NP-complete. Theor. Comput. Sci., 4:237–244, 1977.

[Pap85] C. H. Papadimitriou. An algorithm for shortest-path motion in three dimensions. Inform. Process. Lett., 20:259–263, 1985.

[SS90] J. T. Schwartz and M. Sharir. Algorithmic motion planning in robotics. In J. van Leeuwen, editor, Handbook of Theoretical Computer Science. Vol. A, pages 391–430. Elsevier Science Publishers, B. V., 1990.

[SS86] M. Sharir, and A. Schorr. On shortest paths in polyhedral spaces. SIAM J. Comput., 15:193– 215, 1986.