The condition number of Riemannian approximation problems
Paul Breiding Nick Vannieuwenhoven


Multigrid deflation for Lattice QCD
    [
    Accepted 2019 September 24; received 2019 September 10; in original form 2019 July 12.
==========================================================================================




 
We consider the local sensitivity of least-squares formulations of inverse problems. The sets of inputs and outputs of these problems are assumed to have the structures of Riemannian manifolds.
The problems we consider include the approximation problem of finding the nearest point on a Riemannian embedded submanifold from a given point in the ambient space.
We characterize the first-order sensitivity, i.e., condition number, of local minimizers and critical points to arbitrary perturbations of the input of the least-squares problem.
This condition number involves the Weingarten map of the input manifold, which measures the amount by which the input manifold curves in its ambient space.
We validate our main results through experiments with the n-camera triangulation  problem in computer vision.



 Riemannian least-squares problem, sensitivity, condition number, local minimizers, Weingarten map, second fundamental form



 90C31, 53A55, 65H10, 65J05, 65D18, 65D19






§ INTRODUCTION
 
A prototypical instance of the Riemannian optimization problems we consider here arises from the following inverse problem (IP): Given an input x∈ I, compute y∈ O such that F(y)=x. Herein, we are given a system modeled as a smooth map F between a Riemannian manifold O, the output manifold, and a Riemannian embedded submanifold I⊂^n, the input manifold. In practice, however, we are often given a point a∈^n in the ambient space rather than x∈I, for example because a was obtained from noisy or inaccurate physical measurements. This leads to the Riemannian least squares optimization problem associated to the foregoing IP:

    PIP_y ∈O1/2 F(y) - a ^2.

In <cit.> we called this type of problem a  parameter identification problem (PIP).
Application areas where PIPs originate are data completion with low-rank matrix and tensor decompositions <cit.>, geometric modeling <cit.>, computer vision <cit.>, and phase retrieval problems in signal processing <cit.>.

The manifold O is called the output manifold because its elements are the outputs of the optimization problem. The outputs are the quantities of interest; for instance, they could be the unknown system parameters that one wants to identify based on a given observation a and the forward description of the system F.
In order to distinguish between the inputs x∈I of the IP and the inputs a∈^n of the optimization formulation <ref>, we refer to the latter as ambient inputs.
In this paper, we study the sensitivity of the output of the optimization problem <ref> (and an implicit generalization thereof) when the ambient input a∈^n is perturbed by considering its condition number κ at global and local minimizers and critical points. 






 §.§ The condition number
 
Condition numbers are one of the cornerstones in numerical analysis<cit.>.
Its definition was introduced for Riemannian manifolds already in the early days of numerical analysis. In 1966, Rice <cit.> presented the definition in the case where input and output space are general nonlinear metric spaces; the special case of Riemannian manifolds is <cit.>.
Let I and O be Riemannian manifolds and assume that the map

f: I→ O models a computational problem. Then, Rice's definition of the condition number of f at a point x∈I is

    κ[f](x) := 
    lim_ϵ→ 0 sup_x' ∈I,
    _I(x,x') ≤ϵ _O(f(x),f(x'))/_I(x,x'),

where _I(·,·) denotes the distance given by the Riemannian metric on I and likewise for O.
It follows immediately that the condition number yields an asymptotically sharp upper bound on the perturbation of the output of a computational problem when subjected to a perturbation of the input x. Indeed, we have

    _O( f(x), f(x') ) ≤κ[f](x) ·_I( x, x' ) + o( _I( x, x' ) ).

The number κ[f](x) is an absolute condition number. On the other hand, if I is a submanifold of a vector space equipped with the norm ·_I and analogously for O and ·_O, then the relative condition number is κ_rel[f](x):= κ[f](x)  ‖ x ‖_I‖ f(x)‖_O. We focus on computing the absolute condition number, as the relative condition number can be inferred from the absolute one.



 §.§ Condition number of approximation problems

Our starting point is Rice's classic definition <ref>, which is often easy to compute for differentiable maps, rather than Demmel and Renegar's inverse distance to ill-posedness, which is generally more difficult to compute. Nevertheless, Rice's definition cannot be applied unreservedly to measure how sensitive a solution of <ref> is to perturbations of the ambient input a∈^n. At first glance, it seems that we could apply <ref> to

    ϕ : ^n →O,  a ↦_y∈O‖ F(y) - a ‖^2.

However, this is an ill-posed expression, as ϕ is not a map in general! Indeed, there can be several global minimizers for an ambient input a∈^n and their number can depend on a. Even if we could restrict the domain so that ϕ becomes a map, the foregoing formulation has two additional limitations.
First, when <ref> is a nonconvex Riemannian optimization problem, it is unrealistic to expect that we can solve it globally. Instead we generally find local minimizers whose sensitivity to perturbations of the input is also of interest. For our theory to be practical, we need to consider the condition number of every local minimizer separately. Exceeding this goal, we will present the condition number of each critical point in this paper. 
Second, requiring that the systems in <ref> can be modeled by a smooth map F : O→I excludes some interesting computational problems. For example, computing the eigenvalues of a symmetric matrix cannot be modeled as a map from outputs to inputs: there are infinitely many symmetric matrices with a prescribed set of eigenvalues. 

The solution to foregoing complications consists of allowing multivalued or set-valued maps g : I⇉O in which one input x∈I can produce a set of outputs Y⊂O <cit.>. Such maps are defined implicitly by the graph 

    G = { (x,y) ∈I×O| y ∈ g(x) }.
 
In numerical analysis, Wilkinson <cit.> and Woźniakowski <cit.> realized that the condition number of single-valued localizations or local maps of G could still be studied using the standard techniques.
Blum, Cucker, Shub, and Smale <cit.> developed a geometric framework for studying the condition number of single-valued localizations when the graph G is an algebraic subvariety. More recently, Bürgisser and Cucker <cit.> considered condition numbers when the graph G is a smooth manifold. Because these graphs are topologically severely restricted, the condition number at an input–output pair (x,y)∈ G can be computed efficiently using only linear algebra.

The heart of this paper consists of characterizing the condition number of critical points of PIPs (and an implicit generalization thereof) by employing the geometric approach to conditioning from <cit.>. Therefore, we summarize its main points next.




 §.§ Implicit smooth computational problems

Following <cit.>, we model an inverse problem with input manifold I and output manifold O implicitly by a solution manifold:

    S⊂I×O.

Given an input x∈I, the problem is finding a corresponding output y ∈ O such that (x,y) ∈S.
As S models a computational problem, we can consider its condition number.
Now it is defined at (x,y)∈S and not only at the input x, because the sensitivity to perturbations in the output depends on which output is considered.

The condition number of (the computational problem modeled by) S is derived in <cit.>. Let π_I:S→I and π_O:S→O be the projections on the input and output manifolds, respectively. Recall that <cit.> assume that
dim S =dim I,
because otherwise almost all inputs x have infinitely many outputs (dim S > dim I) or no outputs (dim S < dim I). Furthermore, we will assume that π_I is surjective, which means that every input x∈I has at least one output.
When the derivative[See <ref> for a definition.] π_I(x,y) is invertible, the inverse function theorem <cit.> implies that π_I is locally invertible. Hence, there is a local smooth solution map π_O∘π_I^-1, which locally around (x,y) makes the computational problem explicit. This local map has condition number κ[π_O∘π_I^-1](x) from <ref>. Moreover, because the map is differentiable, Rice gives an expression for <ref> in terms of the spectral norm of the differential in <cit.>. In summary, <cit.> conclude that the condition number of S at (x,y)∈S is

    κ[S](x,y)   =
    κ[π_O∘π_I^-1](x)      if π_I(x,y) is invertible,
    ∞    otherwise
       =  (π_O(x,y)) (π_I(x,y))^-1_I→O,

where the latter spectral norm is defined below in <ref>.
Note that if S is the graph of a smooth map f:I→O, then κ[S](x,f(x))=κ[f](x).


 By definition, <ref> evaluates to ∞ if π_I(x,y) is not invertible.


The points in 
W := {(x,y)∈S|κ[S](x,y) < ∞} are called well-posed tuples. Points in S\ W are called ill-posed tuples. An infinitesimal perturbation in the input x of an ill-posed tuple (x,y) can cause an arbitrary change in the output y. It is therefore natural to assume that W≠∅, because otherwise the computational problem is always ill-posed. If there exists (x,y)∈ W, it follows from the inverse function theorem that W is an open submanifold of S, so S =  W. We also assume that W is dense. Otherwise, since S is Hausdorff, around any data point outside of W there would be an open neighborhood on which the computational problem is ill-posed; we believe these input-output pairs should be removed from the definition of the problem.



In summary, we make the following assumption throughout this paper.


The m-dimensional input manifold I is an embedded smooth submanifold of ^n equipped with the Euclidean inner product inherited from ^n as Riemannian metric. The output manifold O is a smooth Riemannian manifold. The solution manifold S⊂I×O is a smoothly embedded m-dimensional submanifold. The projection π_I : S→I is surjective. The well-posed tuples W form an open dense embedded submanifold of S.


This assumption is typically not too stringent. It is satisfied for example if S is the graph of a smooth map f : I→O or F : O→I (as for IPs) by <cit.>. Another example is when S is the smooth locus of an algebraic subvariety of I×O, as most problems in <cit.>.
Moreover, the condition of an input–output pair (x,y) is a local property. Since every immersed submanifold is locally embedded <cit.>, taking a restriction of the computational problem S to a neighborhood of (x,y), in its topology as immersed submanifold, yields an embedded submanifold. The assumption and proposed theory always apply in this local sense.



 §.§ Related work

The condition number <ref> was computed for several computational problems. Already in the 1960s, the condition numbers of solving linear systems, eigenproblems, linear least-squares problems, and polynomial rootfinding, among others were determined <cit.>. An astute observation by Demmel <cit.> in 1987 paved the way for the introduction of condition numbers in optimization. He realized that many condition numbers from linear algebra can be interpreted as an inverse distance to the nearest ill-posed problem. This characterization generalizes to contexts where an explicit map as required in <ref> is not available.  Renegar <cit.> exploited this observation to determine the condition number of linear programming instances and he connected it to the complexity of interior-point methods for solving them. This ignited further research into the sensitivity of optimization problems, leading to (Renegar's) condition numbers for linear programming <cit.>, conic optimization <cit.>, and convex optimization <cit.> among others. In many cases, the computational complexity of methods for solving these problems can be linked to the condition number; see <cit.>. In our earlier work <cit.>, we connected the complexity of Riemannian Gauss–Newton methods for solving <ref> to the condition number of the exact inverse problem F(y)=x. We anticipate that the analysis from <cit.> could be refined to reveal the dependence on the condition number of the optimization problem <ref>. 


Beside Renegar's condition number, other measures of sensitivity are studied in optimization. The sensitivity of a set-valued map g : I⇉O at (x,y) with y∈ g(x) is often measured by the Lipschitz modulus <cit.>. This paper only studies the sensitivity of a single-valued localization ĝ of the set-valued map g, in which case the Lipschitz modulus K of g at (x,y) is a Lipschitz constant of ĝ <cit.>. Consequently, Rice's condition number κ[ĝ](x), which we characterize in this paper for Riemannian approximation problems, is always upper bounded by K; see, e.g., <cit.>.



 §.§ Main contributions

We characterize the condition number of computing global minimizers, local minimizers, and critical points of the following implicit generalization of <ref>:

    IPIPπ_O∘_(x,y)∈S1/2 a - π_I(x,y) ^2,

where a ∈^n is a given ambient input and the norm is the Euclidean norm. We call this implicit formulation of <ref> a Riemannian implicit PIP.

Our main contributions are <ref>. The condition number of critical points of  the Riemannian approximation problem min_x∈O a - x ^2, is treated in <ref>. 
<Ref> contains the condition number of <ref> at global minimizers when a is sufficiently close to the input manifold. Finally, the most general case is <ref>: it gives the condition number of <ref> at all critical points including local and global minimizers.

Aforementioned theorems will show that the way in which the input manifold I curves in its ambient space ^n, as measured by classic differential invariants, enters the picture.
To the best of our knowledge, the role of curvature in the theory of condition of smooth computational problems involving approximation was not previously identified. We see this is as one of the main contributions of this article.



 §.§ Approximation versus idealized problems

The results of this paper raise a pertinent practical question: If one has a computational problem modeled by a solution manifold S that is solved via an optimization formulation as in <ref>, which condition number is appropriate? The one of the idealized problem S, where the inputs are constrained to I? Or the more complicated condition number of the optimization problem, where we allow ambient inputs from I's ambient space ^n?
We argue that the choice is only an illusion.

The first approach should be taken if the problem is defined implicitly by a solution manifold S⊂ I× O and the ambient inputs are perturbations a∈^n of inputs x∈I. The idealized framework will still apply because of <ref>.
Consider for example the n-camera triangulation problem, discussed in <ref>, where the problem S_MV is defined by computing the inverse of the camera projection process.
If the measurements a were taken with a perfect pinhole camera for which the theoretical model is exact, then even when a∉I_MV (e.g., errors caused by pixelization), as long as it is a sufficiently small perturbation of some x∈ I_MV, the first approach is suited.

On the other hand, if it is only postulated that the computational problem can be well approximated by the idealized problem S, then the newly proposed framework for Riemannian approximation problems is the more appropriate choice. In this case, the “perturbation” to the input is a modeling error.
An example of this is applying n-camera triangulation to images obtained by cameras with lenses, i.e., most real-world cameras.
In this case, the computational problem is not S_MV, but rather a more complicated problem that takes into account lens distortions. When using S_MV as a proxy, as is usual <cit.>, the computational problem becomes an optimization problem. Curvature should be taken into account then, as developed here.








 §.§ Outline

Before developing the main theory, we briefly recall the necessary concepts from differential geometry and fix the notation in the next section. A key role will be played by the Riemannian Hessian of the distance from an ambient input in ^n to the input manifold I. Its relation to classic differential-geometric objects such as the Weingarten map is stated in <ref>. The condition number of critical points of <ref> is first characterized in <ref> for the special case where the solution manifold represents the identity map. The general case is studied in <Ref>; the former deriving the condition number of global minimizers, while the latter treats critical points.
The expressions of the condition numbers in <ref> are obtained from implicit formulations of the Riemannian IPIPs. 
For clarity, in the case of local minimizers, we also express them as condition numbers of global minimizers of localized Riemannian optimization problems in <ref>. 
<Ref> explains a few standard techniques for computing the condition number in practice along with their computational complexity. Finally, an application of our theory to the triangulation problem in multiview geometry is presented in <ref>.



§ PRELIMINARIES
  

We briefly recall the main concepts from Riemannian geometry that we will use. The notation introduced here will be used throughout the paper.
Proofs of the fundamental results appearing below can be found in <cit.>.



 §.§ Manifolds
 
By smooth m-dimensional manifold M we mean a C^∞ topological manifold that is second-countable, Hausdorff, and locally Euclidean of dimension m. The tangent space pM is the m-dimensional linear subspace of differential operators at p. A differential operator at p∈M is a linear map v_p from the vector space C^∞(M) over  of smooth functions f : M→ to  that satisfies the product rule:
v_p (f · g) = (v_p f) · g(p) + f(p) · (v_p g)
for all f, g ∈ C^∞(M). If M⊂^n is embedded, the tangent space can be identified with the linear span of all tangent vectors d/dt|_t=0γ(t) where γ⊂M is a smooth curve passing through p at 0. 



A differentiable map F : M→N between manifolds M and N induces a linear map between tangent spaces Fp:pM→F(p)M, called the derivative of F at p. If v_p∈pM, then w_F(p) := (Fp)(v_p)
is the differential operator w_F(p)(f) := v_p(f∘ F) for all f∈ C^∞( N).

The identity map 1_M : M→M will be denoted by 1, the space on which it acts being clear from the context.



 §.§ Tubular neighborhoods

Let M^m ⊂^n be an embedded submanifold. The disjoint union of tangent spaces to M^m is the tangent bundle of M:

    TM := ∐_p∈MpM = { (p, v) | p ∈M and  v ∈pM}.

The tangent bundle is a smooth manifold of dimension 2m. The normal bundle is constructed similarly. The normal space of M in ^n at p is the orthogonal complement in ^n of the tangent space pM, namely N_p M := ( p^n )^⊥. The normal bundle is a smooth embedded submanifold of ℝ^n ×ℝ^n of dimension n. 
Formally, we write

    NM := ∐_p ∈MN_pM = { (p, η) | p ∈M and η∈N_pM}.

Let δ : M→_+ be a positive, continuous function. Consider the following open neighborhood of the normal bundle:

    V_δ = { (p, η_p) ∈NM|η_p < δ(p) },

where the norm is the one induced from the Riemannian metric; in our setting it is the standard norm on ^n. The map E : NM→^n, (p,η) ↦ p + η is smooth and there exists a δ such that the restriction E|_V_δ becomes a diffeomorphism onto its image <cit.>. Consequently, T = E(V_δ) is an n-dimensional, open, smooth, embedded submanifold of ^n that forms a neighborhood of M. The submanifold T is called a tubular neighborhood of M and δ is called its height (function).



 §.§ Vector fields

A smooth vector field is a smooth map X : M→TM from the manifold to its tangent bundle taking a point p to a tangent vector v_p. By the notation X|_p we mean X(p). A smooth vector field X on a properly embedded submanifold M⊂N can be extended to a smooth vector field X on N such that it agrees on M: X|_M = X. 


A smooth frame of M^m is a tuple of m smooth vector fields (X_1, …, X_m) that are linearly independent: (X_1|_p, …, X_m|_p) is linearly independent for all p ∈M. If these tangent vectors are orthonormal for all p, then the frame is called orthonormal.

Let X be a smooth vector field on M. An integral curve of X is a smooth curve γ⊂M such that γt = X|_γ(t) for all t ∈ (-1,1). For every p ∈M the vector field X generates an integral curve γ with starting point p = γ(0). There is always a smooth curve realizing a tangent vector v_p∈pM, i.e., a curve γ with p = γ(0) and v = γ0.





























 §.§ Riemannian manifolds


A Riemannian manifold is a smooth manifold M equip­ped with a Riemannian metric g: a positive definite symmetric bilinear form g_p:pM×pM→ that varies smoothly with p. The metric g_p induces a norm on pM denoted ‖ t_p ‖_ M := √(g_p(t_p, t_p)).
The only Riemannian metric we will explicitly use in computations is the standard Riemannian metric of ^n, i.e., the Euclidean inner product g_p(x,y) = ⟨ x, y⟩ = x^T y.

If M and N are Riemannian manifolds with induced norms ‖·‖_ M and ‖·‖_ N, and if F :M→N is a differentiable map, then the spectral norm of Fp is

    Fp_M→N := sup_t_p ∈pM∖{0}‖ (Fp)(t_p) ‖_ N/‖  t_p‖_ M.

If the manifolds are clear from the context we can drop the subscript and write Fp.



 §.§ Riemannian gradient and Hessian


Let M be a Riemannian manifold with metric g, and let f : M→ be a smooth function. The Riemannian gradient of f at p∈M is the tangent vector ∇ f|_p ∈pM that is dual to the derivative fp : pM→ under g, i.e., g_p(∇ f|_p, t) = (fp)(t) for all t∈pM. It is known that ∇ f is a vector field, hence explaining the notation ∇ f|_p for its value at p∈M.

When M⊂^n is a Riemannian embedded submanifold with the inherited Euclidean metric, the Riemannian Hessian <cit.> can be defined without explicit reference to connections, as in <cit.>:

    Hess_p (f) : pM→pM, η_p ↦P_pM( (∇ fp)(η_p) ),

where ∇ f now has to be viewed as a vector field on ^n and P_pM projects orthogonally onto the tangent space pM. In other words, the Hessian takes η_p to P_pM( d/dt|_t=0∇ f|_γ(t))
where γ(t) ⊂M is a curve realizing the tangent vector η_p.

The Hessian Hess_p(f) is a symmetric bilinear form on the tangent space pM whose (real) eigenvalues contain essential information about the objective function f at p. If the Hessian is positive definite, negative definite, or indefinite at a critical point p of f (i.e., ∇ f|_p = 0 or, equivalently, fp=0), then p is respectively a local minimum, local maximum, or saddle point of f.




§ MEASURING CURVATURE WITH THE WEINGARTEN MAP
 

A central role in this paper is played by the Riemannian Hessian of the 
distance function

    d : I→,  x ↦1/2x - a^2,

which measures the Euclidean distance between a fixed point a ∈^n and points on the embedded submanifold I⊂^n. It is the inverse of this Riemannian Hessian matrix that will appear as an additional factor to be multiplied with π_O(x,y) (π_I(x,y))^-1 from <ref> in the condition numbers of the problems studied in <ref>.

The Riemannian Hessian of d contains classic differential-geometric information about the way I curves inside of ^n <cit.>. Indeed, from <cit.> we conclude that it can be expressed at x∈I as

    Hess_x (d) = 1_xI + _x( P_xM (x-a) ),

where P_xM is the projection onto the normal space, and
_x : xM→ S^2( xM ) is a classical differential invariant called the second fundamental form <cit.>, which takes a normal vector in xM and sends it to a symmetric bilinear form on the tangent space xM. If we let η_x := a - x, then at a critical point of d we have η_x ∈xI, so that the Hessian of interest satisfies

    H
     H_η := Hess_x(d) = 1 - S_η,  where  S_η := _x(η_x).

Note that we dropped the base point x∈I from the notation in H_η and S_η. In the above, S_η is interpreted as a symmetric endomorphism on xI and it is classically called the shape operator or Weingarten map <cit.>.

The Weingarten map measures how much the embedding into the ambient space ^n bends the manifold I^m in the following concrete way <cit.>.
For η≠ 0, let w:=η‖η‖ and let c_1,…,c_m be the eigenvalues of S_w = 1‖η‖ S_η. As S_w is self-adjoint, all c_i's are all real. The eigenvalues of S_η are then c_1‖η‖ ,…,c_m‖η‖. The c_i's are called the principal curvatures of I at x in direction of η. They have the following classic geometric meaning <cit.>: If u_i∈T_x I is a unit length eigenvector of the eigenvalue c_i,
then locally around x and in the direction of u_i the manifold I contains an infinitesimal arc of a circle with center x+c_i^-1w =x+(c_i‖η‖)^-1η through x. This circle is called an osculating circle. The radii r_i := | c_i|^-1 of those circles are called the critical radii of I at x in direction η. See also <ref>.

For several submanifolds I⊂^n relevant in optimization applications, the Weingarten map S_η was computed. Absil, Mahoney, and Trumpf <cit.> provide expressions for the sphere, Stiefel manifold, orthogonal group, Grassmann manifold, and low-rank matrix manifold. Feppon and Lermusiaux <cit.> further extended these results by computing the Weingarten maps and explicit expressions for the principal curvatures of the isospectral manifold and “bi–Grassmann” manifold, among others.
Heidel and Schultz <cit.> computed the Weingarten map of the manifold of higher-order tensors of fixed multilinear rank <cit.>.




§ THE CONDITION NUMBER OF CRITICAL POINT PROBLEMS
 


Prior to treating <ref> in general, we focus on an important and simpler special case.
This case arises when the solution manifold S⊂I×O is the graph of the identity map 1 : I→I, so O =I. Given an ambient input a ∈^n, the problem <ref> now comprises finding an output x ∈I so that the distance function d from <ref> is minimized. This is can be formulated as the standard Riemannian least-squares problem 

    min_x ∈I1/2 x - a ^2.

Computing global minimizers of this optimization problem is usually too ambitious a task when I is a complicated nonlinear manifold. Most Riemannian optimization methods will converge to local minimizers by seeking to satisfy the first-order necessary optimality conditions of <ref>, namely a-x∈N_x I. We say that such methods attempt to solve the critical point problem (CPP) whose implicit formulation is

    CPPS_CPP := { (a,x) ∈ℝ^n ×I|  a - x ∈N_xI}.


We call S_CPP the critical point locus. It is an exercise to show that S_CPP is linearly diffeomorphic to the normal bundle NI so that we have the following result.

 
S_CPP is an embedded submanifold of ^n×I of dimension n.


Consequently, the CPP falls into the realm of the theory of condition from <ref>. We denote the coordinate projections of S_CPP by Π_^n: S_CPP→ℝ^n and
 Π_I: S_CPP→I for distinguishing them from the projections associated to S. The corresponding submanifold of well-posed tuples is

    W_CPP := { (a,x) ∈S_CPP|Π_^n(a,x) is invertible}.

The fact that it is a manifold is the following result, proved in <ref>.

 W_CPP is an open dense embedded submanifold of S_CPP.


We show in <ref> that the condition number of this problem can be characterized succinctly as follows.[We take the convention that whenever an expression ·_I→O is ill-posed it evaluates to ∞.]


Let (a,x)∈S_CPP and η:=a-x. Then,

    κ[S_CPP](a,x) = ‖ H_η^-1‖_I→I,

where H_η is given by <ref>. Moreover, if (a,x) ∈S_CPP∖W_CPP, then H_η is not invertible.



κ[S_CPP] : S_CPP→∪{+∞} is continuous and κ[S_CPP](x,y) is finite if and only if (x,y)∈W_CPP.


Curvature has entered the picture in <ref> in the form of the Riemannian Hessian H_η from <ref>.
It is an extrinsic property of I, however, so that the condition number depends on the Riemannian embedding of I into ^n. Intuitively this is coherent, as different embeddings give rise to different computational problems each with their own sensitivity to perturbations in the ambient space.

Let us investigate the meaning of the foregoing theorem a little further. Recall, from <ref>, the definition of the principal curvatures c_i as the eigenvalues of the Weingarten map S_w where w = η/η and η=a-x. Then, <ref> implies that

    κ[S_CPP](a,x) = max_1≤ i≤ m1/| 1-c_i‖η‖|,

so its condition is completely determined by the amount by which I curves inside of ^n and the distance of the ambient input a to the critical point x∈I.
We see that κ[S_CPP](a,x) = ∞ if and only if I contains an infinitesimal arc of the circle with center x+η = a (so that there is an i with c_i ‖η‖ = 1). Consequently, if η lies close to any of the critical radii r_i = |c_i|^-1 then the CPP is ill-conditioned at (a,x).

The formula from <ref> makes it apparent that for some inputs we may have that κ[S_CPP](a,x)<1. In other words, for some ambient inputs a∈ℝ^n the infinitesimal error ‖Δ a‖ can shrink! <Ref> gives a geometric explanation of when the error is amplified (panel (a)) and when it is shrunk (panels (b) and (c)).



We also have an interpretation of κ[S_CPP](a,x) as a normalized inverse distance to ill-posedness: for (a,x)∈S_CPP, let η=a-x be the normal vector that connects x to a.
The locus of ill-posed inputs for the CPP is Σ_CPP = S_CPP∖W_CPP. The projection Π_^n(Σ_CPP) is well studied in the literature. Thom <cit.> calls it the target envelope of S_CPP, while Porteous <cit.> calls it the focal set of I, and for curves in the plane it is called the evolute; see Chapter 10 of <cit.>. 

The points in the intersection of {x}× (x+η) with Σ_CPP are given by the multiples of the normal vector η whose length equals one of the critical radii r_1,…,r_m; this is also visible in <ref>. Comparing with <ref> we find

    1/κ[S_CPP](a,x) =  min_η_0∈η:   (x,a+η_0) ∈Σ_CPPη_0/| ‖η_0‖ - ‖η‖ |,
     where η = a-x.

This is an interpretation of Rice's condition number in the spirit of Demmel <cit.> and Renegar <cit.>.


If I is a smooth algebraic variety in ^n, then solving S_CPP reduces to solving a system of polynomial equations. For almost all inputs a∈^n this system has a constant number of solutions over the complex numbers, called the Euclidean distance degree <cit.>. <Ref> thus characterizes the condition number of real solutions of this special system of polynomial equations.





§ THE CONDITION NUMBER OF APPROXIMATION PROBLEMS
 
Next, we analyse the condition number of global minimizers of <ref>. This setting will be crucial for a finer interpretation of <ref> in <ref>.

By suitably choosing the height function δ <cit.>, we can assume for each ambient input a in a tubular neighborhood T⊂^n of I that there is a unique point x on I that minimizes the distance from I to a. 





In this case, the computational problem that consists of finding global minimizers of <ref> can be modeled implicitly with the graph

    APS_AP := {(a,y) ∈ T ×O| (P_I(a),y)∈S},

 where P_I :  T →I, a↦argmin_x∈I‖ x-a‖ is the nonlinear projection that maps a ∈T to the unique point on I that minimizes the distance to a. We refer to the problem S_AP as the approximation problem (AP) to distinguish it from the problem of finding critical points of <ref>, which is treated in <ref>.

If S_AP were a smooth submanifold of T ×O, its condition number would be given by <ref>. However, this is not the case in general. We can nevertheless consider the subset of
well-posed tuples for <ref>:

    W_AP := {(a,y) ∈ T ×O| (P_I(a),y)∈W}.

In <ref>, we prove the following result under <ref>.

 
W_AP is a smooth embedded submanifold of T×O of dimension n. Moreover, W_AP is dense in S_AP.


The condition number of <ref> is now given by <ref> on W_AP. For the remaining points, i.e., S_AP∖W_AP, the usual definition <ref> still applies. From this, we obtained the following characterization, which is proved in <ref>.

 
Let the ambient input a ∈T⊂^n be sufficiently close to I, specifically lying in a tubular neighborhood of I, so that it has a unique closest point x ∈I. Let η=a-x. Then,

    κ[S_AP](a,y) =  (π_O(x,y))(π_I(x,y))^-1 H_η^-1_I→O,

where π_I and π_O are the projections onto the input and output manifolds respectively, and H_η is as in <ref>. Moreover, if (a,y)∈S_AP∖W_AP then π_I(x,y) is not invertible.



κ[S_AP] : T→∪{+∞} is continuous and κ[S_AP](a,y) is finite if and only if (a,y)∈W_AP.


 
If f : I→O is a smooth map and S is the graph of f, then the foregoing specializes to

    (π_O(x,y))(π_I(x,y))^-1 = fx,   so  κ[S_AP](x,y) =  (fx) H_η^-1_I→O

Similarly, if F : O→I is a smooth map with graph S, then 

    (π_O(x,y))(π_I(x,y))^-1 = (Fy)^-1,   so  κ[S_AP](x,y) =  (Fy)^-1 H_η^-1_I→O

by the inverse function theorem.
These observations apply to <ref> as well.


The contribution (π_O(x,y))(π_I(x,y))^-1 to the condition number originates from the solution manifold S; it is intrinsic and does not depend on the embedding. On the other hand, the curvature factor H_η^-1 is extrinsic in the sense that it depends on the embedding of I into ^n, as in <ref>.

The curvature factor in <ref> disappears when S_η = 0. This occurs when the ambient input a lies on the input manifold, so η=0. In this case, H_0^-1=1, so we have the following result.

Let (x,y)∈ S. Then, we have κ[S](x,y)=κ[S_AP](x,y).





In other words, if the ambient input, given by its coordinates in ^n, lies exactly on I,
the condition number of <ref> equals the condition number of the computational problem S, even though there can be many more directions of perturbation in ^n than in I!
By continuity, the effect of curvature on the sensitivity of the computational problem can essentially be ignored for small ‖η‖, i.e., for points very close to I. This is convenient in practice because computing the Weingarten map is often more complicated than obtaining the derivatives of the projection maps π_O and π_I.




§ THE CONDITION NUMBER OF GENERAL CRITICAL POINT PROBLEMS
 
We now consider the sensitivity of critical points of <ref>. 
As mentioned before, our motivation stems from solving <ref> when a closed-form solution is lacking. Applying Riemannian optimization methods to them, one can in general only guarantee to find points satisfying the first-order optimality conditions <cit.>. Consequently, in theory, they solve a generalized critical point problem: Given an ambient input a ∈^n, produce an output y ∈O such that there is a corresponding x ∈I that satisfies the first-order optimality conditions of <ref>.
We can implicitly formulate this computational problem as follows:

    GCPPS_GCPP := {(a,x,y) ∈^n ×I×O| (a,x)∈S_CPP and   (x,y)∈S}.



The graph S_GCPP is defined using three factors, where the first is the input and the third is the output. The second factor I encodes how the output is obtained from the input. Triples are ill-posed if either (x,y)∈S is ill-posed for the original problem or (a,x)∈S_CPP is ill-posed for the <ref>. Therefore, the well-posed locus is

    W_GCPP := {(a,x,y) ∈^n ×I×O| (a,x)∈W_CPP and   (x,y)∈W}.

We prove in <ref> that the well-posed tuples form a manifold.

 
W_GCPP is an embedded submanifold of ^n ×W of dimension n.
Moreover, W_GCPP is dense in S_GCPP.


Consequently, on W_GCPP we can use the definition of condition from <ref>. 
Our final main result generalizes <Ref>; it is proved in <ref>.



Let (a,x,y)∈S_GCPP and η=a-x. Then, we have

    κ[S_GCPP](a,x,y) = ‖ (π_O(x,y) )   (π_I(x,y))^-1 H_η^-1‖_I→O,

where π_I:S→I and π_O:S→O are coordinate projections, and
H_η is given by <ref>. Moreover, if (a,x,y)∈S_GCPP∖W_GCPP, then either π_I(x,y) or H_η is not invertible.



κ[S_GCPP] : S_GCPP→∪{+∞} is continuous and the condition number κ[S_GCPP](a,x,y) is finite if and only if (a,x,y)∈W_GCPP.



Since I inherits the Euclidean metric from ^n, the following result is immediately obtained from <ref>.

 
Let (a,x,y)∈S_GCPP and η=a-x. We have

    κ[S](x,y)/max_1≤ i≤ m| 1-c_i‖η‖|≤κ[S_GCPP](a,x,y)
    ≤κ[S](x,y)/min_1≤ i≤ m| 1-c_i‖η‖|,

where c_1,…,c_m are the principal curvatures of I^m at x in direction η=a-x.


Since κ[S](x,y) = ‖ (π_ O(x,y)) (π_ I(x,y))^-1‖_I→O, the submultiplicativity of the spectral norm both yields
κ[S](x,y)
≤κ[S_GCPP](a,x,y) ‖ H_η‖
and also
κ[S_GCPP](a,x,y)
≤κ[S](x,y)‖ H_η^-1‖,
concluding the proof.




§ CONSEQUENCES FOR RIEMANNIAN OPTIMIZATION


In <ref>, we established a connection between the CPP and the Riemannian optimization problem <ref>: S_CPP is the manifold of critical tuples (a,x) ∈^n ×I, such that x is a critical point of the squared distance function from I to a.
For the general implicit formulation in <ref>, however, the connection to Riemannian optimization is more complicated. For clarity, we therefore also characterize the condition number of local minimizers of <ref> as the condition number of the unique global minimizer of a localized Riemannian optimization problem.

The solutions of the GCPP were defined as the generalization of the AP to critical points. However, we did not show that they too can be seen as the critical tuples of a distance function optimized over in a Riemannian optimization problem. This connection is established next and proved in <ref>.

 
(a,x,y) ∈W_GCPP if and only if (a,(x,y)) is a well-posed critical tuple of the function optimized over in the Riemannian optimization problem

    min_(x,y) ∈W1/2 a - π_I(x,y) ^2.



Our main motivation for studying critical points instead of only the global minimizer as in <ref> stems from our desire to predict the sensitivity of outputs of Riemannian optimization methods for solving approximation problems like <ref>. When applied to a well-posed optimization problem, these methods in general only guarantee convergence to critical points <cit.>. However, in practice, convergence to critical points that are not local minimizers is extremely unlikely <cit.>. The reason is that they are unstable outputs of these algorithms: Applying a tiny perturbation to such critical points will cause the optimization method to escape their vicinity, and with high probability, converge to a local minimizer instead.


One should be careful to separate condition from numerical stability. The former is a property of a problem, the latter the property of an algorithm.
We do not claim that computing critical points of GCPPs other than local minimizers are ill-conditioned problems. We only state that many Riemannian optimization methods are unstable algorithms for computing them via <ref>. This is not a bad property, since the goal was to find minimizers, not critical points!



For critical points of the GCPP that are local minimizers of the problem in <ref> we can characterize their condition number as the condition number of the unique global minimizer of a localized problem. This result is proved in <ref>.

 
Let the input manifold I, the output manifold O, and well-posed loci W and W_GCPP be as in <ref>. Assume that (a^⋆,x^⋆,y^⋆) ∈W_GCPP and x^⋆ is a local minimizer of d_a^⋆ :I→, x ↦1/2 a^⋆ - x ^2.
Then, there exist open neighborhoods A_a^⋆⊂^n around a^⋆ and  N_(x^⋆,y^⋆)⊂W around (x^⋆,y^⋆), so that

    ρ_(a^⋆,x^⋆,y^⋆) : A_a^⋆→O,   a ↦π_O∘_(x,y) ∈N_(x^⋆,y^⋆)1/2 a - π_I(x,y) ^2

is an <ref> and the condition number of this smooth map satisfies 

    κ[ρ_(a^⋆,x^⋆,y^⋆)](a) 
    
       = κ[S_GCPP](a, P_I_x^⋆(a), ρ_(a^⋆,x^⋆,y^⋆)(a)),

where I_x^⋆=π_I(N_(x^⋆,y^⋆)) is the projection onto the first factor, and P_I_x^⋆(a) = _x∈I_x^⋆ a - x ^2 is the nonlinear projection onto I_x^⋆.


Taking I=O and S as the graph of the identity map 1, the next result follows.


Let the input manifold I, and well-posed loci W and W_CPP be as in <ref>. Assume that (a^⋆,x^⋆) ∈W_CPP is well-posed and that x^⋆ is a local minimizer of d_a^⋆:I→,  x ↦1/2 a^⋆ - x ^2.
Then, there exist open neighborhoods A_a^⋆⊂^n around a^⋆ and I_x^⋆⊂ I around x^⋆, so that 

    ρ_(a^⋆,x^⋆): A_a^⋆→I,    a ↦_x ∈I_x^⋆1/2 a - x ^2

is an <ref> and its condition number is 
κ[ρ_(a^⋆,x^⋆)](a) 
= κ[S_CPP](a,ρ_(a^⋆,x^⋆)(a)).


An interpretation of these results is as follows. 
The condition number of a generalized critical point (a^⋆,x^⋆,y^⋆) that happens to correspond to a local minimum of <ref> equals the condition number of the map ρ_(a^⋆,x^⋆,y^⋆). This smooth map models a Riemannian AP on the localized computational problem N_(x^⋆,y^⋆)⊂W which has a unique global minimum for all ambient inputs from A_a^⋆. Consequently, every local minimum of <ref> is robust in the sense that it moves smoothly according to ρ_(a^⋆,x^⋆,y^⋆) in neighborhoods of the input a^⋆ and (x^⋆,y^⋆)∈W. In this open neighborhood, the critical point remains a local mimimum.


The foregoing entails that if a local minimum (x^⋆,y^⋆)∈W is computed for <ref> with ambient input a^⋆ using Riemannian optimization, then computing a local minimum (x, y) for <ref> with perturbed input a^⋆ + Δ and starting point (x^⋆,y^⋆) results in dist_O(y,y^⋆) ≤κ[S_GCPP](a^⋆,x^⋆,y^⋆) Δ + o(Δ), where dist_O is as in <ref>, provided that Δ is sufficiently small.




§ COMPUTING THE CONDITION NUMBER
 
Three main ingredients are required for evaluating the condition number of <ref> numerically: 

 
  (i) the inverse of the Riemannian Hessian H_η from <ref>,
 
  (ii) the derivative of the computational problem (π_O(x,y))(π_I(x,y))^-1, and 
 
  (iii) the spectral norm <ref>.

The only additional item for evaluating the condition number of <ref> relative to the idealized inverse problem “given x ∈I, find a y∈O such that (x,y)∈S” is the curvature term H_η^-1 from item (i). Therefore, we focus on evaluating (i) and only briefly mention the standard techniques for (ii) and (iii).

The derivative in (ii) is often available analytically as in <ref>. If no analytical expression can be derived, it can be numerically approximated when O⊂^q is an embedded submanifold by computing the standard Jacobian of S, viewed as the graph of a map ^n→^q, and composing it with the projections onto xI and yO. 

Evaluating the spectral norm in (iii) is a linear algebra problem. Let I⊂^n be embedded with the standard inherited metric and let the Riemannian metric of O be given in coordinates by a matrix G. If G = R^T R is its Cholesky factorization <cit.>, then A _I→O =  R A _2 = σ_1(RA), where A : xI→A(x)O is a linear map and ·_2 is the usual spectral norm, i.e., the largest singular value σ_1(·). The largest singular value can be computed numerically in 𝒪( min{ m^2 p, p^2 m } ) operations, where m=I and p=O, using standard algorithms <cit.>. If the linear map LA can be applied in fewer than 𝒪(m p) operations, a power iteration can be more efficient <cit.>.

The Riemannian Hessian H_η can be computed analytically using standard techniques from differential geometry for computing Weingarten maps <cit.>, or it can be approximated numerically <cit.>. Numerical approximation of the Hessian with forward differences is implemented in the Matlab package Manopt <cit.>. Its inverse can in both cases be computed with 𝒪(m^3) operations.

We found two main analytical approaches helpful to compute Weingarten maps of embedded submanifolds I⊂^n. The first relies on the following characterization of the Weingarten map S_η_x in the direction of a normal vector η_x = a - x. It is defined in <cit.> as

    S_η_x: xI→xI,  v_x ↦P_xI( - (∇_v_x N)|_x ),

where P_xI denotes the orthogonal projection onto xI, N is a smooth extension of η to a normal vector field on an open neighborhood of x in I, and ∇_v_x is the Euclidean covariant derivative <cit.>. In practice, the latter can be computed as

    (∇_v_x Y)|_x = d/d t Y|_γ(t),

where Y is an arbitrary vector field on I and γ(t) is a smooth curve realizing v_x <cit.>. In other words, S_η_x maps η_x to the tangential part of the usual directional derivative of N in direction of v_x ∈xI⊂^n.



The second approach relies on the second fundamental form of I, which is defined as the projection of the Euclidean covariant derivative onto the normal space N_x I:

    _x(X,Y) := ( (X,Y) )|_x :=  P_xI( (∇_X|_x Y)|_x ),

where now both X and Y are vector fields on I. The second fundamental form is symmetric in X and Y, and _x(X,Y) only depends on the tangent vectors X|_x and Y|_x. It can be viewed as a smooth map _x : xM×xM→xM. Contracting the second fundamental form with a normal vector η_x ∈xI yields an alternative definition of the Weingarten map: For all v, w ∈xI we have
⟨ S_η_x (v), w ⟩ = ⟨_x(v,w), η_x ⟩,
where ⟨·,·⟩ is the Euclidean inner product.
Computing the second fundamental form is often facilitated by pushing forward vector fields through a diffeomorphism F : M→N.
In this case, there exists a vector field Y on N that is F-related to a vector field X on M: Y|_p = (F_* X)|_p := (FF^-1(p))(X|_F^-1(p)). The integral curves generated by X and Y = F_* X are related by Proposition 9.6 of <cit.>. This approach will be illustrated in the next section.

As a rule of thumb, we find that computation of the condition number of <ref> is approximately as expensive as one iteration of a Riemannian Newton method for solving this problem. Typically, the cost is not worse than 𝒪(m^3 + n m^2 + min{m^2 p, p^2 m}) operations, where the terms correspond one-to-one to items (i)–(iii), and where I^m ⊂^n and p=O.



§ TRIANGULATION IN COMPUTER VISION
 

A rich source of approximation problems whose condition can be studied with the proposed framework is multiview geometry in computer vision; for an introduction to this domain see <cit.>. The tasks consist of recovering information from r ≥ 2 camera projections of a scene in the world ^3. We consider the pinhole camera model. In this model the image formation process is modeled by the following transformation, as visualised in <ref>:

    μ_r: y↦[ A_ℓy + b_ℓ/c_ℓ^T y + d_ℓ ]_ℓ=1^r,

where A_ℓ∈^2× 3, b_ℓ∈^2, c_ℓ∈^3, and d_ℓ∈; see <cit.>. Clearly, μ_r is not defined on all of ^3.  We clarify its domain in <ref> below.

The vector x = μ_r(y) ∈^2r obtained by this image formation process is called a consistent point correspondence.
Information that can often be identified from consistent point correspondences include camera parameters and scene structure <cit.>.

As an example application of our theory we compute the condition number of the triangulation problem in computer vision <cit.>. In this computational problem we have r ≥ 2 stationary projective cameras and a consistent point correspondence x = (x_1, x_2, …, x_r) ∈^2r. The goal is to retrieve the world point y∈^3 from which they originate. Since the imaging process is subject to noisy measurements and the above idealized projective camera model holds only approximately <cit.>, we expect that instead of x we are only given a = (a_1, a_2, …, a_r) close to x. Thus, x is the true input of the problem, y is the output and a is the ambient input.

According to <cit.> , the “gold standard algorithm” for triangulation solves the (Riemannian) optimization problem

    min_y∈^31/2a - μ_r(y) ^2.

Consequently, the triangulation problem can be cast as a GCCP provided that two conditions hold: (i) the set of consistent point correspondences is an embedded manifold, and (ii) intersecting the back-projected rays is a smooth map from aforementioned manifold to ^3. We verify both conditions in the following subsection.








 §.§ The multiview manifold

The image formation process can be interpreted as a projective transformation from the projective space ^3 to ^2 in terms of the 4 × 3 camera matrices P_ℓ=[[   A_ℓ   b_ℓ; c_ℓ^T   d_ℓ ]] <cit.>.
It yields homogeneous coordinates of x_ℓ = ( z^ℓ_1 / z^ℓ_3, z^ℓ_2 / z^ℓ_3 ) ∈^2 where z^ℓ = P_ℓ[[ y; 1 ]].
Note that if z_3^ℓ = 0 then the point y has no projection[Actually, it has a projection if we allow points at infinity, i.e., if we consider the triangulation problem in projective space.] onto the image plane of the ℓ-th camera, which occurs precisely when y lies on the principal plane of the camera <cit.>. This is the plane parallel to the image plane through the camera center c_ℓ; see <ref>. It is also known that points on the baseline of two cameras, i.e., the line connecting the camera centers visualised by the dashed line in <ref>, all project to the same two points on the two cameras, called the epipoles <cit.>. Such points cannot be triangulated from only two images. For simplicity, let B⊂^3 be the union of the principal planes of the first two cameras and their baseline, so B is a 2-dimensional subvariety. The next result shows that a subset of the consistent point correspondences outside of B forms a smooth embedded submanifold of ^2r.


Let O_MV = ^3 ∖B with B as above.
The map μ_r : O_MV→^2r
is a diffeomorphism onto its image I_MV = μ_r(O_MV).


Clearly μ_r is a smooth map between manifolds (c_ℓ^T y + d_ℓ = z_3^ℓ 0). It only remains to show that it has a smooth inverse.
Theorem 4.1 of <cit.> states that μ_2's projectivization is a birational map, entailing that μ_2 is a diffeomorphism onto its image. Let π_1:4 : ^2r→^4 denote projection onto the first 4 coordinates. Then, μ_2^-1∘π_1:4 is a smooth map such that (μ_2^-1∘π_1:4) ∘μ_r = μ_2^-1∘μ_2 = 1_O_MV, having used that the domain O_MV is the same for all r. For the right inverse, we see that any element of I_MV can be written as μ_r(y) for some y∈O_MV. Hence,

    (μ_r ∘ (μ_2^-1∘π_1:4))(μ_r(y)) = (μ_r ∘1_O_MV)(y) = μ_r(y),

so it is the identity on I_MV. This proves that μ_r has μ_2^-1∘π_1:4 as smooth inverse. 


As μ_r^-1 is required in the statement of the Weingarten map, we give an algorithm for computing μ_r^-1 = μ_2^-1∘π_1:4.


Assume that we are given x∈I_MV⊂^2r in the r-camera multiview manifold and let its first four coordinates be (x_1, y_1, x_2, y_2). Then, by classic results <cit.>, the unique element in the kernel of

    [ (x_1 e_3^T - e_1^T) P_1; (y_1 e_3^T - e_2^T) P_1; (x_2 e_3^T - e_1^T) P_2; (y_2 e_3^T - e_1^T) P_2;                         ]∈^4 × 4

yields homogeneous coordinates of the back-projected point in ^3.

The solution manifold S_MV⊂I_MV×O_MV is the graph of μ_r^-1. Therefore, it is a properly embedded smooth submanifold; see, e.g., <cit.>.




 §.§ The second fundamental form

Vector fields and integral curves on I_MV can be viewed through the lens of μ_r:

We construct a local smooth frame of I_MV by pushing forward a frame from O_MV by μ_r. Then, the integral curves of each of the local smooth vector fields are computed, after which we apply the Gauss formula for curves <cit.> to compute the second fundamental form.

Because of <ref> a local smooth frame for the multiview manifold I_MV⊂^2r is obtained by pushing forward the constant global smooth orthonormal frame ( e_1, e_2, e_3 ) of ^3 by the derivative of μ_r.
Its derivative is

    μ_ry : yO_MV→μ_r(y)I_MV,   ẏ↦[ A_ℓẏ/c_ℓ^T y + d_ℓ - (c_ℓ^T ẏ) A_ℓy + b_ℓ/(c_ℓ^T y + d_ℓ)^2 ]_ℓ=1^r,

and so a local smooth frame of I_MV is given by

    E_i : I_MV→TI_MV, μ_r(y) ↦[ A_ℓe_i/c_ℓ^T y + d_ℓ - c_ℓ,iA_ℓy + b_ℓ/(c_ℓ^T y + d_ℓ)^2 ]_ℓ=1^r,    i=1,2,3,

where c_ℓ,i=c_ℓ^T e_i. It is generally neither orthonormal nor orthogonal.

We compute the second fundamental form of the multiview manifold by differentiation along the integral curves generated by the smooth local frame (E_1,E_2,E_3). The integral curves through μ_r(y) generated by this frame are the images of the integral curves passing through y generated by the e_i's due to <cit.>. The latter are seen to be g_i(t) = y + t e_i by elementary results on linear differential equations. Therefore, the integral curves generated by E_i are

    γ_i(t)
     = μ_r(g_i(t))
     = [ A_ℓ(y + t e_i) + b_ℓ/c_ℓ^T (y + t e_i)  + d_ℓ ]_ℓ=1^r.

The components of the second fundamental form at x = μ_r(y) ∈I_MV are then

    _x( E_i, E_j )
        = P_NI_MV( ddt E_j|_γ_i(t)) 
       = P_NI_MV( ddt[ A_ℓe_j/c_ℓ^T (y + te_i) + d_ℓ - c_ℓ,jA_ℓ (y + t e_i) + b_ℓ/(c_ℓ^T (y + t e_i) + d_ℓ)^2 ]_ℓ=1^r
     ) 
       = P_NI_MV( [ - c_ℓ,i/α^2_ℓ(y) A_ℓe_j - c_ℓ,j/α^2_ℓ(y) A_ℓe_i + 2 c_ℓ,i c_ℓ,j/α_ℓ^3(y) (A_ℓy + b_ℓ) ]_ℓ=1^r  ),

where α_ℓ(y) = c_ℓ^T y + d_ℓ and i,j = 1, 2, 3.



 §.§ A practical algorithm
 

The Weingarten map of I_MV in the direction of the normal vector η∈N_xI_MV is obtained by contracting the second fundamental form with η; that is, S_η = ⟨_x(E_i, E_j), η⟩. This can be computed efficiently using linear algebra operations. Partitioning η = [ η_ℓ ]_ℓ=1^r with η_ℓ∈^2, the symmetric coefficient matrix of the Weingarten map relative to the frame E = (E_1, E_2, E_3) becomes

    S_η =
     ∑_ℓ=1^r [ 
     2 c_ℓ,i c_ℓ,j/α^3_ℓ(y)η_ℓ^T (A_ℓy + b_ℓ) - c_ℓ,i/α^2_ℓ(y)η_ℓ^T A_ℓe_j - c_ℓ,j/α^2_ℓ(y)η_ℓ^T A_ℓe_i ]_i,j=1^3∈^3 × 3,

where y = μ_r^-1(x).


To compute the spectrum of the Weingarten map using efficient linear algebra algorithms, we need to express it with respect to an orthonormal local smooth frame by applying Gram–Schmidt orthogonalization to E. This is accomplished by placing the tangent vectors E_1, E_2, E_3 as columns of a 2r × 3 matrix J and computing its compact QR decomposition Q R = J. The coefficient matrix of the Weingarten map expressed with respect to the orthogonalized frame Q=(Q_1, Q_2, Q_3) is then S_η = R^-TS_η R^-1. 


Since μ_r : O_MV→I_MV is a diffeomorphism, (π_O_MV(x,y)) (π_I_MV(x,y))^-1 equals the inverse of the derivative of μ_r by <ref>. As R is the matrix of μ_ry : yO_MV→xI_MV expressed with respect to the orthogonal basis (Q_1,Q_2,Q_3) of I_MV and the standard basis of ^3, we get that κ[S_GCPP](x+η,x,y) equals

    (π_O_MV(x,y)) (π_I_MV(x,y))^-1 H_η^-1_I→O
     =  R^-1 (I - S_η)^-1_2 
     = 1/σ_3( (I - S_η)R ),

where I is the 3 × 3 identity matrix, and σ_3 is the third largest singular value.

The computational complexity of this algorithm grows linearly with the number of cameras r. Indeed, the Weingarten map S_η can be constructed in 𝒪(r) operations, the smooth frame E is orthogonalized in 𝒪(r) operations, the change of basis from S_η to S_η requires a constant number of operations, and the singular values of a 3×3 matrix adds another constant.




 §.§ Numerical experiments

The computations below were implemented in Matlab R2017b <cit.>. The code we used is provided as supplementary files accompanying the arXiv version of this article. It uses functionality from Matlab's optimization toolbox. The experiments were performed on a computer running Ubuntu 18.04.3 LTS that consisted of an Intel Core i7-5600U CPU with 8GB main memory.


We present a few numerical examples illustrating the behavior of the condition number of triangulation. The basic setup is described next. We take all 10 camera matrices P_i ∈^3 × 4 from the “model house” data set of the Visual Geometry Group of the University of Oxford <cit.>. These cameras are all pointing roughly in the same direction. In our experiments we reconstruct the point y := p + 1.5^10v∈^3, where
p= (-1.85213, -0.532959, -5.65752)∈^3 is one point of the house and v=(-0.29292,-0.08800, -0.95208) is the unit-norm vector that points from the center of the first camera p.
Given a data point a ∈^2r, it is triangulated as follows. First a linear triangulation method is applied, finding the right singular vector corresponding to the least singular value of a matrix whose construction is described (for two points) in <cit.>. We then use it as a starting point for solving optimization problem <ref> with Matlab's nonlinear least squares solver . For this solver, the following settings were used:  and  set to 10^-28,  equal to 10^-14, and both  and  set to 10^4. We provided the Jacobian to the algorithm, which it can evaluate numerically.

Since the multiview manifold I_MV is only a three-dimensional submanifold of ^2r and the computational complexity is linear in r, we do not focus on the computational efficiency of computing κ[S_GCPP]. In practical cases, with r ≪ 1000, the performance is excellent. For example, the average time to compute the condition number for r=1000 cameras was less than 0.1 seconds in our setup. This includes the time to compute the frame E, the Weingarten map expressed with respect to E, the orthogonalization of the frame, and the computation of the smallest singular value.




  
Experiment 1
In the first experiment, we provide numerical evidence for the correctness of our theoretical results. We project y to the image planes of the r=10 cameras: x = μ_r(y). A random normal direction η∈xI_MV is sampled by taking a random vector with i.i.d. standard normal entries, projecting it to the normal space, and then scaling it to unit norm. We investigate the sensitivity of the points along the ray a(t) := x - t η.

The theoretical condition number of the triangulation problem at a(t) can be computed numerically for every t using the algorithm from <ref>. We can also estimate the condition number experimentally, for example by generating a large number of small perturbations a(t) + E(t) with E(t)≈0, solving <ref> numerically using  and then checking the distance to the true critical point y. However, from the theory we know that the worst direction of infinitesimal perturbation is given by the left singular vector of (I - S_-tη)R corresponding the smallest (i.e. third) singular value, where R is as in <ref>. Let u(t) ∈^3 denote this vector, which contains the coordinates of the perturbation relative to the orthonormal basis Q=(Q_1,Q_2,Q_3); see <ref>. Then, it suffices to consider only a small (relative) perturbation in this direction; we took E(t) := 10^-6a(t) Q u(t). We solve <ref> with input a(t) + E(t) and output y_est(t) using  with the exact solution of the unperturbed problem y as starting point. The experimental estimate of the condition number is then κ_est(t) = y - y_est(t) / E(t).



In <ref> we show the experimental results for the above setup, where we chose t = 10^i x with 100 values of i uniformly spaced between -3 and 2. It is visually evident that the experimental results support the theory. Excluding the three last values of t, the arithmetic and geometric means of the condition number divided by the experimental estimate are approximately 1.006664 and 1.006497 respectively. This indicates a near perfect match with the theory.

There appears to be one anomaly in <ref>, however: the three crosses to the right of the singularity at t^⋆ = 79.64416 x. What happens is that past the singularity, y is no longer a local minimizer of the optimization problem min_y∈O_MV1/2 a(t) - μ_10(y) ^2. This is because x(t) is no longer a local minimizer of min_x ∈I_MV1/2 a(t) - x ^2.
Since we are applying an optimization method, the local minimizer happens to move away from the nearby critical point and instead converges to a different point.
Indeed, one can show that the Riemannian Hessians (see <cit.>) of these problems are positive definite for small t, but become indefinite past the singularity t^*.









  
Experiment 2
The next experiment illustrates how the condition number varies with the distance t along a normal direction t η∈NI_MV. We consider the setup from the previous paragraph. The projection of the point y onto the image planes of the first k cameras is x_k := μ_k(y)=(P_1 y, …, P_k y). A random unit-norm normal vector η_k in N_x_kI_MV is chosen as described before. We investigate the sensitivity of the perturbed inputs a_k(t) := x_k + t η_k for t ∈. The condition number is computed for these inputs for k=2,3,5,10 and ± t/x_k = 10^i for 10^4 values of i uniformly spaced between -3 and 4. The results are shown in <ref>. We have 3 peaks for a fixed k because I_MV = 3. The two-camera system has only two singularities in the range |t| < 10^4 x_2. These figures illustrate the continuity of κ_GCPP from <ref>.

Up to about |t| ≤ 10^-2x_k we observe that curvature hardly affects the condition number of the GCPP for the triangulation problem. Beyond this value, curvature plays the dominant role: both the singularities and the tapering off for high relative errors are caused by it. Moreover, the condition number decreases as the number of cameras increases. This indicates that the use of minimal problems <cit.> in computer vision could potentially introduce numerical instability.

 



§ ACKNOWLEDGEMENTS

The authors thank Nicolas Boumal for posing a question on unconstrained errors, which eventually started the project that led to this article. We are greatly indebted to Carlos Beltrán who meticulously studied the manuscript and asked pertinent questions about the theory.
Furthermore, the authors would like to thank Peter Bürgisser for hosting the second author at TU Berlin in August 2019, and Joeri Van der Veken for fruitful suggestions and pointers to the literature.

We thank two anonymous reviewers whose remarks enticed us to immensely improve the presentation of our results, specifically in <ref>, and to include a discussion of computational aspects in <ref>. We also thank the editor, Anthony So, for adeptly handling the reviewing process. 



§ PROOFS OF THE TECHNICAL RESULTS
 
We recall two classic results that we could not locate in modern references in differential geometry.
A common theme in the three computational problems considered in <ref> is that they involve the projection
Π_I : NI→I,    (x,η) ↦ x
from the normal bundle NI to the base of the bundle, the embedded submanifold I⊂^n. 
This projection is smooth <cit.> and its derivative is characterized by the following classical result, which dates back at least to Weyl <cit.>.


  Let (x,η)∈N I and (ẋ, η̇)∈T_(x,η)N I. Then,
  P_T_x I(ẋ+η̇) = H_ηẋ,
  where H_η is as in <ref>.































A useful consequence is that it allows us to characterize when H_η is singular.

 
 Let (a,x) ∈S_CPP, η = a - x, and H_η be as in <ref>. Then, H_η is invertible if and only if (a,x)∈W_CPP.


 Let (ȧ,ẋ) ∈(a,x)S_CPP. Then, we have ẋ = (Π_I(a,x))(ȧ,ẋ) and ȧ = (Π_^n(a,x))(ȧ,ẋ), so that <ref> yields the equality of operators

    P_xI Π_^n(a,x) = H_η Π_I(a,x).

 
On the one hand, if (a,x) ∈W_CPP then the map Π_^n(a,x) is invertible. By multiplying with its inverse we get P_xI = H_η (Π_I(a,x)) (Π_^n(a,x))^-1. As H_η is an endomorphism on xI, which is the range of the left-hand side, the equality requires H_η to be an automorphism, i.e., it is invertible.

On the other hand, if (a,x)∉W_CPP, there is some (ȧ,ẋ)≠ 0 in the kernel of Π_^n(a,x). Since ȧ = (Π_^n(a,x))(ȧ,ẋ) = 0, we must have ẋ≠ 0. Consequently, applying both sides of <ref> to (ȧ,ẋ) gives 0 = H_ηẋ, which implies that H_η is singular. This finishes the proof. 


The next result follows almost directly from <ref> and appears several times as an independent statement in the literature. The earliest reference we could locate is Abatzoglou <cit.>, who presented it for C^2 embedded submanifolds of ^n in a local coordinate chart. 


 
Let I⊂^n be a Riemannian embedded submanifold. Then, there exists a tubular neighborhood T of I such that the projection map P_I : T→I,  a ↦_x ∈I a - x ^2 is a smooth submersion with derivative

    P_Ia = H_η^-1P_xI,

where x = P_I(a), η = a-x, and H_η is given by <ref>.








In the next subsections, proofs are supplied for the main results.



 §.§ Proofs for <ref>
 









Recall that Π_^n is a smooth map, so if its derivative is nonsingular at a point (a,x)∈S_CPP, then there exists an open neighborhood of (a,x) where this property remains valid. Therefore, applying the inverse function theorem <cit.> at (a,x) ∈W_CPP, the fact that Π_^n(a,x) is invertible implies that there is a neighborhood N_(a,x) of (a,x) in W_CPP and a neighborhood of a in ^n such that N_(a,x) is diffeomorphic to (an open ball in) ^n. Hence, W_CPP is open in S_CPP, and, hence, it is an open submanifold of dimension n.

To show that W_CPP is dense in S_CPP, take any (a,x) ∈S_CPP∖W_CPP. Let η = a-x. Then, (x,x) + (αη,0) ∈S_CPP because η∈xI. By <ref>, (x,x) + (αη,0) ∈W_CPP with α∈(0,∞) if and only if H_αη=1 - S_αη = 1 - α S_η is singular. 

This occurs only if 1/α equals one of the eigenvalues of the Weingarten map S_η. Consequently, (a,x) can be reached as the limit of a sequence (x,x)+ (α_n η,0), which is wholly contained in W_CPP by taking α_n outside of the discrete set.

Applying <cit.>, we can conclude that W_CPP is even an embedded submanifold. This concludes the argument.




Let (a,x) ∈W_CPP be arbitrary.
On W_CPP, the coordinate projection Π_^n: S_CPP→ℝ^n has an invertible derivative. Consequently, by the inverse function theorem <cit.> there exist open neighborhoods X⊂S_CPP of (a,x) and Y of a⊂^n
such that Π_^n|_X:  X→ Y has a smooth inverse function that we call ϕ_(a,x). Its derivative is ϕ_(a,x)a = (Π_^n(a,x))^-1. Consider the smooth map Φ_(a,x) := Π_I∘ϕ_(a,x), where Π_I: S_CPP→I is the coordinate projection. The CPP condition number at (a,x)∈W_CPP is κ[Φ_(a,x)](a) = κ[W_CPP](a,x)= ‖Φ_(a,x)a‖_^n→I. Since Φ_(a,x)a = (Π_I(a,x)) (Π_^n(a,x))^-1, we have that
κ[W_CPP](a,x) = ‖ (Π_I(a,x)) (Π_^n(a,x))^-1‖_^n→I. Now it follows from <ref> that (Π_I(a,x)) (Π_^n(a,x))^-1 = H_η^-1P_xI. 
As the metrics on I and ^n are identical, H_η^-1P_xI_^n→I = H_η^-1_I→I.
This concludes the first part.

The second part is <ref>.





 §.§ Proofs for <ref>
 


By <ref> the projection P_I : T→I is a smooth map on T. Let us denote the graph of this map by V := { (a, P_I(a) ) ∈T×I| a∈ T}. It is a smooth embedded submanifold of dimension T= N I = n; see, e.g., <cit.>. Let x=P_I(a), η = a-x and (ȧ, ẋ) ∈(a,x)V be a nonzero tangent vector.
Because of <ref> it satisfies ẋ = H_η^-1P_xIȧ. Hence,
(Π_^n(a,x))(ȧ, H_η^-1P_xIȧ ) = ȧ 0,
so that its kernel is trivial. It follows that V⊂W_CPP. Since their dimensions match by <ref>, V is an open submanifold. As W_CPP is embedded by <ref>, V is embedded as well by <cit.>.
Moreover, by construction, we have
 
    W_AP = (V×O) ∩ (^n ×W) ⊂^n ×I×O.

The first part of the proof is concluded by observing that the proof of <ref> applies by replacing W_CPP with V and W_GCPP with W_AP.

Finally, we show that W_AP is dense in S_AP. Let (x+η,y)∈S_AP with x∈I and η∈N_x I be arbitrary. Then P_I(x+η) = x and (x,y)∈S.
As W is an open dense submanifold of S, there exists a sequence (x_i,y_i)∈W such that lim_i→∞ (x_i, y_i) = (x,y). Consider (x_i + η_i, y_i) ∈W_AP with η_i ∈N_x_iI and η_i →η chosen in such a way that the first component lies in the tubular neighborhood; this is possible as T is an open submanifold of ^n. Then, the limit of this sequence is (x+η,y). Q.E.D.



For (a,y) ∈W_AP, let x=P_I(a) and η=a-x. By <ref>, κ[W_AP](a,y) = ‖(π_O∘π_I^-1∘P_ℐa)‖_^n→O.
Since (x,y) ∈W, we know that π_I(x,y) is invertible. Moreover, by <ref> the derivative of the projection P_I :  T→I at a is P_Ia = H_η^-1P_xI.
Therefore, 

    κ[W_AP](a,y) = ‖ (π_O(x,y)) (π_I(x,y) )^-1H_η^-1‖_I→O.


For (a,y) ∈S_AP∖W_AP, we have (x,y)∉W with x=P_I(a). By definition of W, this entails that π_I(x,y) is not invertible, concluding the proof.













 §.§ Proofs for <ref>
 


Combining <ref>, it is seen that W_GCPP is realized as the intersection of two smoothly embedded product manifolds:
 
    W_GCPP = (W_CPP×O) ∩ (^n ×W) ⊂^n ×I×O.

Note that codim(W_CPP×O) = I, by <ref>, and codim(^n ×W) = O, because W = S and <ref> stating that S = I. Thus, it suffices to prove that the intersection is transversal <cit.>, because then W_GCPP would be an embedded submanifold of codimension I + O.

Consider a point (a,x,y) ∈W_GCPP. Transversality means that we need to show

    (a,x,y)(W_CPP×O) + (a,x,y)(^n ×W)=(a,x,y)(^n ×I×O).

Fix an arbitrary (ȧ, ẋ, ẏ) ∈a^n×xI×yO. Then, it suffices to show there exist ( (ȧ_1, ẋ_1), ẏ_1) ∈(a,x)W_CPP×yO and ( ȧ_2, (ẋ_2, ẏ_2))∈a^n×(x,y)W such that
(ȧ, ẋ, ẏ) = (ȧ_1 + ȧ_2, ẋ_1 + ẋ_2, ẏ_1 + ẏ_2).
It follows from <ref> that ȧ_1 = H_ηẋ_1 + η̇,
where η=a-x and where η̇∈N_x I. The tangent vectors ẏ_1 ∈yO, ȧ_2 ∈a^n and  ẋ_1 ∈xI can be chosen freely.
Choosing
ẋ_1 = ẋ and ẋ_2 = 0; ẏ_1 = ẏ and ẏ_2 = 0 and ȧ_1 =
H_ηẋ_1 and ȧ_2 = ȧ - ȧ_1 yields (ȧ, ẋ, ẏ). This concludes the first part of the proof.

It remains to prove W_GCPP is dense in S_GCPP. For this, we recall that

    S_GCPP = (S_CPP×O) ∩ (^n ×S).

Let (a,x,y)∈S_GCPP and let η = a-x∈N_x I. Then, (a,x) ∈S_CPP and (x,y)∈S.
As W is an open dense submanifold of S by <ref>, there exists a sequence (x_i,y_i)∈W such that lim_i→∞ (x_i, y_i) = (x,y). Moreover, there exists a sequence η_i →η with η_i ∈x_iI so that (x_i + η_i, x_i) ∈S_CPP, because Π_ I(S_CPP) =  I. Since W_CPP is open dense in S_CPP
by <ref>, only finitely many elements of the sequence are not in W_CPP. We can pass to a subsequence (x_j+η_j,x_j)∈W_CPP. Now (x_j + η_j, x_j, y_j) ∈W_GCPP and moreover its limit is (x+η,x,y) ∈S_GCPP. Q.E.D.



For a well-posed triple (a,x,y)∈W_GCPP we obtain from <ref> that κ[W_GCPP](a,x,y) = ‖ (Π_O(a,x,y)) (Π_^n(a,x,y))^-1‖_^n→O. We compute the right-hand side of this equation.
The situation appears as follows:

    W_CPP^Π_^n[r]_Π_I[d]    ^n    W_GCPP^Π_O[d]_Π_^n[l]   
    I   W^π_I[l]_π_O[r]    O

We have in addition the projections (1×Π_I):W_GCPP→W_CPP,  (a,x,y)↦ (a,x) and (Π_I×1):W_GCPP→W,  (a,x,y)↦ (x,y). With this we have

    (Π_^n∘(1×π_I))(a,x,y) = Π_^n(a,x,y)   and   (π_O∘(Π_I×1))(a,x,y) = Π_O(a,x,y).

Taking derivatives we get

    (Π_^n(a,x)) (1×π_I)(a,x,y)  = Π_^n(a,x,y) and 
    (π_O(a,x)) (Π_I×1)(a,x,y)  = Π_O(a,x,y).

Since (x,y)∈W, the derivative π_I(x,y) is invertible, and so (1×π_I)(a,x,y) is also invertible. On other hand, since (a,x)∈W_CPP, the derivative Π_^n(a,x) is invertible.
Altogether we see that Π_^n(a,x,y) is invertible,
and that

    (Π_O(a,x,y)) (Π_^n(a,x,y))^-1
    
       =  (π_O(x,y)) ( (Π_I×1)(a,x,y)) ((1×π_I)(a,x,y))^-1  (Π_^n(a,x))^-1.

The derivatives in the middle satisfy

    (ẋ, ẏ)    =(Π_I×1)(a,x,y) ( (1×π_I)(a,x,y))^-1 (ȧ, ẋ),

while on the other hand 
(ẋ, ẏ) = (π_I(x,y))^-1 (Π_I(a,x)) (ȧ, ẋ)
 .
This implies that the two linear maps on the right hand side in the previous equations are equal.
Hence,

    (Π_O(a,x,y)) (Π_^n(a,x,y))^-1   = (π_O(x,y)) (π_I(x,y))^-1 (Π_I(a,x)) ( Π_^n(a,x) )^-1
       = (π_O(x,y)) (π_I(x,y))^-1 H_η^-1P_T_x I,

where we used <ref>. Taking spectral norms on both sides completes the proof for points in W_GCPP.
Finally, for points outside W_GCPP the proof follows from the definition of W_GCPP combined with <ref>.






 §.§ Proofs for <ref>







Taking the derivative of the objective function at a point (x,y)∈W we see that the critical points satisfy ⟨ a - x, π_I(x,y)(ẋ, ẏ) ⟩ = 0
for all (ẋ, ẏ) ∈(x,y)W. By the definition of W, the derivative  π_I(x,y) is surjective, which implies that (π_I(x,y)) ((x,y)W) = xI. Therefore, the condition of being a critical point is equivalent to a - x∈N_x I.




By assumption (a^⋆,x^⋆,y^⋆)∈W_GCPP. In particular, this implies (a^⋆,x^⋆)∈W_GCPP, so that Π_^n(a^⋆,x^⋆):T_(a^⋆,x^⋆)W_CPP→^n is invertible. By the inverse function theorem <cit.>, there is an open neighborhood U_(a^⋆,x^⋆)⊂W_CPP of (a^⋆,x^⋆) such that Π_^n restricts to a diffeomorphism on

    A_a^⋆=Π_^n(U_(a^⋆,x^⋆)).

Let Π_^n^-1 denote the inverse of Π_^n on A_a^⋆.
Combining <ref>, the derivative of Π_ I∘Π_^n^-1 is given by
(Π_I(a,x)) (Π_^n(a,x))^-1 
= H_η^-1P_xI
and thus has constant rank on  A_a^⋆ equal to dim I.

This implies that

    (1× (π_O(x,y)) (π_I(x,y))^-1)(Π_I(a,x)) (Π_^n(a,x))^-1: T_a  A_a^⋆→T_(x,y)W

has constant rank on A_a^⋆
equal to dim I =dim W, by <ref>. Consequently, there is a smooth function
Φ:  A_a^⋆→W
with Φ(a^⋆) = (x^⋆,y^⋆). The above showed that its derivative has constant rank equal to W, so Φ is a smooth submersion.
By Proposition 4.28 of <cit.>, Φ(A_a^⋆) ⊂W is an open subset, so it is a submanifold of dimension I = W. Similarly, the derivative π_I(x,y) at (x,y)∈W has maximal rank equal to I, by definition of W. Hence, π_I∘Φ is also a smooth submersion.
After possibly passing to an open subset of A_a^⋆, we can thus assume

    N_(x^⋆,y^⋆) :=Φ( A_a^⋆)⊂W  and   I_x^⋆ = π_ I( N_(x^⋆,y^⋆)) ⊂ I

are open submanifolds each of dimension dim  I and that

 
  (i) x^⋆ is the global minimizer of the distance to a^⋆ on I_x^⋆, and hence,
 
  (ii) all x∈I_x^⋆ lie on the same side of x^⋆I_x^⋆, considering the tangent space x^⋆I_x^⋆ as an affine linear subspace of ^n with base point x^⋆.

Now let a ∈A_a^⋆ be a point in the neighborhood of a^⋆.
The restricted squared distance function d_a :  N_(x^⋆,y^⋆)→,  (x,y) ↦1/2 a - π_I(x,y) ^2 that is minimized on N_(x^⋆,y^⋆) is smooth.
Its critical points satisfy

    ⟨ a - π_I(x,y), π_I(x,y)(ẋ,ẏ) ⟩ = 0    for all   (ẋ,ẏ) ∈(x,y)N_(x^⋆,y^⋆).

Since (x^⋆,y^⋆) ∈N_(x^⋆,y^⋆)⊂W, the derivative π_I(x^⋆,y^⋆) of π_I : W→I is invertible.
Then, for all (x,y)∈N_(x^⋆,y^⋆) the image of π_I(x,y) is the whole tangent space xI. Consequently, the critical points must satisfy a - x ⊥xI.


By construction, for every a ∈A_a^⋆ there is a unique (x,y) ∈N_(x^⋆,y^⋆) so that (a,x,y) ∈W_GCPP, namely (x,y)=Φ(a).
This implies that (x,y) is the unique critical point of the squared distance function restricted to N_(x^⋆,y^⋆).



It only remains to show that for all a∈A_a^⋆, the minimizer of d_a is attained in the interior of N_(x^⋆,y^⋆). This entails that the unique critical point is also the unique global minimizer on N_(x^⋆,y^⋆).
We show that there exists a δ > 0 such that restricting A_a^⋆ to the open ball B_δ(a^⋆) of radius δ centered at a^⋆ in ^n yields the desired result.

Here is the geometric setup for the rest of the proof, depicted in <ref>: let a∈ B_δ(a^⋆) and let (x,y)∈ N_(x^⋆,y^⋆) be the minimizer for d_a in the closure of N_(x^⋆,y^⋆). Then, x is in the closure of I_x^⋆ by continuity of π_ I. Furthermore, let x'=P(x), where P := P_x^⋆I_x^⋆ is the orthogonal projection onto the tangent space x^⋆I^⋆ considered as an affine linear subspace of ^n with base point x^⋆.

The image of B_δ(a^⋆) under P is an open ball B_δ(x^⋆) ∩x^⋆I_x^⋆⊂x^⋆I_x^⋆. The distance function h : I_x^⋆→,   x ↦x - P( x ) measures the “height” of x from the tangent space x^⋆I_x^⋆.
Now, we can bound the distance from a to x as follows.

    a - x    =  a - a^⋆ + a^⋆ - x^⋆ + x^⋆ - x' + x'  - x
       ≤ a - a^⋆ +  a^⋆ - x^⋆ +  x^⋆ - x'  +  x- x' 
       ≤ a^⋆ - x^⋆ + 2 δ + h(x) 
       ≤ a^⋆ - x^⋆ + 2 δ + max_z ∈P_δ(x^⋆) h(z),

where the first bound is the triangle inequality, and where P_δ(x^⋆) ⊂I_x^⋆ is the preimage of B_δ(x^⋆)∩x^⋆I_x^⋆ under the projection P restricted to I_x^⋆.

Now let x∈∂I_x^⋆ be a point in the boundary of I_x^⋆ and, again, x'=P(x).  We want to show that the distance ‖ a-x‖ is strictly larger than <ref>. This would imply that local minimizer for d_a is indeed contained in the interior of N_(x^⋆,y^⋆).
We have

    a - x    =  a - a^⋆ + a^⋆ - x^⋆ + x^⋆ - x '+ x' - x
       ≥ a^⋆ - x^⋆ + x^⋆ - x' + x' - x  -  a - a^⋆.

Note that x^⋆ = P(a^⋆). Hence we have the orthogonal decomposition

    a^⋆ - x^⋆ + x^⋆ - x' + x' - x = (1 - P)(a^⋆ - x) + P(a^⋆ - x) = a^⋆ - x.

Plugging this into the above, we find

    a - x   ≥ (1 - P)(a^⋆ - x) + P(a^⋆ - x)  - δ
       = √((1 - P)(a^⋆ - x)^2 +  P(a^⋆ - x) ^2 ) - δ
       ≥√((1 - P)(a^⋆ - x)^2+ ϵ^2 ) - δ,

where ϵ = min_z∈ P(∂I_x^⋆)‖ x^⋆ - z‖ > 0.

Observe that a^⋆ - x^⋆ points from x^⋆ to a^⋆ and likewise x' - x points from x to x'. Therefore, ⟨ a^⋆ - x^⋆, x' - x⟩≥ 0 because both point away from the manifold I_x^⋆, which lies entirely on the side of x by assumption (ii) above; see also <ref>. Consequently,
(1 - P)(a^⋆ - x)  = ‖ a^⋆ - x^⋆ +x' - x  ‖≥ a^⋆ - x^⋆,
which implies

    a - x    ≥√(a^⋆ - x^⋆^2 + ϵ^2 ) - δ = a^⋆ - x^⋆ + ϵ^2/2 a^⋆ - x^⋆ - δ + 𝒪(ϵ^4).

Note that h is a continuous function with lim_x→ x^⋆ h(x) = 0. Hence, as δ→0 the maximum height tends to zero as well. Moreover, ϵ and a^⋆ - x^⋆ are independent of δ. All this implies that we can choose δ >0 sufficiently small so that

    ϵ^2/2 a^⋆ - x^⋆ - δ + 𝒪(ϵ^4)≥ 2 δ + max_z ∈P_δ(x^⋆) h(z).

It follows that ‖ a-x‖ is lower bounded by  <ref>, so that the minimizer for d_a must be contained in the interior of N_(x^⋆,y^⋆).



The foregoing shows that on B_δ(a^⋆) ∩A_a^⋆, the map Π_I∘Π_^n^-1 equals the nonlinear projection P_I_x^⋆ onto the manifold I_x^⋆.
Putting everything together, we obtain

    ρ_(a^⋆,x^⋆,y^⋆)(a)
     = (Π_O∘Π_^n^-1) (a)
     = (π_O∘π_I^-1∘Π_I∘Π_^n^-1) (a)
     = (π_O∘π_I^-1∘P_I_x^⋆) (a);

the second equality following <ref>. The condition number of this smooth map is given by <ref>. Comparing with <ref> concludes the proof.

 
siam



