On a Geometrical Interpolation Problem

Octavian Stănăşilă

Abstract

Given a finite set of points A in ℝ^{n} and a geometrical pattern Φ we define a type of distance between A and Φ and study how to find a pattern Φ among a parametrized family of geometrical objects such that such a distance is minimum.

Mathematics Subject Classification: 41A29 Key words: pattern, distance, point–pattern, least squares interpolation problem, Newton–Raphson method.

1 Preliminaries

Many bidimensional ”patterns” have a concatenation of plane curve arcs as boundaries; similarly, 3-dimensional ”patterns” have as boundary a concatenation of surface pieces. For this reason, it could be useful to introduce and to study some new metrical properties, which permit to formulate and solve some interpolation problems.

Let U ⊂ ℝ^{n} be an open set and fix an integer k ≥ 1. Consider a map f : U → ℝ^{k} of class C^{1}(U), f = (f_{1}, ..., f_{k})^{T}; for any x ∈ U, denote J_{f}(x) = (\frac{∂f_{i}}{∂x_{j}} ); 1 ≤ i ≤ k, 1 ≤ j ≤ n, the Jacobian matrix. The most important case in what follows will be when k < n; in this case, a point x ∈ U is said regular for f if rank(J_{f}(x)) = k. The set of type f^{−1}(0) generalize the plane curves, the surfaces, space curves etc. Definition 1.1. We call a pattern in U any set of the form Φ = ⋃_{i=1}^{M} Φ_{i}, such that there are C^{1} functions g_{i} : U → ℝ^{k}, 1 ≤ i ≤ M and Φ_{i} ⊂ g_{i}^{−1}(0) have all their points regular.

If k = 1, Φ ⊂ g^{−1}(0), where g = g_{1} · g_{2}...g_{M}; if x is regular for g_{i}, 1 ≤ i ≤ M, it could be not regular for g.

Given a point a ∈ U, a = (a_{1}, ..., a_{n})^{T} a pattern Φ as above, it could be useful to define a suitable distance d(a, Φ) between them. If so, take A_{1}, A_{2}, ..., A_{N} ∈ U as N distinct points (”points of surveillance”); one knows that even the simple interpolation Lagrange problem has some obstructions to be solved, instead it could be useful to determine, under some conditions, a function f such that the sum ∑_{i=1}^{N} d(A_{i}, f^{−1}(0))^{2} is minimum. Such problems will be treated in what follows.

Balkan Journal of Geometry and Its Applications, Vol.1, No.2, 1996, pp. 97-103 \begin{matrix}c \\ ◯ \end{matrix}Balkan Society of Geometers, Geometry Balkan Press

On the other hand, recall that if A ∈ M_{k,n}(ℝ), B ∈ M_{k,1}(ℝ) are two matrices, then A^{+} ∈ M_{n,k} (ℝ) means the Penrose pseudoinverse of A, [4] and by putting B^{+} = AA^{+}B, the vector ξ = A^{+}B is unique in ℝ^{n} such that Aξ = B^{+}, the latter being just the orthogonal projection of B on ImA (that is ∥ B − B^{+} ∥ = minimum); ξ is called the pseudosolution of the linear system AX = B. If k ≤ n and rankA = n, then A^{+} = A^{T} · (AA^{T})^{−1}.

2 A distance between a point and a pattern

Definition 2.1. Let U ⊂ ℝ^{n} be open and f : U → ℝ^{k}, f = (f_{1}, ..., f_{k})^{T} be a C^{1}-map; for any point a ∈ U define the distance from a to the pattern f^{−1}(0) as being

(1)  δ(a, f^{−1}(0)) =∥ J_{f}(a)^{+} · f(a) ∥ (Euclidean norm)

Example. In the case when k = 1 and f is linear nonnull, f(x) = ∑_{i=1}^{n} c_{i}x_{i}, put c = (c_{1}, ..., c_{n})^{T}; then for any a ∈ ℝ^{n}, a = (a_{1}, ..., a_{n})^{T}, one has f(a) = c^{T} · a, J_{f}(a) = c^{T} and J_{f}(a)^{+} = \frac{c}{∥c∥^{2}} . Put p = a − J_{f}(a)^{+} · f(a) hence p = a − \frac{1}{∥c∥^{2}} (c · c^{T} · a). Then f(p) = c^{T} · p = 0 and the vector a − p is normal to the hyperplane f^{−1}(0). In this case, ∥ a − p ∥= δ(a, f^{−1}(0)), justifying thus the definition 2.1. So δ extends the Euclidean distance. Proposition 2.2. Suppose that 1 ≤ k < n and a ∈ U is a regular point for a C^{1}-map f : U → ℝ^{k}. Then

(2)  δ(a, f^{−1}(0)) = (f(a)^{T} · (J_{f}(a) · J_{f}(a)^{T} )^{−1} · f(a))^{\frac{1}{2}} .

Proof. Let J = J_{f}(a) hence rankJ = k (maximum); in this case, the symmetrical matrix J · J^{T} is invertible and moreover, J^{+} = J^{T} · (J · J^{T})^{−1}. Then by the definition 2.1.,

δ(a, f^{−1}(0))^{2} =∥ J^{+} ·f(a) ∥^{2}=< J^{+} ·f(a), J^{+} ·f(a) >= f(a)^{T} ·(J^{+})^{T} ·J^{+} ·f(a).

Since J^{+} = J^{T} (J · J^{T})^{−1}, one obtains:

δ(a, f^{−1}(0))^{2} = f(a)^{T} · (J · J^{T})^{−1} · J · J^{T} · (J · J^{T})^{−1} · f(a).

But J · J^{T} · (J · J^{T})^{−1} = I_{n} therefore δ(a, f^{−1}(0))^{2} = f(a)^{T} · (J · J^{T})^{−1} · f(a), whence the proposition. Corollary 2.3. Fix a ∈ U. The map f → δ(a, f^{−1}(0)), restricted to C^{1}– functions on U for which a is regular, is continuous. Corollary 2.4. Let f_{i} : ℝ^{n} → ℝ, f_{i}(x) = ∑_{j=1}^{n} a_{ij}x_{j}, where ∑_{j=1}^{n} a_{ij}^{2} = 1 for 1 ≤ i ≤ k < n and the k × n-matrix A = (a_{ij}) has the rank k. Consider the hyperplanes

H_{i} = {x ∈ ℝ^{n}|f_{i}(x) = b_{i}}, 1 ≤ i ≤ k

and put B = (b_{1}, ..., b_{k})^{T}. Then the pseudosolution of the linear system A·X = B minimizes the sum ∑_{i=1}^{k} δ(x, H_{i})^{2}. Proof. Define f : ℝ^{n} → ℝ^{k}, f(x) = (f_{1}(x) − b_{1}, ..., f_{k}(x) − b_{k}) hence by proposition 2.2 δ(x, f^{−1}(0)) = f(x)^{T} · (AA^{T})^{−1} · f(x) for any x ∈ ℝ^{n}. Put B^{+} = AA^{+}B = (b_{1}^{+}, ..., b_{k}^{+})^{T}; then ξ = A^{+}B is the pseudosolution of the system AX = B and Aξ = AA^{+}B = B^{+}. We have

δ(x, H_{i}) = d(x, H_{i}) = |f_{i}(x) − b_{i}|

and

δ(ξ, H_{i}) = d(ξ, H_{i}) = |f_{i}(ξ) − b_{i}|,

hence

∑_{i=1}^{k} δ(x, H_{i})^{2} = d(f(x), b)^{2} ≥ d(Imf, b)^{2} ≥∥ B − B^{+} ∥^{2}=

= ∑_{i=1}^{k} (b_{i}^{+} − b_{i})^{2} = ∑_{i=1}^{k} |f_{i}(ξ) − b_{i}|^{2} = ∑_{i=1}^{k} δ(ξ, H_{i})^{2}.

The distance given in the definition 2.1 has some convenient geometrical properties. Under obvious hypothesis, one directly proves: Proposition 2.5. Let a ∈ U and a C^{1}-function f : U → ℝ^{k}.

1) If A ∈ M_{k}(ℝ) is nonsingular and g = A · f, then δ(a, f^{−1}(0)) = δ(a, g^{−1}(0));

2) If Q ∈ M_{n}(ℝ) is an orthogonal matrix and T : ℝ^{n} → ℝ^{n}, T x = Qx + c is an isometry, then δ(a, f^{−1}(0)) = δ(T a, T f^{−1}(0));

3) For any λ > 0, δ(λa, g^{−1}(0)) = λδ(a, f^{−1}(0)), where g(x) = f (\frac{x}{λ}).

3 A least-square interpolation problem and an algorithm

Let A = {A_{1}, ..., A_{N}} ⊂ U (U ⊂ ℝ^{n}open) be a fixed set of distinct ”points of surveillance”. For any pattern Φ = f^{−1}(0), where f : U → ℝ^{k} is a C^{1}-function such that each A_{i} is regular for f, we put

(3)  δ(A, Φ) = ∑_{i=1}^{N} δ(A_{i}, f^{−1}(0))^{2}.

One can formulate the following problem:

(Π). Given A, determine a function f which minimizes δ(A, f^{−1}(0)). Such a problem could have some applications in Pattern Recognition [3]. In fact it is a geometrical nonlinear variant of the least squares method and the solution is generally not unique. In [1] one proves that whenever F ⊂ ℝ^{n} is closed and f : ℝ^{n} → ℝ is continuous, with F = f^{−1}(0), then for any continuous function ϵ : ℝ^{n} → (0, ∞) there exists a C^{∞}-function g : ℝ^{n} → ℝ such that g^{−1}(0) = F and ∀x ∈ ℝ^{n}, |f(x) − g(x)| < ϵ(x). This shows that for any pattern there are smooth functions which define it.

In what follows, we present the case when the unknown belongs to a parametrized family of functions (e.g. quadrics, cubics, spline-functions etc).

Let P ⊂ ℝ^{m} be an open subset in a parameters space and F : U × P → ℝ^{k} be a map such that for any p ∈ P, p = (p_{1}, ..., p_{m})^{T}, F determines a C^{1}-function f : U → ℝ^{k}, x → F (x, p); suppose there is a bijective correspondence between such functions and parameters (this happens for instance in the case of the

polynomial functions of degree at most d, where P = ℝ^{n} and m = ( \begin{matrix}n + d \\ n \end{matrix} );

in such a case, the function δ(A, Φ) becomes a function of the coefficients of the polynomials). In the case of the linear dependence on parameters, we also can impose supplementary relations on parameters, which do not modify the solution; for instance, the following condition

(4)  ∑_{i=1}^{N} J_{f}(A_{i}) · J_{f}(A_{i})^{T} = I_{k}.

Indeed, since the matrices J_{f}(A_{i}) · J_{f}(A_{i})^{T} are symmetrical, positively semidefinite and nonsingular (A_{i} being supposed regular points of f), these matrixes result positively definite and the same is true for the matrix C = ∑_{i=1}^{N} J_{f}(A_{i}) · J_{f}(A_{i})^{T}. So there is an orthogonal matrix Q ∈ M_{k}(ℝ) such that Q^{T} · C · Q = I_{k} and therefore Qf will verify (4).

Suppose now that the points A_{1}, ..., A_{N} are ”sufficiently near” of f^{−1}(0), in the sense that δ(A_{i}, f^{−1}(0)) ≃∥ f(A_{i}) ∥, for each i. By (3), δ(A, f^{−1}(0)) ≃ ∑_{i=1}^{N} ∥ f(A_{i}) ∥^{2}. Finally suppose fixed r C^{2}-functions φ_{1}, ..., φ_{r} : U → ℝ^{k}, linearly independent and consider that the unknown C^{2}-function f has the form f = (f_{1}, ..., f_{k})^{T}, where f_{i} = ∑_{j=1}^{r} p_{ij}φ_{j}, 1 ≤ i ≤ k, with p_{ij} unknown real constant (as parameters); in matricial writing, f = P^{T} · φ, where P = (p_{ij}) ∈ M_{r,k}(ℝ). Check the function f, that is the matrix P, such that the sum ∑_{i=1}^{N} ∥ f(A_{i}) ∥^{2} is minimum, with the restriction (4). In this case,

∑_{i=1}^{N} ∥ f(A_{i}) ∥^{2}= ∑_{i=1}^{N} f(A_{i})^{T} · f(A_{i}) =

= ∑_{i=1}^{N} φ(A_{i})^{T} · P · P^{T} · φ(A_{i}) = tr(P^{T} · A · P),

where A = ∑_{i=1}^{N} φ(A_{i}) · φ(A_{i})^{T} ∈ M_{r}(ℝ) is a known symmetrical, positively semi-definite matrix. On the other part, since

J_{f}(A_{i}) · J_{f}(A_{i})^{T} = P^{T} · J_{φ}(A_{i}) · J_{φ}(A_{i})^{T} · P,

by putting

B = ∑_{i=1}^{N} J_{φ}(A_{i}) · J_{φ}(A_{i})^{T} ,

with B ∈ M_{r}(ℝ) known, the relation (4) becomes P^{T} ·B·P = I_{k}. Thus, the above least-squares problem reduces to the following matricial one: check a matrix P ∈ M_{r,k}(R) of rank k such that P^{T} · B · P = I_{k} and tr(P^{T} · A · P) be minimum.

In order to solve this, we first remark that P^{T} ·A·P can be assumed diagonal [indeed, P^{T} · A · P is symmetrical and positively semi-definite, hence there is an orthogonal matrix Q ∈ M_{k}(ℝ), i.e., Q · Q^{T} = I_{k}, such that Q^{T} (P^{T} · A · P)Q is diagonal; but tr(P^{T} ·A·P) = tr((P Q)^{T} ·A·P ·Q) and moreover, P^{T} ·B ·P = I_{k} if and only if (P Q)^{T} ·B ·P ·Q = I_{k}. Thus, one can substitute P by P Q]. Denote by x_{1}, ..., x_{k} the column vectors of the matrix P; then tr(P^{T} ·A·P) = ∑_{i=1}^{k} x_{i}^{T} Ax_{i} and the relation P^{T} ·B·P = I_{k} can be written x_{i}^{T} ·B·x_{j} = δ_{ij} for any 1 ≤ i, j ≤ k. Consider the Lagrangean

L(x_{1}, ..., x_{k}, L) = ∑_{i=1}^{k} x_{i}^{T} Ax_{i} − ∑_{u=1}^{k} ∑_{v=1}^{k} λ_{uv}(x_{u}^{T} Bx_{v} − δ_{uv}),

where L = (λ_{uv}); 1 ≤ u, v ≤ k is the matrix of the Lagrange multipliers; the necessary (here also sufficient) extreme condition is ∇L = 0 and this leeds to AP − BP L = 0. But P^{T} AP = D (diagonal) hence L = I_{k}L = (P^{T} BP)L = P^{T} (BP L) = P^{T} AP = D. Therefore the matrix L will be diagonal and much more, positively semi-definite (since A is so). Put L = diag(α_{1}, ..., α_{k}) and the relation AP − BP L = 0 leeds to (A − α_{i}B)x_{i} = 0, 1 ≤ i ≤ k, hence x_{i} is an eigenvector of the matrix-bundle A − αB, with an eigenvalue α_{i}. One can apply different methods to determine the vectors x_{i} and, by this, the matrix P.

All the above can be shortly concentrated in the following. Proposition 3.1. (Algorithm to solve the problem Π). Let U ⊂ ℝ^{n} be an open set. Fix r linearly independent C^{2}-functions φ_{1}, ..., φ_{r} : U → ℝ^{k} and N ”points of surveillance” A_{1}, ..., A_{N} ∈ U. The problem is to look for a function f : U → ℝ^{k} [whose components f_{1}, ..., f_{k} are linear combinations of φ_{1}, ..., φ_{r}, namely (f_{1}, ..., f_{k}) = (φ_{1}, ..., φ_{r})P, with P ∈ M_{r,k}(ℝ) a matrix to be determined], such that the pattern f^{−1}(0) is the nearest to the set A = {A_{1}, .., A_{N}}.

Step I

Determine the r×r-matrices A = ∑_{i=1}^{N} φ(A_{i})·φ(A_{i})^{T} and B = ∑_{i=1}^{N} J_{φ}(A_{i})· J_{φ}(A_{i})^{T}, where φ = (φ_{1}, ..., φ_{r})^{T}.

Step II

The column-vectors x_{1}, ..., x_{k} of the looked for matrix P (of rank k) are just the eigen vectors of the matrix-bundle A−αB and moreover verify P^{T} BP = I_{k}. One determines thus x_{1}, ..., x_{k} and P.

As well as the matrices A, B are known (that depending on the choise of the functions φ_{1}, ..., φ_{r}), this algorithm requires O(r^{3}) operations.

4 A Newton-Raphson type result

In the paper [2] it was proved the following result:

Let D ⊂ ℝ^{n} be a nonempty convex bounded set, f : D → D a C^{2}-map such that any point of f^{−1}(0) is nonsingular for f. For u ∈ D fixed, define the function h_{u} : D × [0, 1] → ℝ^{n}, h_{u}(x, t) = f(x) + (t − 1)f(u). Then for a.e. u ∈ D, the set

{(x, t) ∈ D × [0, 1]|0 ≤ t ≤ 1, h_{u}(x, t) = 0}

either consists of a finite number of closed curves in D ×[0, 1], or a finite number of arcs in D × (0, 1) with their ends in D × {1} or D × {0}, or a finite number of arcs which start from D ×{0} and end in D ×{1}; all these three kinds of curves are disjoint, of class C^{1}. Therefore, one can found a solution for the equation f(x) = 0 by following the curve h_{u}^{−1}(0) which starts from (u, 0) for some u ∈ D; such a curve will attain a solution ξ of the equation f(x) = 0 as soon as t = 1 is touched.

Let f : U → ℝ^{n} be a C^{2}-map (U ⊂ ℝ^{n} open) and h(x, t), h : U × ℝ → ℝ^{n} a C^{2}-map such that for any x ∈ U, h(x, 1) = f(x) and the equation h(x, 0) = 0 has a solution u ∈ U; for instance, take h(x, t) = f(x) + (t − 1)f(u) as above, or h(x, t) = (1 − t)(x − u) + tf(x). Suppose that (γ) = h^{−1}(0) is a curve in ℝ^{n+1} which joins (u, 0) and a point (ξ, 1) such that f(ξ) = 0. Let a ∈ (γ) be fixed; choose a tangent vector τ in a at (γ) such that J_{h}(a) · τ = 0, ∥ τ ∥= 1

and det ( \begin{matrix}J_{h}(a) \\ τ^{T} \end{matrix} ) > 0. Then choose a step p > 0 sufficiently small such that if

b = a+pτ, then h(b) is near to 0. Put c = b−J_{h}(b)^{+} ·h(b). By the definition 2.1, δ(b, h^{−1}(0)) =∥ b − c ∥; c is just near the point where the hyperplane, passing by b and orthogonal to b, intersects (γ). If u is a point of simple bifurcation for k, then in the neighborhood of u, h^{−1}(0) represents the union of two curves γ_{1}, γ_{2}; if we take the arc length s as parameter and γ_{1}(0) = u, γ_{2}(0) = u, then

det ( \begin{matrix}J_{h}(s) \\ γ̇(s)^{T} \end{matrix} ) changes its sign in s = 0 for γ = γ_{1} and γ = γ_{2} and conversely;

both the curves are suitable for the next algorithm. Proposition 4.1. Suppose f, h satisfy the conditions of the above formulated result of [2]. Take a solution u ∈ D of the equation h(x, 0) = 0. Take a_{0} = (u, 0) and apply the described scheme which yields the sequence (a_{k})_{k≥0}, where a_{k+1} = c and a_{k} = a as above. This sequence converges to a point (ξ, 1) such that f(ξ) = 0.

The proof uses a typical reasoning for the Newton-Raphson method. The new thing is that b can be singular for h and for this reason we use the pseudoinverse; in fact by the Sard theorem almost all points of ℝ^{n} are regular values for the map h.

References

[1] F.Broglia, A.Tognoli, Approximation of C^{∞} functions without changing their zero-set, Annales Inst.Fourier, 39, fasc.3, 611-632, 1989.

[2] S.N.Chow, J.Mallet-Paret, J.A.Yorke, Finding zeros of maps by homotopy methods, Math.Programming (A. Bachem ed), Springer Verlag, 15-56, 1983.

[3] V.Neagoe, O.Stănăşilă, Teoria recunoaşterii formelor, Ed.Academiei Romane, 1992.

[4] R.Penrose, A generalized inverse for matrices, Proc. Cambridge Phil.Soc., 51, 406-413, 1955.

[5] G.Taubin, Nonplanar curve and surface estimation in 3-space, Proc. IEEE Conf.Roboties Automation, 1988.

University Politehnica of Bucharest Department of Mathematics II Splaiul Independentei 313 77206 Bucharest, Romania