Detecting and approximating decision boundaries in low dimensional spaces

A method for detecting and approximating fault lines or surfaces, respectively, or decision curves in two and three dimensions with guaranteed accuracy is presented. Reformulated as a classification problem, our method starts from a set of scattered points along with the corresponding classification algorithm to construct a representation of a decision curve by points with prescribed maximal distance to the true decision curve. Hereby, our algorithm ensures that the representing point set covers the decision curve in its entire extent and features local refinement based on the geometric properties of the decision curve. We demonstrate applications of our method to problems related to the detection of faults, to Multi-Criteria Decision Aid and, in combination with Kirsch's factorization method, to solving an inverse acoustic scattering problem. In all applications we considered in this work, our method requires significantly less pointwise classifications than previously employed algorithms.


Introduction
Let us consider a piecewise constant function f : R m ⊃ Ω → {1, 2, . . ., n} with Ω being compact, simply connected and equipped with a piecewise smooth boundary, m ∈ {2, 3}.Such f subdivides Ω into mutually disjoint subsets Ω i := f −1 (i) with Ω = ∪ n i=1 Ω i and Ω i ∩ Ω j = ∅, if i = j.We assume that each Ω i features a piecewise smooth boundary.We are interested in approximating Γ i,j := Ω i ∩ Ω j relying on as few evaluations of f as possible and present in this work an algorithm for this task.We choose this quantity as a measure of efficiency because evaluating f can be arbitrarily costly in applications and dominates the runtime in such a case.Our problem can be immediately understood as a classification problem, such that we identify Ω i with a class i and interpret the curves or surfaces of discontinuity Γ i,j as decision curves or surfaces.
One field of application is economics and operations research.Multicriteria Decision Aid (MCDA) methods can help a decision maker to choose the best one from a finite number of alternatives based on different, even conflicting criteria.MCDA methods assume that a decision depends on quantifiable parameters (x 1 , . . ., x m ) ("input factors") and is drawn deterministically.For an overview over various MCDA approaches, applications and case studies, we refer to [14,25] among many others and the references cited therein.In this context, Ω i is the set of all input factors that lead to the decision for the i-th alternative in the MCDA method.Analysing the decision process with respect to the input factors means consistently describing all Ω i based upon evaluating the MCDA method for arbitrary combinations of input factors.This can be achieved by approximating all Γ i,j .Computing a reconstruction of an obstacle in three dimensions is a field of application in acoustic scattering theory.More precisely, one wants to determine the support of an inhomogeneous object (its boundary to be exact) from measured far-field data which typically is a desired task in non-destructive testing.The far-field data are obtained for different incident waves and measured points on the unit sphere.Several reconstruction algorithms to find the boundary of the unknown inhomogeneity are available such as iterative methods [22], decomposition methods [31] and sampling/probe methods (see [26] for a detailed overview).The latter ones can be further categorised into the linear sampling method, the generalized linear sampling method, the factorization method, the probe method and variants of it (refer to [12,6,20,18,26], respectively).However, here we will focus on the classical factorization method for the acoustic transmission problem, refer also to [4], with which one can decide if a given point is located inside or outside the obstacle.Therefore, the factorization method transfers the reconstruction of an obstacle to a classification problem and thus into a field of application of our method.Note that the far-field data within [4] has also been used in [8] and [17].
Finding and approximating the sets Γ i,j , sometimes called fault lines, is important in the exploration of natural resources.The presence and the location of Γ i,j can provide useful insights for exploring and later on exploiting of oil reservoirs [16] and play a significant role in some geophysical applications [15].The underlying mathematical problem is closely related to ours, albeit not the same, as usually, the function f considered does not provide integer values as in our case.Therefore, an additional algorithm for detecting fault lines is required then, and the classification of a single point may be not trivial anymore as it is in our case.Moreover, algorithms for fault detection may need to deal with noisy data (e.g.[10]), whereas we consider certain data only.
Many algorithms for detecting and approximating the sets Γ i,j have been proposed, like [5,16,15,1,10] among many others, which all feature strengths and weaknesses.However, the vast majority of these approaches restrict to the 2D case, whereas we present a method for 2D and 3D.The algorithm proposed in [1] and the work cited therein was the starting point for our research.Classification is one of the standard problems in Machine Learning.There are a lot of powerful and versatile algorithms available which could readily applied to our problem; we refer to [9] for an overview.These methods are however designed for uncertain data in high dimensional spaces, whereas we consider secure data in low dimensional spaces, a completely different use case.
Our method approximately describes the Γ i,j by providing a set of points with a guaranteed maximal normal distance to Γ i,j .These points are intended for constructing a polygon (2D) or a surface triangulation in 3D.While there are more sophisticated and elegant ways of describing these sets, it allows us to (approximately) replace an actual classification by a simple and fast point-in-polygon test.
This article is organised as follows: We describe our algorithm for 2D and 3D in Section 2 and elaborate on the direct and inverse acoustic scattering problem, one of our applications, in Section 3. We present results and applications of our algorithm to the detection of faults, to decision modeling and to inverse scattering in Section 4 and finally conclude in Section 5.In this section, we present our algorithm for approximating Γ i,j = ∅ for m = 2 by sufficiently many well distributed and ordered points sufficiently close Γ i,j .This implicitly constitutes a polygonal description of Γ i,j .For convenience, we provide flow charts (Figs. 1 and 2) and describe selected building blocks in detail.In the flow charts, actual building blocks are typed in monospace lettering.Algorithm 2.1 (initialise).This algorithm aims at providing initial approximations to any Γ i,j ; we refer to Fig. 1 (right) for an overview.From now on, we assume Γ i,j = ∅.In Building block initialset, we sample f on Ω rather coarsely and obtain an initial point set X along with the corresponding classification information.Building block barymeans creates additional sampling points in the vicinity of any Γ i,j .Following Allasia et al. [1], we employ a k near -nearest neighbour approach: For any x ∈ X, let N (x) be the set of the k near -nearest neighbours of x in We compute the barycentres b c i of all N c i (x), 1 ≤ i ≤ r, and then their arithmetic means y c i ,c j = 0.5(b c i +b c j ), 1 ≤ i < j ≤ r.Let M (x) be the set of all y c i ,c j generated from N (x).If N (x) ⊂ Ω for some , we set M (x) = ∅.We end up with M = x∈X M (x).By definition of M, there are no duplicate points; however, a practical implementation requires removing duplicates.After classifying the points in M, we repeat barymeans on M obtaining sets M 2 (x) for any x ∈ M, then M 2 = x∈M M 2 (x) and ultimately a further enriched set of sampling points X = X ∪ M ∪ M 2 .Building block iniapprox computes initial approximations for all Γ i,j .For any x ∈ (M ∪ M 2 ) ∩ Ω i , we search the nearest point x ∈ X ∩ Ω j , j > i.We use x and x as starting points for a bisection algorithm on the line x x.If the bisection algorithm is successful, we end up with a point pair x (i) ∈ Ω i and x (j) ∈ Ω j with x (i) − x (j) ≤ 2ε b , where ε b is a user-prescribed threshold.Then, the distance of x (i,j) = 0.5(x (i) + x (j) ) to Γ i,j is at most ε b .We subsume the bisection process up to ε b and computing x (i,j) from its results in Building block bisection.From now on, we consider point triplets for approximating Γ i,j only; for any such triplet x, the superscript (i) denotes the point in Ω i , the superscript (j) its counterpart in Ω j and the superscript (i, j) the arithmetic mean of the two points.We end up with a set of triplets S i,j .We moreover set may occur that some triplets in S i,j are tightly clustered.We thin such clusters as they add to complexity but not to accuracy by removing appropriate triplets (Fig. 3).After cluster removal, initialise provides sets of triplets S i,j .
Remark 2.2.Building block initialise detects Γ i,j = ∅ by S i,j = ∅.However, depending on X, it may happen that S i,j = ∅ even if Γ i,j = ∅.Reliably detecting all non-empty Γ i,j depends on a sufficiently large X and thus ultimately on the user.
Remark 2.5.The bisection-based approach in Building block iniapprox implicitly assumes that x x intersects Γ i,j and therefore may fail if this does not hold (Fig. 5).However, such failure modes occur only rarely in practical computations, and we implemented fallbacks in that case.
Algorithm fill (for an overview, we refer to the flow chart in Fig. 2) provides triplets with a maximal user-prescribed distance ε b to Γ i,j and maximal mutual distance ε gap based upon      6: Sorting due to Allasia may fail (left; sorting indicated in blue), whereas our method succeeds in the present situation (right).We obtain two ordered subsets (direction of sorting is indicated by arrows which after combination yields the correct ordering.The dashed line represents a connection to the nearest neighbour rejected due to angle in our approach.S i,j , i.e. the result of Algorithm 2.1 (initialise).As we assume that ε b ε gap , we define the distance of two triplets x, y ∈ S i,j as x (i) − y (i) and the distance of x to Γ i,j as the distance of x (i,j) to Γ i,j .Following a bottom-up approach, we start with discussing selected Building blocks employed in fill.
Building block 2.6 (sort).This Building block is to sort a set of triplets S according to their position along Γ.We omit the indices i and j for clarity.Let us assume that Γ is piecewise smooth and fulfills an inner cone condition with angle β angle .We first search a triplet x start ∈ S closest to the boundary of Ω and assume x start to be the first triplet in the sorted set.We initially set Š = {x start }.Let now the triplets in Š ⊂ S be already sorted with x r being the last of those, r > 1.We consider the k sort nearest triplets y 1 , . . ., y ksort to x r in S \ Š, sorted by increasing distance to x r .If for s = 1, ∠(x r − x r−1 , x i − y s ) < β angle , we set x r+1 = y s and add x r+1 to Š; otherwise, we repeat with s + 1.If s > k sort , we store Š, set S = S \ Š, and repeat the sorting procedure until S = ∅ ending up with a finite number of disjoint ordered subsets.At the end, we combine all sorted subsets based upon the Euclidean distance between first and last points of the subsets reversing the order of a subset if necessary.
Remark 2.7.Allasia et al. [1] present a simpler sorting method than ours as they do not enforce the condition ∠(x r − x r−1 , x i − y s ) < β angle .However, it may fail if x start is not the true starting point and if additionally the points are unevenly distributed along Γ (Fig. 6) in contrast to ours.Of course, our sorting method can fail as well, but due to our experience, it is more reliable than Allasia's method and works sufficiently well.There are many more sophisticated sorting methods based e.g.upon graph theory [2,3,13,24], which are more reliable than our approach, but more time-consuming and much harder to implement.
For describing Building block inipairs, which is part of fill, we introduce another two Building blocks, which will be used in adapt as well.
Building block 2.8 (estcurv).Let be S loc = {x 1 , . . ., x r }, r > 2, a set of ordered points near Γ up to ε b .This Building block estimates the curvature c of Γ in x , 1 < < r by leastsquares fitting an approximation using Gaussian radial basis functions (RBFs) and then c by the curvature of that approximation in x .We hereby assume that after shifting and suitable rotation, Γ can be locally represented as a graph of an unknown function.As the points in X loc are located on Γ only up to ε b , we penalise the second derivative of the RBF approximation subject to a maximal residual of ε b .This coincides with the maximal deviation in the value at an approximation point.We employ Tikhonov regularization with parameter estimation using Morozov's discrepancy principle.
If Γ cannot be considered a graph of a function even after rotation, we draw as a fallback a circle through the points with indices − 1, and + 1 and use the inverse of its radius for estimating c .We estimate c in the first or last point of S loc by drawing a circle through the first or last three points in S loc .
Building block 2.9 (esterror).This Building block estimates the maximal deviation δ of a smooth curve from a straight line between two points on the curve with distance d.A straightforward calculation reveals that where c denotes the maximal curvature of the curve between the two points.For S loc as in Building block 2.8, we estimate the maximal deviation δ of Γ from the straight line between consecutive points x and x +1 by replacing c with max{c , c +1 } from Building block 2.8 (estcurv).Hereby, we rely on (1) and neglect higher order terms.
Now we are prepared to discuss fill.
Algorithm 2.10 (fill).Building block 2.6 (sort) sorts all triplets in S i,j according to their position near Γ i,j .We detect gaps in the representation of Γ i,j by S i,j by considering subsequent triplets x and x +1 .If the distance d of x to x +1 is larger than a user-prescribed threshold ε gap , we consider this a gap and aim to equidistantly add R = d /γ triplets near Γ i,j between x and x +1 .To do so, Building block inipairs places new points z ,r , +1 and computes from these initial point pairs Here, n denotes the (estimated) outer normal unit vector of Ω i near x ,r .Applying Building block 2.9 (esterror) to with user-prescribed safety factors ε safemax and ε safemin .
Remark 2.11.Having a local RBF approximation of Γ at hand when computing c in estcurv, it seems to be straightforward for efficiency reasons to choose points z ,r on that RBF approximation instead of just subdividing a straight line.Numerical experiments did not show any significant advantage of that approach compared to ours for fill.This could be related to the uneven distribution of points on Γ near gaps to fill, which may decrease the quality of approximation.Therefore, we stick to the easier approach presented here.However, estimating the curvature using Building block 2.8 (estcurv) is sufficiently reliable for efficiently computing starting pairs.
However, the pairs of starting points (aka starting pairs) obtained from inipairs are not necessarily valid.We call a starting pair for approximating Γ i,j valid, if one of its points belongs to Ω i and the other one to Ω j .Therefore, we introduce Building block 2.12 (startpairs).
Building block 2.12 (startpairs).This Building block obtains a valid starting pair from a pair of points (x + ,r , x − ,r ).If the starting pair is already valid, we return it as the result.If x + ,r or x − ,r belongs to a third class, we stop without result.If both points belong to the same class, we reflect x + ,r on z ,r obtaining x ,r (Fig. 7).If both (x + ,r , x ,r ) and (x − ,r , x ,r ) still belong to the same class, we repeat this process with changing roles and escalating distances at most k rep times (typically k rep = 3) and stop, if any of the resulting point pairs is valid or one of the points belongs to a third class.
) is not, as both points belong to the same class.However, (x − ,2 , x ,2 ) is valid.1.We iterate the process of filling gaps in fill, as the arclength of Γ between two subsequent triplets may considerably exceed the length of the straight line between them, such that after a first pass, the mutual distance of subsequent triplets may still exceed ε gap in some cases.
2. If only two or even less triplets are known on Γ i,j , estimating curvature is impossible with estcurv, and computing α by (2) fails.For this case, we implemented fallbacks.
Using bisection, fill creates new triplets, yielding new sets S i,j .
Remark 2.15.If Ω i is not simply connected, some Γ i,j may consist of several components.We detect this using fill.If there are some significant gaps which can not be filled, it indicates the presence of several components.We then subdivide S i,j correspondingly and proceed with every subset separately.
Fig. 8 indicates that even with filling gaps, S i,j may not appropriately represent Γ i,j , as parts of Γ i,j before the first known triplet x 1 and after the last known may be neglected.Algorithm expand is used to expand S i,j to a representation of the complete curve Γ i,j .It relies on several building blocks, which we discuss first.
Building block 2.16 (extrapolate).For a given ordered set S with n ≥ 2 triplets close up to ε b to a curve Γ and with average distance d avg , we fit a polynomial with degree n − 1 in local coordinates.We compute these coordinates by least-squares-fitting a line to S. As points in S are located on Γ up to ε b only, we do not interpolate, but penalise the second derivative in a least-squares approximation.Following Morozov's discrepancy principle, we regularise such that the maximal residual is approximately ε b .
Building block 2.17 (stepsize).Provided that Γ extends sufficiently far before x 1 ∈ S, it seems to be straightforward to seek for a new triplet with distance ε gap to x 1 .However, we limit the step size for extrapolation based upon the curvature c of Γ in the vicinity of x 1 .Extrapolating far is unreliable in case of large curvature c, and adapt will insert additional points afterwards in that region anyway for accuracy reasons, such that it is much more efficient to adjust the step length to the local properties of Γ beforehand.Let us assume that a polygonal final approximation of Γ may deviate at most by ε err from Γ.We compute the step size which would lead to a deviation of ε err from a straight line segment between x 1 and the new triplet yet to compute.This is a natural upper bound for the step size l extra in extrapolation.
Rearranging (1) and neglecting higher order terms leads to for the maximal admissible step length l max .However, evaluating (3) is numerically unstable if cε err is small.We set (cε err Algorithm 2.18 (expand).This algorithm aims at finding triplets near Γ i,j beyond the first or last known in S i,j until the start or end of Γ i,j is reached or Γ i,j turns out to be a closed curve.
For the sake of simplicity, we refer in what follows to finding triplets before the first one in S i,j .Finding triplets near Γ i,j beyond the last one in S i,j works analogously.We add a new triplet before x 1 by extrapolating an approximation γ of Γ i,j with Building block 2.16 (extrapolate) , . . ., x (i,j) kextra with x ∈ S i,j and choose the step size according to Building block 2.17 (stepsize).This way, we obtain an extrapolating curve γ and some s 0 such that x 1 − γ(s 0 ) ≈ l extra .We then create a point pair (x + s 0 , x − s 0 ) based on γ(s 0 ) in a similar way as in fill.Computing a valid starting pair with Building block 2.12 (startpairs) and subsequent bisection yields a new triplet in S i,j .We repeat this process until we reach the true starting point of Γ i,j .As heuristic criterion for γ(s 0 ) exceeding this starting point, we consider (Fig. 9).In this case we employ Building block 2.19 (reducestepsize) for obtaining a valid pair of points (x + s , x − s ), which represents the starting point of Γ i,j .After adding it to S i,j , expand (displayed as blue dots), we construct γ by extrapolate.As (x + s 0 , x − s 0 ), displayed as cyan-coloured dots, fulfills (4), we apply bisection with respect to the parameter s until a valid starting pair based upon γ(s) can be constructed (displayed in magenta), from which we compute the final triplet in S 1,2 by bisection.terminates.If Γ ij is closed, expanding S i,j as described above would lead to an endless loop.Therefore, we start expanding S i,j , but detect after every addition of a triplet, if Γ i,j is closed.If so, we stop expanding and resort.We skip the details to keep the presentation uncluttered.Algorithm expand yields approximating sets Si,j .
Building block 2.19 (reducestepsize).Using the parameter value s 1 corresponding to x (i,j) 1 as lower bound and s 0 as upper bound, we obtain s, which leads to a valid starting pair (x + s , x − s ) based upon γ(s) and fulfils γ(s) − γ(s ) < ε b by bisection with respect to condition (4).Here, s denotes the second to last parameter in the bisection process (compare Fig. 9).Algorithm 2.20 (adapt).Based on esterror, this algorithm inserts a triplet (approximately) halfway between consecutive triplets x and x +1 , if esterror indicates an error larger than ε err and removes a triplet, if the estimated error of both line segments the triplet belongs to is smaller than ε coarse .In contrast to fill, we employ the local RBF approximation from estcurv for computing an initial point x new between x and x +1 when refining.With x ± new = x new ± α n, we then proceed as in fill.We compute α according to (2), but replace δ by δ = 1/16c 3 d 4 , as the error due to (1) refers to an approximation by line segments and is overly pessimistic for an RBF approximation.For robustness, we never delete consecutive triplets in one pass of the adaptive loop.After at most k adap refinement and coarsening sweeps, we end up with final sets Ŝi,j .

Approximating Γ i,j in 3D
For approximating Γ i,j in three dimensions, we stick to the general procedure for approximating these sets in two dimensions (Fig. 1).While initialise remains unchanged, fill, expand and adapt differ from their 2D counterparts, as these exploit ordering of the points on a decision curve or fault line.There is however no straightforward ordering of points on a surface.Because fill and adapt rely on esterror as in two dimensions, we discuss error estimation first.
Building block 2.22 (esterror).This Building block aims at estimating the maximal error e of a linear triangulation of Γ based upon a finite set of points S near Γ up to ε b .Let T be a Delaunay triangulation of S loc ⊂ S and let us assume that Γ can be represented on the support of T by an unknown function g after appropriate change of coordinates.For some triangle T ∈ T , it holds according to [29, Theorem 4.1] Here, R describes the radius of the circumcircle of T , d the distance from its center to T and I T g the Lagrange interpolant of g on T ; with |g| 2,∞,T , we denote the L ∞ -norm of the second derivative of g on T .It remains to estimate |g| 2,∞,T .As g is unknown, we employ an RBF approximation ϕ as in Building block 2.8 (estcurv) instead and approximate |g| 2,∞,T by evaluating its second derivative in the vertices of T and its center.The maximum of these four values yields the desired approximation φ of |g| 2,∞,T and therefore Algorithm 2.23 (fill).For any triplet x in S i,j , we search the k near nearest neighbours x 1 , . . ., x knear and switch to a local 2D coordinate system by computing the optimal fitted plane in the sense that the sum of the squared distances between the points and the plane is minimal (see [27]).We then compute in local coordinates a Delaunay triangulation of the patch S loc = {x, x 1 , . . ., x knear }.If the maximal edge length of a triangle in this triangulation exceeds ε gap , we mark the centre of that triangle as starting point for constructing a valid starting pair similar as in Building inipairs, as long as the triangle is not too anisotropic, i.e. does not contain angles approaching 0 • or 180 • .However, large gaps in S i,j are not covered by this local approach (Fig. 12).If there is a gap in S i,j and if the current point is at the boundary of this gap, it will be at the boundary of the current patch.We detect this by exploiting the local coordinate system.We assume that a patch in the vicinity of a large gap is at least slightly elongated in tangential direction to the patch boundary.Therefore, we consider the first local coordinate.If the first local coordinate of the current point is almost the minimum or maximum of all respective first coordinates of the patch, then it is at its boundary, and we assume a large gap to fill.In this case, we construct x new on the elongated line between the centre of gravity of the patch and x, enrich S local by x and continue as described above.Looping over all triplets in S i,j , removing duplicate starting points and starting points which are extremely close enables us to compute triplets using Building block bisection similar as in Algorithm 2.10 (fill).This yields an enriched representation of Γ i,j .Repeating this procedure with the enriched set until no gaps are detected anymore leads to S i,j .
Algorithm 2.25 (expand).When expanding a representation S i,j of Γ i,j , we distinguish between expansion towards inner boundaries Γ i,j ∩ Γ k, and expansion towards outer boundaries Γ i,j ∩ ∂Ω.We consider expanding to inner boundaries first.Let S I i,j ⊆ S i,j be all triplets closer than 0.75ε gap to a triplet in another S k, .For some x ∈ S I i,j , let x * be the triplet in S k, closest to x.We estimate the normal vector n of Γ i,j in x and project xx * on the corresponding tangential plane at x, yielding t.Let x ∈ S i,j be the nearest triplet to x fulfilling ∠(−t, xx ) < α expand and E the plane which contains x and is spanned by xx and n.We reduce the expansion to the two-dimensional case by expanding the curve Γ i,j ∩ E with Algorithm 2.18 (expand) from Section 2.1, using x and x as initial set of triplets (compare Fig. 12, right).
For expanding to outer boundaries, we construct S B i,j as follows: For each coordinate direction i, we select n expand,i points in S i,j which are minimal or maximal with respect to where n outer stands for the outer normal vector of the assigned boundary facet.Figure 15 illustrates the purpose of condition (6).For any triplet in S B i,j , we proceed as for expanding to inner boundaries, but with t = n outer .All new triplets in the same facet F of ∂Ω constitute the set S F i,j .As this set represents the curve Γ i,j ∩ F in the plane F , we are confronted with approximating a decision curve or fault line in two dimensions.Therefore, we apply Algorithms 2.10 (fill) and 2.18 (expand) to S F i,j .After adding the new triplets on the boundary of Γ, we apply Algorithm 2.23 again and end up with an enlarged set Si,j .
Algorithm 2.26 (adapt).In contrast to the two-dimensional case, we do not adaptively coarse the set of triplets.While this is possible and does not harm accuracy, it complicates computing a surface mesh from the set of triplets representing Γ i,j .
For adaptive refinement, we do not rely on a global triangulation of Si,j .Instead, for a given triplet x ∈ Si,j , let Sloc consist of the k near nearest triplets in Si,j to x.We least-squares fit a plane E to Sloc i,j as in Algorithm 2.23 (fill) and create a Delaunay triangulation T of Sloc i,j projected to E. For each non-degenerated triangle in T , we estimate the error applying esterror to S(i,j) loc .If the estimated error exceeds ε err , we employ the center of the triangle as a starting point for adding a new triplet.Looping over all x ∈ Si,j naturally leads to many duplicate or very close starting points which need to eliminated before adding triples via bisection.However, these duplicates do not harm the efficiency of our method, as no function evaluations are required before computing triplets from starting points.We enrich Si,j with the new triples created and repeat the refinement procedure at most k adap times or until no more starting points for computing triplets have been created.

Direct and inverse acoustic scattering problem
In this section, we briefly explain the direct and inverse scattering for the acoustic transmission problem.Most of the description can be found in [4] and we refer the reader to [4] for more details.
Let the scatterer D ⊂ R 3 be a given bounded domain with boundary ∂D assumed to be of class C 2,α .The normal unit vector ν points into the exterior E = R 3 \D of the scatterer.The exterior E assumed to be simply-connected is an infinite homogeneous isotropic non-absorbing acoustic medium which is characterised by the mass density e , mean compressibility κ e , and sound speed c e = 1/ √ κ e e .Likewise, the interior of D is characterised by i , κ i , and The given scatterer is excited by a time-harmonic acoustic incident plane wave of the form where k e = ω/c e is the wave number of the acoustic wave in the host medium, ω > 0 the angular frequency, and d ∈ S 2 the direction of incidence with S 2 = {x ∈ R 3 : x = 1} the unit sphere, where • denotes the standard Euclidean norm in R 3 .Note that the incident field also depends on k e , but this dependence is suppressed.The incident wave (7) interferes with the penetrable scatterer and creates two waves.The first wave is the scattered field u sca (x; d) defined for x ∈ E propagating outward and the second wave is the transmitted field u int (x; d) defined for x ∈ D. The total field in E denoted by u ext (x; d) is the superposition of u int (x; d) and u sca (x; d) each of which satisfies the Helmholtz equation (the reduced wave equation) in E with wave number k e .Likewise, the transmitted field satisfies the Helmholtz equation in D with wave number k i .Precisely, we have Due to the continuity of the acoustic pressure and the normal component of the particle velocity across ∂D yields the transmission boundary conditions where τ = i / e > 0 is the mass density ratio of the two media.To ensure a well-posed boundary value problem, the scattered field u sca (x; d) needs to satisfy the Sommerfeld radiation condition lim r→∞ r ∂ r u sca (x; d) − ik e u sca (x; d) = 0 with r = x .The classical acoustic transmission problem reads: find the functions

The direct acoustic transmission problem
Given the incident field (i.e. the direction of incidence d and the wave number k e ), the scatterer D (hence also its boundary ∂D), the wave number k i , and the parameter τ , one has to solve (8) -( 12) for u int (x; d) and u sca (x; d).In the direct problem, one is only interested in the far-field u ∞ ( x; d) of u sca (x; d) which is given by uniformly in all directions x ∈ S 2 .The far-field can be found by evaluating an integral equation over ∂D given two density functions determined by first solving a 2 × 2 system of integral equation of the second kind over ∂D (see [4,Section 4.1]).Of course, the integral equations at hand cannot be solved analytically and have to be solved numerically for example by the boundary element collocation method (see [21,Chapter 5] for more details).
To sum up, in the direct acoustic transmission problem one is interested in u ∞ ( x; d) for x ∈ S 2 given the scatterer's boundary ∂D and the direction of incidence d ∈ S 2 .The parameters k e , k i , and τ are given.

The inverse acoustic transmission problem
The parameters k e , k i , and τ are given.In the inverse acoustic transmission problem one tries to find/reconstruct the domain's boundary from the knowledge of the far-field patterns u ∞ ( x; d) for all x, d ∈ S 2 .This can be achieved with the factorization method originally invented by Kirsch (see [20]).The theoretical justification of the factorization method for the acoustic transmission problem is given in [4, Chapter 3] and shown to work practically in [4,Chapter 4].We briefly outline the algorithm: Assume that the far-field data are given for xi and dj with i, j ∈ {1, . . ., m} stored in the matrix A ∈ C m×m .First, compute a singular decomposition of A = U ΛV * with Λ = diag(λ 1 , . . ., λ m ).For a given point z ∈ R 3 compute the expansion coefficient of with respect to V by which is a matrix-vector multiplication (z) = V r z .Finally, we compute and plot the isosurfaces of z → W (z). The values of W (z) should be much smaller for z / ∈ D than those lying within D. The threshold value can be approximately determined by a scatterer such as the unit sphere and then reused for other scatterers as well.Note that until now, an equidistant set of points N within a predefined box have been used to find the values of W (z) leading to an amount of N × N × N function evaluations for such a "sampling" method (see also [8] for other sampling methods).This can be considerably reduced as shown in the next section.

Test cases related to the detection of faults
In this section, we consider test cases which have been defined in the context of the detection of fault lines or fault surfaces.
We compute all three test problems starting with X from Example 2.4 and the same parameters used there and in the subsequent examples in Section 2.1.Our results (Figs. 18, 19, and 20) demonstrate a successful approximation of all Γ i,j .We provide details in Tables 2 and 3 combined with the results for Test problem 2.3 from the preceding section.Algorithm iniapprox is the most demanding Building block in terms of function evaluations in all examples considered.Averaged over Test problems 2.3, 4.2, and 4.3, it takes on average 4.0 classifications for adding one triplet in fill.However, it takes 12.2 classifications per triplet added for Test problem 4.1.This is due to the fact that Γ 1,3 consists of two components, which is detected by a failed attempt of filling the large gap between these two components (Remark 2.15).This process adds classifications, but no triplets.In average over the four examples, expand requires 5.5 classifications per triplet added.A significant part of the classifications is required for the very first and very last step, as finding them requires bisection on the extrapolation curve γ, which involves at least one classification per iteration step.Considering the same ratio for adapt would be misleading, as triplets are added and removed.
Allasia et al. present numbers of function evaluations for their method of surface reconstruction applied to these test problems in [1].They report 5185 evaluations of f for Test problem 4.1 (our method: 2213), 3479 for Test Problem 4.2 (our method: 871), and 2099 for Test problem 4.3 (our method: 450).It turns out that the number of function evaluations is significantly higher than for our algorithm.However, this comparison is limited by several factors.In contrast to us,

Decision analysis and modeling
MCDA methods depend on several parameters called "input factors".Almost all common MCDA models require to set up a performance matrix P = (p i,j ), where p i,j encodes the objective benefit or cost of the i-th alternative with respect to the j-th criterion.Both finding suitable criteria for modeling the decision process and setting up P requires expertise (e.g.[28] among many others) and is beyond the scope of this work.Instead, we assume that P is given and exactly known, although this assumption may be questioned in many practical applications.
In addition to P , the decision maker needs to provide non-negative weights w = (w 1 , . . .w c ) which reflect the importance of the criteria from his point of view.Although originally intended as a decision support tool, MCDA methods have been used for some time for decision analysis, e.g. to predict decisions of actors under changed framework conditions, modelled by a changed performance matrix.In the context of decision analysis, w is in many practical applications not exactly known and hard to obtain, as e.g.surveys are time-consuming and prone to bias due to socially accepted answers and other effects.Moreover, the restriction to fixed weights ignores possible diversity within actors.All of the above motivates a robustness analysis of the decision predicted, which in many cases focuses on examining the robustness of the decision with respect to perturbations in w.
In what follows, we apply our algorithm for computing Γ i,j to such a robustness analysis and consider as example an application from the scientific monitoring of the mobility turnaround in Germany.Ball et al. [7] investigate car users' attitudes towards the purchase of hybrid (HEV) and electric vehicles (BEV) versus conventional cars with internal combustion engine (ICE).For this purpose, they identify 13 criteria and divide them into 5 categories ("Ecological","Economic", "Social", "Comfort", "Other").The authors weight both the criteria within the categories and the categories themselves.Taking the former weights as given, we obtain the performance matrix in Table 4 from the data in [7] [7], and weightings from the car users' point of view prior to normalisation (last row).
Weighting) [11] as in [7], we use the more complex MCDA method SIR-TOPSIS [30] here, which includes the very widespread MCDA methods Promethee II [19] and SAW as special cases.We omit all details on how TOPSIS works for the sake of brevity and refer instead to [30].SIR-TOPSIS requires as many other MCDA methods that the non-negative weights are normalized, i.e. n i=1 w i = 1.Therefore, the set of normalized admissible weights is the standard simplex in R n .For visualisation, we consider m = 3 or m = 4 of the weights to be variable, and the rest to be fixed.Let be w = w v + w f , where w v consists of the variable weights and zeroes elsewhere, and w f of the fixed ones, correspondingly.Therefore, normalisation implies w i,v = 1 − w i,f := c f such that the set of variable weights corresponds to a downscaled standard simplex in R m , which we embed in R m−1 by appropriate translation and rotation.This yields the equilateral triangle (for m = 3) and the regular tetrahedron (for m = 4) shown in Fig. 21.
For a 2D-visualisation of the decision space, we consider the weights corresponding to "Ecological", "Economic" and "Comfort" variable (Fig. 21) and colour all weights leading to a decision in favour of ICEs black, for HEVs blue and for BEVs green.We divide the weights w proposed in [7] as a representation of car user's mindset in 2020 in a fixed and variable part, setting w = wv + wf and display wv as bright red dot in Fig. 21.In contrast to [7], SIR-TOPSIS seems to predict a shift to HEVs under today's conditions.For measuring robustness, we consider the largest sphere around wv still fully contained in the set of weightings leading to HEVs W HEV and propose its radius ρ as simple measure of robustness.As we have a polygonal approximation of W HEV at hand, iteratively approximating ρ boils down to an intersection test of polygons, if we approximate the circle by a sufficiently fine polygon.We end up with ρ ≈ 0.06, which reconciles the findings of Ball et al. and ours: As ρ is rather small, the decision in favour of HEVs is not very robust, as even small perturbations of the weightings may lead to a decision in favour of ICEs.In [7], the decision towards ICEs was found not be very robust.However, both ours and Ball et al.'s results indicate that a broad shift towards BEVs is unlikely to happen under current circumstances.
For m = 3 and m = 4, we have c f = 0.9942.The downscaled standard simplex is rotated and translated to the equilateral triangle with vertices c f (0.4082, −0.7071) , c f (0.4082, 0.7071), and c f (−0.8165, 0).Therefore, we set Ω = c f [−0.9, 0.4082] × c f [−0.9, 0.9] and employ an initial point set X consisting of 100 Halton-distributed points, where points far away from the triangle have been discarded, as they cannot aid approximating Γ i,j (Fig. 22, left).We choose the values of all algorithm-related parameters as for the examples in Sections 2.1 and 4.1, except of ε gap = c f • 0.05.For m = 4, we start with X consisting of 500 Halton-distributed points in the vicinity of the tetrahedron and take all values for the parameters of the algorithm from Section 2.2, except of k adap = 3.For the final sets Ŝi,j , we refer to Fig. 22 (right and middle) and display the number of function evaluations in Table 5.These sets have been used for approximating ρ; Fig. 21 was generated with the proposed algorithm, albeit using a finer initial point set and modified parameter settings, as analytical descriptions of the decision curves and surfaces, resp., are unknown.

Conclusions and Outlook
In this article, we presented a method for approximating manifolds of discontinuity of a function f in 2D and 3D and demonstrated successful applications of our method to Multi-criteria Decision Aid, to generic test cases connected with the detection of faults and to an inverse acoustic scattering problem.In all cases, our method requires significantly fewer evaluations of f than previously existing algorithms we compared our method to.At least for the inverse acoustic scattering problem, this leads to a significant acceleration of the overall reconstruction, as computing the factorization method dominates the computational time even in our computations, where the number of such computations could be reduced by a factor of approx.7 to 12. Our algorithm could be easily enhanced with classification algorithms as in [1] in order to consider more general functions than we did, provided the classification is certain.This enables for tackling fault detection problems and for applying our method for interpolation of piecewise smooth functions assuming that f is smooth on Ω 1 , . . ., Ω n with Ω = Ω 1 ∪ . . .∪ Ω n , but globally discontinuous.In [23], the authors propose interpolation with radial basis functions on each Ω i in this setting.This however requires knowledge about the boundaries of each Ω i which could be obtained using our method.For that purpose, [23] provides a fault detection algorithm based on local approximation properties.Combining these two approaches is subject of our current research.

Figure 3 :
Figure 3: Cluster removal by removeclusters in Example 2.4 for S 1,2 (greyed out); points in S (1) 1,2 are displayed in red.This figure is an excerpt of Fig. 4, right.

Figure 4 :
Figure 4: Left: Partition of Ω = [0, 1] 2 according to Example 2.4 indicated by grey solid lines with the initial point set X displayed as grey crosses.We moreover show M (black points) and M 2 (blue points); Right: S (1) 1,2 (red points), S (blue points), and S (cyan-coloured points).

Figure
Figure6: Sorting due to Allasia may fail (left; sorting indicated in blue), whereas our method succeeds in the present situation (right).We obtain two ordered subsets (direction of sorting is indicated by arrows which after combination yields the correct ordering.The dashed line represents a connection to the nearest neighbour rejected due to angle in our approach.

Figure 8 :
Figure 8: Algorithm 2.10 (fill) for Test problem 2.4.Previously existing points are greyed out in S (i) i,j .Otherwise, we stick to the colouring scheme of Fig. 4.

) 2 =
v and search for the roots v min and v max of v 2 + 4v − 16cd.According to Vieta, v min = − 2 + √ 4 + 16cd and v max = −16cd/v min .Resubstituting v yields l max = 4 d/(−cv min ).Some safeguarding of this result leads to a step length of l extra = min{ε gap , β growth d avg , l max } , where we estimate c using Building block 2.8 (estcurv) applied to S loc = {x i,j) 1 , . . ., x (i,j) kextra } with a user-defined parameter k extra .The term β growth d avg increases robustness, as extrapolation is reliable only near the points to extrapolate, and it may happen that d avg ε gap .

Figure 10 :
Figure 10: Sets S(i) i,j (left) and Ŝ(i) i,j (right) for Test problem 2.4.We stick to the colouring scheme of Fig. 4.

Figure 12 :I 1 , 2 I 1 , 2
Figure 12: Left: Gradually filling large gaps in S i,j (displayed as grey dots).The triplet x is represented in golden colour; the triplets in S local , displayed in red, define a local coordinate system sketched with black arrows.The centre of gravity of S local is shown as a blue dot.Right: Expansion of S 1,2 to the inner boundary Γ 1,2 ∩ Γ 1,3 .Triplets in S 1,2 \ S I 1,2 are displayed as grey dots, triplets in S I 1,2 as blue dots.2D-Expanding {x, x } on E ∩ Γ 1,2 yields a new triplet near Γ 1,2 ∩ Γ 1,3 , displayed as red dot.

Example 2 . 27 .
We compute Test problem 2.24 starting with 200 Halton-distributed points and k near = 10.All other parameters are set as in the computation of Test problem 2.4 in

2 , 3
in the visualisation for the sake of clarity.

Figure 15 :
Figure 15: Suitable (left) and unsuitable (right) triplets for expanding to a domain boundary.

Figure 18 :
Figure 18: S (i) i,j (left), S(i) i,j (middle), and Ŝ(i) i,j (right) for Test problem 4.1.As for the computations of Example 2.4 in Section 2.1, previously existing points are greyed out.

Figure 21 :
Figure21: Visualisation of the decision space for 3 (left) and 4 (right) variable criteria.The car users' weightings according to[7] are displayed as red dot along with the circumsphere with radius ρ = 0.06.

Table 1 :
Number of triplets per Γ i,j and number of function evaluations for Example 2.27.
1,3 for Test problem 2.24.Section 2.1.For the corresponding numbers of triplets and the number of function evaluations, we refer to Table 1.Figs. 13, 14 and 16 provide visualisations of the sets of triplets.

Table 2 :
Number of function evaluations per Building block for reconstructing the subdomains.

Table 3 :
Total number of triplets after the respective Building block.Allasia et al. need to employ an explicit classification algorithm, as they consider more general functions f than we do, which may infer additional function evaluations.The number of function evaluations required heavily depends on the accuracy of the approximation and resolution of the fault line, such that comparing the number of function evaluations given a prescribed tolerance for both algorithms would provide a more sound basis for comparing the efficiency of the two algorithms.Detailed information on accuracy is however lacking in Test problem 4.1.

Table 4 :
. In contrast to SAW (Simple Additive Performance matrix P with corresponding criteria generated from

Table 6 :
24: Reconstructed Scatterers long cylinder, peanut, short cylinder, and sphere.Number of function evaluations per Building block for reconstructing the scatterers.

Table 7 :
Total number of triplets after the respective Building block for reconstructing the scatterers.