Thomas Dietenbeck

Ahlem Othmani

Marco Attene

Jean-Marie Favreau

A Framework for Mesh Segmentation and Annotation using Ontologies

Mesh segmentation and annotation using semantics has received an increased interest with the recent democratisation of 3D reconstruction methods. The common approach is to perform this task in two steps, by first segmenting the mesh and then annotating it. However, this approach does not allow one part to take advantage of the other. In image processing, some methods are combining segmentation and annotation, but they are not generic and require implementation adjustments or rewritings for each modification of the expert knowledge. In this work, we describe an original framework that mixes segmentation and annotation while minimizing the required geometric analysis and we give preliminary results showing its feasability.

Our framework provides a generic ontology describing object feature concepts (geometry, topology, etc.) and algorithms allowing to detect these concepts. This ontology can be enlarged by any expert to formally describe a specific domain. The formalized domain description is then used to automatically perform the joint segmentation and annotation of objects and their features, by selecting at each step the most relevant algorithm given the previously detected semantics. This methodology has several advantages. Firsly it allows to segment and annotate objects without any knowledge in mesh or image processing by simply describing the object features in terms of ontological concepts. Secondly this framework can be easily reused and applied to different contexts by simply building on our generic ontology. Finally performing the joint segmentation and annotation allows to use in an efficient way the expert knowledge, reducing possible segmentation errors and the computation time by always launching the most efficient algorithm.

Introduction

During the last two decades, an important work has been done in the data mining and mesh processing communities to integrate a semantic dimension to their work. One of the main goal is to be able to extract an abstract description of the manipulated data, using some semantic descriptors. Bridging the gap between raw data and semantic concepts is a very complicated task, usually implying a good knowledge about the specific applicative domain the systems are working on. This link between the expert knowledge (i.e. the semantics) and the raw data can be achieved using learning techniques or by designing a deterministic system, expressing the knowledge of the expert in a language of computer science.

Techniques for semantic extraction have been explored more precisely for mesh segmentation, where the objects and their subparts can be described very precisely using semantic terms, to describe the shape, the structure or the functionality. A classical problem is to be able to identify an object given its shape, and to recognize each of its subparts by first segmenting them, then by labelling these subparts using the concepts available on the semantic domain of this object class. Some of these approaches [START_REF] Hudelot | Fuzzy spatial relation ontology for image interpretation[END_REF][START_REF] Hassan | Ontology-guided mesh segmentation[END_REF][START_REF] Fouquier | Sequential model-based segmentation and recognition of image structures driven by visual features and spatial relations[END_REF] are able to question the partial semantic description of the scene and to adjust their behaviour to the context. However, we noticed that all these approaches contain in their implementation, algorithms or procedures very specific to the applicative domain.

In this article, we describe an original framework for mesh segmentation that push up the semantic approach, creating a bridge bewteen an expert knowledge description and the segmentation algorithms. This framework allows an expert in a specific domain to formally describe his own domain in terms of a fundamental ontology, without any skill in geometric algorithms. This formalized domain description is then used by the system to automatically recognize objects and their features within that domain.

The genericity of our approach is insured by a multi-layer ontology modeling the expert knowledge. The first layer corresponds to the basic properties of any object, such as shapes and structures. The next layers are specific to each application, describing the functionalities and possible configurations of the objects in this domain. The segmentation and identification mecanism is hidden behind the concepts of the first layer, which are associated to the segmentation algorithms. We wish to highlight here that this work is exploratory and only studies the feasability of the proposed approach on a simple context.

After an overview of the existing methods in semantic-driven expert systems and for mesh segmentation using semantics, we introduce our framework for semantic-driven mesh segmentation and annotation, starting with the expert knowledge modelization, followed by a synthetic catalog of the segmentation algorithms, and finally with the description of the proposed expert system. The fourth section present some experimental results to illustrate the relevancy and feasability of our approach. Finally, we present possible extensions of this work and future improvements.

Related work

The design of enhanced vision systems has benefit from semantic knowledge such as ontology-driven strategy. The ontology paradigm in information science constitutes one of the most diffused tools to make people from various backgrounds work together [START_REF] Seifert | Combined semantic and similarity search in medical image databases[END_REF]. In addition, ontology-based interfaces are a key component of an ergonomic [START_REF] Seifert | Combined semantic and similarity search in medical image databases[END_REF], adaptive computer system [START_REF] Seifert | Combined semantic and similarity search in medical image databases[END_REF], especially in biological/clinical fields [START_REF] Othmani | Ontology-driven image analysis for histopathological images[END_REF] for which concepts and standards are constantly shifting. Furthermore, reasoning capabilities embedded in the logical framework on which ontology softwares are built up should be a definite bridge between computer vision scientists and knowledge engineers. Even though still brittle and limited, reasoning inferences out of visual data may enhance the vision system experience [START_REF] Othmani | Ontology-driven image analysis for histopathological images[END_REF] as well. Computer graphics is one of the different domains which benefits from the semantic knowledge through ontology strategy and for that, different applications are designed to ensure the articulation between ontologies and mesh processing and to ensure a semantic representation in different domains including anatomy [START_REF] Hassan | Ontology-guided mesh segmentation[END_REF], product design in e-manufacturing [START_REF] Attene | Characterization of 3D shape parts for semantic annotation[END_REF], robotic [START_REF] Albrecht | Matching CAD object models in semantic mapping[END_REF][START_REF] Gurau | Challenges in using semantic knowledge for 3D object classification[END_REF].

In [START_REF] Camossi | Deriving functionality from 3D shapes: Ontology driven annotation and retrieval[END_REF], a system to support a user in the retrieval and the semantic annotation of 3D models of objects in different application contexts is presented. The ontology provides a representation of the knowledge needed to infer object's shape, functionality and behavior. Then, the annotation and the retrieval is performed based on the functional and behavior characteristics of the 3D model.

In [START_REF] Attene | Characterization of 3D shape parts for semantic annotation[END_REF], the ontology was used to characterize and to annotate segmented parts of a mesh using a system called "ShapeAnnotator". To this aim, the ontology is loaded according the type of the input mesh and the user can link segments to relevant concepts expressed by the ontology. While the annotation of the mesh parts is done by a simple link with the "ShapeAnnotator", [START_REF] Gurau | Challenges in using semantic knowledge for 3D object classification[END_REF] and [START_REF] Shi | An approach to semi-automatic semantic annotation on Web3D scenes based on an ontology framework[END_REF] proposed to feed an ontology with a set of user-defined rules (e.g. geometric properties of objects, spatial relationships) and the final annotation is created according to them.

In [START_REF] Hassan | Ontology-guided mesh segmentation[END_REF], an ontology including an approximation of the geometric shape of some anatomical organs was used to guide the mesh segmentation. The parameters needed to segment the input mesh were provided by the ontology. For the case of semantic classification, an ontology was used also in [START_REF] Albrecht | Matching CAD object models in semantic mapping[END_REF] into a SLAM-generated 3D point cloud map. After the reconstruction of the surface planes in the point cloud, the ontology is used to generate hypotheses of possible object locations and initial poses estimation and the final result is a hybrid semantic map, in which all identified objects have been replaced by their corresponding CAD models. Recently, [START_REF] Feng | A unified framework for mesh segmentation and part annotation[END_REF] proposed a unified framework which bridges semantics and mesh processing. The mesh is divided into a fixed number of parts corresponding to the number of concepts in the ontology. The parts are then annotated based on some rules defined in the ontology (e.g. the head is very dissimilar with limbs). sponding ontology, without code modification, assuming that the basic concepts he will need are part of the original core of the framework.

The only bridge between algorithms and semantics remains here on the classification and the description of each algorithm. In a first part, we describe our multi-layer ontology paradigm and how the expert knowledge on a specific domain is implemented on top of elementary concepts, then we give the specifications of the elementary segmentation algorithms. Finally, we will describe the way our framework is using these two modules to address the question of semantic-driven mesh segmentation.

Expert knowledge description

We propose in our framework to model the expert knowledge in 2 steps: 1) the semantic concepts associated to the applicative domain are grouped into a multi-layer ontology; 2) the possible combinations of concepts of a given layer to form concepts of a higher level. These combinations correspond to equivalent concepts in the ontology and will be classified at runtime by the ontology reasoner. Since the core of our framework needs to connect segmentation algorithms with the applicative domain ontology, we designed a core ontology, called S 0 that contains all the elementary concepts required for a segmentation process of objects.

Elementary semantic concepts

A segmentation and semantic labelling process implies that the algorithmic part is able to identify and label regions with specific properties, such as geometrical properties (e.g. stick, board, cube, vertical region), color or texture properties (e.g. color uniformity, reflectance, texture patterns), but also properties linked to the position and configuration of subparts with regards to others (e.g. parallel regions, A is up wrt B, A is between B and C). We will refer to these properties as unary (linking one region of the mesh to one concept, e.g. geometry, color, texture) or n-ary (linking several regions together through one concept, e.g. topology, distance).

In the expert knowledge description contained in our framework, these properties are constituting the first layer S 0 of our ontology. In this modelization, we choose to group separately the geometrical and chromatical unary properties and the topological n-ary properties. More specifically, an object property and a range concept is associated to each unary property in the ontology (the domain being the object parts). The range concept (e.g. shape, orientation, color) is then specialized into elementary semantic concepts (e.g. shape → cube, cylinder, sphere, . . . ) which will correspond to the actual value of the object property. Ontologies cannot directly model n-ary properties and we thus represent them using two object properties and a domain concept. The domain concept (e.g. position, distance) is also specialized into elementary semantic concepts (e.g. distance → connected, close, far, . . . ) which are then connected to object parts through 2 object properties (isReferenceRegion and isTargetRegion).

Fig. 1a shows an unary and a n-ary property while Fig. 1b represents a subpart of a possible core ontology with domain / range concepts and their corresponding elementary semantic concepts.

Multi-layer ontology

This first layer (S 0 , blue blocks in Fig. 2) is part of the core of our framework. It is enriched for each applicative context with specific semantic concepts. In section 3.3 we will describe the effective expert system that make the bridge between semantics and algorithms. It requires that the specific semantic concepts are grouped into two (or more) layers: one layer called S 1 (yellow blocks in Fig. 2), using only references to concepts from the S 0 layer, and that describe all the object configurations that can be combinatorially drawn up, i.e the result of a cartesian product between elementary semantic concepts. The supplementary layers (S 2 , . . . , S n , green, red and purple blocks in Fig. 2) describe the combination rules to populate a complete scene.

Linking two layers: Equivalent concepts

In addition to the semantic concepts of the application domain, another important expert knowledge consists in how the concepts are linked with each other. This can be expressed as a set of equivalent concepts of the S n layer describing possible or impossible combinations of concepts of the S n-1 layer. The work of the expert is thus strongly simplified since he can describe not only positive rules, but also negatives. For example in the furniture ontology, a chair leg can be described as a stick shape and a vertical orientation; on the other hand, the combination of a headrest without backrest is an incompatible configuration.

In practice, incompatible configurations are specialized in specific concepts, one for each incompatibility type. This specialization allows to perform some reasoning and classification on partially annotated individuals.

Once these equivalent concepts are given by the expert, they are used in two ways: either to build a decision tree (which will be detailled in Section 3.3) or to suggest segmentation correction. Indeed, during the segmentation / annotation process, incompatible configuration might appear due to either a segmentation error or a missing equivalent concept. In this case, the reasoner can be asked for the reason of the inconsistency which is then presented to the user for correction. The main advantage of this approach is that it allows us to ask the user to correct errors only in the regions that caused the classification as incompatible instead of having to explore the whole mesh / labelling.

Semantic type signatures of algorithms

In this work, the core idea about the algorithmic part is to split the method into elementary algorithms that are dedicated to one of the elementary semantic concepts. But this consideration cannot be the only guide to produce a synthetic catalog of the algorithms: the type signature of these algorithms cannot be only defined by a single concept.

Each of these algorithms will be involved in the segmentation and labellization process by splitting the given region and by adding semantic descriptors to the subparts. Since a region can be described by more than one concept (for example, a part can be a stick and vertical), and since our goal is to have the more elementary algorithms as possible, we deduce that it also exist algorithms that are not splitting the given region, but only increasing the knowledge on this region. Finally, and because we want to deal with concepts that are not necessarly involving a single region, we need to distinguish functions associated to unary properties, and functions associated to n-ary ones.

The following synthetic catalog is a proposal to identify each segmentation algorithm in the context of semantic labelling: semantical questions starting with find all (SF), with unary concepts (e.g. find all rectangles in region A), semantical questions starting with is it a (SI), with unary concepts (e.g. is B a flat region), topological questions (T), identifying n-ary relations between regions (e.g. are B and C connected), topological questions starting with find all (TF), using an n-ary relation and 1 to n -1 regions, to find regions satisfying the relation wrt the input regions (e.g. find all regions connected to the region B; find all regions between the regions B and C).

Each algorithm gets as input one or more regions of the original object and returns a set of regions, each of them enriched by a semantic description generated by the function, plus a score in [0; 1] to illustrate the matching between this region and the associated concept. This first synthetic catalog of possible semantic type signatures covers all the useful algorithms for mesh segmentation using semantic description, as illustrated in the following sections.

Expert system

In this section, we describe how the expert knowledge is used in our algorithm to efficiently segment and annotate an object. Fig. 3a gives an overview on this framework.

One of the advantage of our approach is that it gives the possibility to easily compute a tree containing the order of the questions to ask to reach the solution in the most efficient way. To build this decision tree, we first use the equivalent concepts of each layer to build the set of possible configurations. For each layer and starting at the S 0 layer, the Cartesian product between all the properties of the layer is performed. The reasoner is then used to classify the instances in equivalent concepts and incompatible candidates are removed. The remaining ones are then used in the above layer as semantic concepts in the Cartesian product to compute the new list of possible configurations.

Once the set of all possible configurations Ω is created, the idea is to split it according to the concept maximizing a criterion C. The choice of this concept gives us the question to ask and thus a node of the tree. For each possible answer, we then get the corresponding subset and look for the next concept maximizing C. This operation is iterated until only one possibility is left in each subset. In the resulting tree, the root thus corresponds to Ω and stores the first question to ask, each leaf is a possible configuration and intermediary nodes are subsets of Ω and store the next question to ask. Note that this step can be performed only once and offline in order to speed up the process. This procedure is illustrated in Fig. 3b.

The first inline step of our expert system is to browse the decision tree: starting from the root, the system gets the question to ask and run the associated elementary algorithm on the input object. The system then selects the child node corresponding to the algorithm's result and the process iterates until a leaf of the tree is reached meaning that the semantics associated to the object is known as only one possible configuration remains.

In some cases, the expert system might reach a leaf before every part of the mesh is segmented or annotated (because some concepts might be inferred from the presence / absence of others). Some algorithms are then run on the mesh to confirm the global semantics and annotate the missing parts. This can again be done very efficiently by questionning the reasoner, which will give the expert system the missing concepts and thus the elementary algorithm to run.

Experiments on Furnitures Segmentation and Annotation

The experiments have been done on furnitures segmentation and annotation, using basic shapes relevant to this domain (see Fig. 4b). In this work, we focus our interest on the expert system in order to study the feasability of the presented framework. Finely selecting and adjusting elementary segmentation algorithms will be one of the future work we mention in section 5. We implemented the expert system detailed in section 3.3 using Java and the OWL API, and designed our prototype such that the purely mesh manipulations are written in C++, using CGAL. The connection between these two parts is done using a client/server paradigm via sockets.

To address the specific context of furnitures, we designed a dedicated ontology (see Fig. 4a) using Protégé1 on top of a simplified version of the elementary concepts introduced in section 3.1.1. The decision tree associated to this ontology is generated following the procedure described in section 3.3 where the criterion C was chosen to be the dichotomy, i.e. we look for concepts allowing to split the set into two subset of same size. Once the decision tree is generated, we can run our expert segmentation system on meshes. Fig. 5 gives the list of questions that are computed in order to segment and recognize the first bench in Fig. 4b. The other images in Fig. 4b are illustrating the segmentation and identification process using the same expert ontology with various meshes.

We have compared our segmentation and identification method with other equivalent approaches, and we summarize the results in Tab. 1. The first column, titled S. then I. corresponds to an approach where a first segmentation preprocessing is done to split regions, then each region is labeled using the semantic concepts. The second column titled Naive S&I corresponds to an approach where the segmentation algorithms dedicated to specific kind of shapes are run independently to identify the regions. The last column corresponds to our approach. For each method, we detailed the number of semantic find all (SF) and semantic is it a (SI) algorithms required to segment and annotate the 5 objects shown in Fig. 4b.

We choose to distinguish the 2 kinds of algorithms, because the complexity of each algorithm family is significantly different: a SF algorithm will require to browse all the given regions (possibly the whole mesh), and will have to extract subparts from it. In comparison, a SI algorithm will only have to validate or not a feature on a given region. Minimizing the number of SF runs is thus the main goal of a segmentation and identification process.

For each region of an object, the expert ontology uses 3 range concepts (shape, position, orientation) that implies 8 elementary semantic concepts. Each object is also characterized by a supplementary range concept (number of seats) that implies 2 elementary semantic concepts per object. The number of SI algorithms in the S. then I. has been estimated counting for each subregion of the mesh one SI per elementary semantic concept (minus 1 per range concept that can be deduced), plus the 1 elementary semantic concepts of the full object. The Naive S&I consist in running all the available SF algorithms. The number of SF and SI of our method comes from the trace of the experimental runs (see an example in Fig. 5).

We first compare our work to a S. then I. approach, where a segmentation preprocessing is applied before the identification. We cannot quantitatively compare this approach with ours, but since we use the expert knowledge to reduce the number of SF algorithms in our approach, we can deduce that our SF computations are almost equivalently expensive as the preprocessing stage of the S. then I. approach2 . The number of SI algorithms is also strongly reduced with our approach.

The second approach for comparison is a Naive S&I approach, where all the SF algorithms are run. The number of SF algorithms is strongly reduced by our approach, and the comparison between the Naive S&I approach can be summarized with a comparison between the complexity of SF and the complexity of SI algorithms.

These first results are illustrating the relevancy of our method with respect to the existing approaches: mixing segmentation and identification steps is a good approach to reduce the complexity of the global algorithm.

Conclusion and future work

In this paper, we presented a new framework for efficient segmentation and annotation of meshes. It is composed of two blocks: a multi-layer ontology gathering the semantics on the application domain and a processing part allowing to detect elementary geometrical, chromatical and topological concepts. The main advantage of our method is that it separates the knowledge on the domain with the processing allowing an expert to segment and annotate an object without knowledge in image or mesh processing. Another advantage is that using the expert knowledge, we are able to build a decision tree to perform an efficient search amongst the set of possible objects while being able to suggest segmentation and annotation corrections to the user if an impossible configuration is reached.

The ontology we designed for this experiment is very basic, and we plan to experiment during the next months a more complete ontology with more concepts. In particular, the n-ary properties will be integrated in order to express more realistic constraints between subparts of the objects.

The simulation we presented in section 4 is using manual segmentations to mimic the result of the algorithms. Our next step will be to select and adjust elementary algorithms from the litterature, like geometrical extractors [START_REF] Mortara | Plumber: a method for a multi-scale decomposition of 3D shapes into tubular primitives and bodies[END_REF][START_REF] Li | Globfit: Consistently fitting primitives by discovering global relations[END_REF] or segmentation approaches driven by functionality features [START_REF] Laga | Geometry and context for semantic correspondence and functionality recognition in manmade 3D shapes[END_REF].

Introducing automatic segmentations will open many related questions we plan to handle in a near future. One of the next challenges will be to introduce better approaches to choose between algorithms than the dichotomical one. A first criterion to consider could be a weighting system that favour algorithms with a small computational time, or to include the exactness of the algorithms. These weights will be introduced in the decision tree computation in order to design an expert system that handle the question of efficiency.

Applying our approach on meshes acquired from low resolution devices will complicate the job of the segmentation algorithms. It will probably generate incoherent subregions, with overlappings or unlabelled parts. One possible approach is to use fuzzy maps to describe the regions, but another question will have to be handled: how to adjust an existing partial segmentation? Our framework is a good candidate to provide a specific answer to this problem, since the expert knowledge contains information about the expected configurations. One possible extension of this work could be to introduce adjustment algorithms for each elementary concepts, that will be able to adjust a first segmentation using the global knowledge of a specific domain.

Finally, a long term extension of this work will be to introduce it into a machine learning approach, where the ontology will be deduced or extended from an existing one, using a set of shapes of the domain. This extended framework will be a possible challenger of the 3D Shape Retrieval Contest (SHREC) organized each year in the mesh segmentation community.

  FIG. 1: (a) Examples of unary and n-ary property. A and B correspond to region of the object. (b) A subpart of the elementary semantic concepts of the core ontology. The number beside some concepts indicate that they link several regions together (e.g. UpFrom relates 2 regions A and B).

  FIG.2: Examples of multi-layer ontologies for furnitures (a) and street (b) segmentation. Note how a more detailled expert knowledge descriptions of the same domain can be achieved by simply adding a layer (e.g. for the concepts of Back or Support in the case of Furnitures)

  FIG.3: Overview of the proposed method.

  FIG.4: Ontology and meshes used for our experiments on furniture segmentation and annotation.

  DownPosition R1: false 4. SI DownPosition R2: true 5. SI DownPosition R3: true 6. SI StickShape R2: false 7. SI StickShape R3: false 8. SI DownPosition R1: false 9. SI BoardShape R1: true 10. SI BoardShape R2: true 11. SI BoardShape R3: true 12. SI BoardShape R0 -{R1, R2, R3}

•

  R4: true FIG.5: List of the elementary algorithms and corresponding answers generated by our expert system to segment and annotate a bench mesh. SF: semantic find all, SI: semantic is it a.

			S. then I.	Naive S&I	Our method
	Object	#subparts	preprocessing	-		-
			# SF	# SI	# SF	# SI	# SF	# SI
	bench 1	4	0	21	9	0	2	10
	bench 2	6	0	31	9	0	3	13
	couch	6	0	31	9	0	3	14
	chair 1	6	0	31		0	2	15
	chair 2	8	0	41	9	0	2	19
	TAB. 1: Number of segmentation steps for a complete identification and segmentation.

Proposed methodSegmenting an object following its geometry is a non trivial problem, as well as associating semantic concepts to each object and its subparts. Both questions need a very specific processing. Previous works[START_REF] Hudelot | Fuzzy spatial relation ontology for image interpretation[END_REF][START_REF] Attene | Characterization of 3D shape parts for semantic annotation[END_REF][START_REF] Hassan | Ontology-guided mesh segmentation[END_REF][START_REF] Fouquier | Sequential model-based segmentation and recognition of image structures driven by visual features and spatial relations[END_REF] on this topic are more and more going in the direction of mixing the two problems, in order to help the segmentation using the already extracted semantics, and by extracting the semantics from the partial segmentations.The framework we present in this section addresses the same goal, with a strong focus on the separation between segmentation algorithms and the semantic reasonings. The benefits of this approach will be experimented in the next section, and the possible extensions discussed in the last one, but we can already underline one direct benefit of this approach: using this framework to address a new applicative domain will only require the user to design the corre-

http://protege.stanford.edu/

The best strategy for a preprocessing can be to choose the smallest range concept, then run for each elementary concept a SF algorithm.

Acknowledgements Marco Attene thanks the EU FP7 Project VISIONAIR for having supported his contributions in this research.


R G Novikov
email: novikov@cmap.polytechnique.fr

Formulas for phase recovering from phaseless scattering data at fixed frequency

Formulas for phase recovering from phaseless scattering data at fixed frequency

Introduction

We consider the equation

-∆ψ + v(x)ψ = Eψ, x ∈ R d , d ≥ 2, E > 0, (1.1)
where ∆ is the Laplacian, v is a scalar potential such that

v ∈ L ∞ (R d ), supp v ⊂ D,
D is an open bounded domain in R d .

(1.2)

Equation (1.1) can be considered as the quantum mechanical Schrödinger equation at fixed energy E.

Equation (1.1) can also be considered as the acoustic equation at fixed frequency ω. In this setting

E = ( ω c 0 ) 2 , v(x) = (1 -n 2 (x)) ( ω c 0 ) 2 , (1.3)
where c 0 is a reference sound speed, n(x) is a scalar index of refraction.

For equation (1.1) we consider the classical scattering solutions ψ + continuous and bounded on R d and specified by the following asymptotics as |x| → ∞: (d, |k|) e i|k||x| |x| (d-1)/2 f (k, |k|

ψ + (x, k) = e ikx + c
x |x| ) + O ( 1 |x| (d+1)/2 ) , x ∈ R d , k ∈ R d , k 2 = E, c(d, |k|) = -πi(-2πi) (d-1)/2 |k| (d-3)/2 , (1.4)
where a priori unknown function f = f (k, l), k, l ∈ R d , k 2 = l 2 = E, arising in (1.4) is the classical scattering amplitude for (1.1).

In order to find ψ + and f from v one can use the following Lippmann-Schwinger integral equation (1.5) and formula (1.7) (see, e.g., [BS], [FM]):

ψ + (x, k) = e ikx + ∫ D G + (x -y, k)v(y)ψ + (y, k)dy,
(1.5) where

G + (x, k) def = -(2π) -d ∫ R d e iξx dξ ξ 2 -k 2 -i0 = G + 0 (|x|, |k|), (1.6) where x ∈ R d , k ∈ R d , k 2 = E,
k ∈ R d , l ∈ R d , k 2 = l 2 = E.
We recall that ψ + describes scattering of the incident plane waves e ikx on the potential v. And the second term of the right-hand side of (1.4) describes the scattered spherical waves.

In addition to ψ + , we consider also the function R + describing scattering of spherical waves generated by point sources. The function

R + = R + (x, x ′ , E), x ∈ R d , x ′ ∈ R d , can be defined as the Schwartz kernel of the standard resolvent (-∆ + v -E -i0) -1 . Note that R + = R + (x, x ′ , E) = -G + 0 (|x -x ′ |, √ E)
for v ≡ 0, where G + 0 is the function of (1.6). Given v, to determine R + one can use, in particular, the following integral equation

R + (x, x ′ , E) = -G + 0 (|x -x ′ |, √ E) + ∫ D G + 0 (|x -y|, √ E)v(y)R + (y, x ′ , E)dy, (1.8) where x ∈ R d , x ′ ∈ R d . The function R + (x, x ′ , E) at fixed x ′ ∈ R d describes scattering of the spherical wave -G + 0 (|x -x ′ |, √ E)
generated by a point source at x ′ . In addition,

R + (x, x ′ , E) = - c(d, √ E) (2π) d e i √ E|x| |x| (d-1)/2 ψ + ( x ′ , - √ E x |x| ) + O ( 1 |x| (d+1)/2 ) as |x| → ∞ at fixed x ′ ,
(1.9)

R + (x, x ′ , E) = R + (x ′ , x, E), (1.10)
where c is the constant of (1.4), ψ + is the function of (1.4), (1.5).

In connection with aforementioned facts concerning R + see, e.g., Section 1 of Chapter IV of [FM].

Let S d-1 r = {k ∈ R d : |k| = r}, r > 0.
(1.11).

We consider the following three types of scattering data for equation (1.1):

(a) f (k, l), where (k, l) ∈ Ω ′ f ⊆ Ω f , Ω f = S d-1 √ E × S d-1 √ E ; (1.12a) (b) ψ + (x, k), where (x, k) ∈ Ω ′ ψ ⊆ Ω ψ , Ω ψ = (R d \(D ∪ ∂D)) × S d-1 √ E , (1.12b) assuming (1.13); (c) R + (x, y, E), where (x, y) ∈ Ω ′ R ⊆ Ω R , Ω R = (R d \(D ∪ ∂D)) × (R d \(D ∪ ∂D)), (1.12c)
assuming (1.13).

In (1.12b), (1.12c) we assume also that

R d \(D ∪ ∂D) is connected. (1.13)
We consider the following inverse scattering problems for equation (1.1) at fixed As regards to acoustic scattering experiments in the framework of the model described by (1.1), (1.3), the complete scattering data f , ψ + , R + of Problems 1.2 can be measured directly. Nevertheless, in some cases it may be more easy to measure the phaseless versions of these data. Therefore, Problems 1.2 are also of interest in the framework of acoustic inverse scattering.

On the other hand, in the literature many more results are given on Problems 1.1 (see [ABR], [Be], [Bu], [BAR], [BSSR], [ChS], [E], [F1], [G], [HH], [HN], [I], [IN], [M], [Na], [N1]- [N6], [R], [S] and references therein) than on Problems 1.2 (see Chapter X of [ChS] and recent works [KR], [N7] and references therein).

The works [K1], [K2], [K3], [KR], [N7] give also results on analogs of Problems 1.2, where E is not fixed. Besides, analogs of Problems 1.2 in dimension d = 1, where E is not fixed, were considered, in particular, in [AS], [KS].

Let

B r = {x ∈ R d : |x| < r}, r > 0. (1.14)
Suppose that, for some r > 0, D ⊆ B r .

(1.15)

In connection with Problems 1.1, under assumption (1.15), it is well known, in particular, that any of the scattering data (a) f on Ω f , (b) ψ + on ∂B r × S d-1 √ E , or (c) R + on ∂B r × ∂B r uniquely and constructively determine two other data; see [Be].

In addition, in [N1] for d ≥ 3 (see also [N3]) and in [Bu] for d = 2 it was shown that f on Ω f at fixed E uniquely and constructively determines v on R d , under assumption (1.2). For related exact stability estimates, see [S], [HH], [IN], [I]. Besides, for approximate but efficient methods for solving Problem 1.1a, see [N4], [N5], [ABR], [BAR], [N6].

The main results of the present work can be summarized as follows.

First, we give explicit asymptotic formulas for finding

f (k, l) at fixed (k, l) ∈ Ω f , k ̸ = l, from |ψ + (x, k)| 2 for x = sl/|l|, s ∈ Λ = [r 1
, +∞[ for arbitrary large r 1 ≥ r (assuming, e.g., (1.15)); see Theorem 2.1 and Corollary 2.1 of Section 2.

In addition, we have the asymptotic formula (2.11) for finding |ψ +∞[ (assuming, e.g., (1.15)).

+ (x ′ , k)| 2 at fixed (x ′ , k) ∈ Ω ψ from |R + (x, x ′ , E)| 2 for x = -sk/|k|, s ∈ Λ for any unbounded Λ ⊂ [r,
The aforementioned formulas give explicit reductions of Problems 1.2b, 1.2c to Problem 1.1a for appropriate Ω ′ ψ , Ω ′ R and Ω ′ f . In connection with reductions of Problems 1.2b, 1.2c to Problem 1.1a we give also additional global uniqueness results summarized in Theorem 2.2 of Section 2.

Second, we give global uniqueness results for Problem 1.2b for the case when Ω ′ ψ is an open subset of Ω ψ and for Problem 1.2c for the case when Ω ′

R is an open subset of Ω R ; see Theorem 2.3 of Section 2. In this connection we recall also that for Problem 1.1a in its initial formulation there is no uniqueness, in general; see [N7].

Actually, as soon as Problems 1.2b, 1.2c are reduced to Problem 1.1a, one can use all known results for exact or approximate solving Problem 1.1a; see [N6] and other works cited above in connection with Problems 1.1.

Finally, we indicate some possible generalizations and extentions of results of the present work; see Remarks 2.1-2.4 at the end of Section 2. In particular, in a subsequent work we plan to consider phaseless inverse scattering in dimension d = 1 when E is not fixed, using an analog of Theorem 2.1 for d = 1.

To our knowledge, no exact general result on phase recovering from phaseless scattering data for equation (1.1) at fixed E was given in the literature before the present work.

The main results of the present work are presented in detail in the next section.

Main results

We represent f and c of (1.4) as follows: iβ(d,|k|) .

f (k, l) = |f (k, l)|e iα(k,l) , c(d, |k|) = |c(d, |k|)|e
(2.1)

We consider

a(x, k) = |x| (d-1)/2 (|ψ + (x, k)| 2 -1), (2.2) a 0 (x, k) = 2Re ( c(d, |k|)e i(|k||x|-kx) f ( k, |k| x |x| )) , (2.3) δa(x, k) = a(x, k) -a 0 (x, k) (2.4) for x ∈ R d \{0}, k ∈ R d \{0}
, where ψ + are the scattering solutions of (1.4), (1.5), f is the scattering amplitude of (1.4), (1.7), c is the constant of (1.4).

For real-valued potential v satisfying (1.2) one can show that

|δa(x, k)| ≤ δ 0 (|x|, |k|), x ∈ R d , k ∈ R d \{0}, (2.5) δ 0 (r, ρ) = O(r -1/2 ) for r → +∞, d = 2, δ 0 (r, ρ) = O(r -1 ) for r → +∞, d ≥ 3, (2.6)
at fixed ρ > 0 (where δ 0 depends also on v). Estimates (2.5), (2.6) are proved in Section 3.

The key result of the present work consists in the following theorem:

Theorem 2.1. Let real-valued potential v satisfy (1.2), d ≥ 2, and f , a be the functions of (1.4), (2.2). Let (k, l) ∈ Ω f of (1.12a), k ̸ = l, and

T = 2π ( E 1/2 ( 1 - kl E )) -1 .
(2.7)

Then the following formulas hold:

|f | ( cos α sin α ) = (2|c| sin(2πT -1 (s 1 -s 2 ))) -1 × ( -sin(2πT -1 s 2 + β) sin(2πT -1 s 1 + β) -cos(2πT -1 s 2 + β) cos(2πT -1 s 1 + β) ) × (( a((s 1 + nT )l/|l|, k) a((s 2 + nT )l/|l|, k) ) - ( δa((s 1 + nT )l/|l|, k) δa((s 2 + nT )l/|l|, k)
))

, n ∈ N,

(2.8)

s 1 , s 2 ∈ [0, T ], s 1 ̸ = s 2 (mod T /2), α = α(k, l), |f | = |f (k, l)|, |c| = |c(d, |k|)|, β = β(d, |k|),
where α, β are the angles of (2.1), c is the constant of (1.4), δa is defined by (2.4),

|δa(s j + nT )l/|l|, k)| ≤ δ 0 (s j + nT, |k|) = { O ( (nT ) -1/2 ) , n → ∞, d = 2, O ( (nT ) -1 ) , n → ∞, d ≥ 3, (2.9)
where j = 1, 2, δ 0 is the function of (2.5).

Theorem 2.1 is proved in Section 3.

In connection with formula (2.8) we consider

Λ = ∪ +∞ j=1 (s 1 + n j T ) ∪ (s 2 + n j T ), n j ∈ N, n j → +∞, j → +∞, (2.10) for fixed T > 0 and s 1 , s 2 ∈ [0, τ ], s 1 ̸ = s 2 (mod T /2).
Corollary 2.1. Let real-valued potential v satisfy (1.2), d ≥ 2, and ψ + , f be the scattering functions of (1.4). Let (k, l) ∈ Ω f of (1.12a), k ̸ = l, and T be defined by (2.7). Then:

|ψ + (sl/|l|, k)| 2 , s ∈ Λ of (2.10), uniquely determines f (k, l) via (2.1), (2.8), (2.9).

In addition, due to (1.9), we have that

|ψ + (x ′ , k)| 2 = (2π) d |c(d, |k|)| -1 s (d-1)/2 |R + (-sk/|k|, x ′ , k 2 )| 2 + O(s -1 ) (2.11) for s → +∞ at fixed x ′ ∈ R d , k ∈ R d \{0}.
Formulas (2.2), (2.8),(2.9), (2.11) give explicit reductions of Problems 1.2b, 1.2c to Problem 1.1a for appropriate Ω ′ ψ , Ω ′ R and Ω ′ f . In this connection we have also the uniqueness results of the following theorem:

Theorem 2.2. Let real-valued potential v satisfy (1.2), (1.13), (1.15), d ≥ 2, and ψ + , f , R + be the scattering functions of (1.4), (1.5), (1.7), (1.8), (1.9). Then:

(1) the scattering amplitude

f (k, l) for fixed (k, l) ∈ Ω f , k ̸ = l, is uniquely determined by the phaseless scattering data |ψ + (sl/|l|, k)| 2 , s ∈ [r 1 , r 2 ], r ≤ r 1 < r 2 ,
(2) the scattering amplitude f (k, l) for fixed (k, l) ∈ Ω f , k ̸ = l, is uniquely determined by the phaseless scattering data

|R + (-sk/|k|, s ′ l/|l|, E)| 2 , (s, s ′ ) ∈ [r 1 , r 2 ] × [r ′ 1 , r ′ 2 ], r ≤ r 1 < r 2 , r ≤ r ′ 1 < r ′ 2 ,
where E and r are the parameters of (1.12a), (1.15);

(3) the phaseless scattering data |ψ + | 2 on a fixed open subset Ω ′ ψ of Ω ψ uniquely determine the scattering amplitude f on Ω f at fixed E, (4) the phaseless scattering data

|R + | 2 on a fixed open subset Ω ′ R of Ω R uniquely determine the scattering amplitude f on Ω f at fixed E.
Theorem 2.2 is proved in Section 3. This proof involves real-analytic continuations.

As a corollary of results of [N1], [N3], [Bu] and items (3) and (4) of Theorem 2.2, we have also the following global uniqueness results on Problems 1.2b, 1.2c:

Theorem 2.3. Let real-valued potential v satisfy (1.2), (1.13), d ≥ 2. Then: • the phaseless scattering data |ψ + | 2 on a fixed open subset Ω ′ ψ of Ω ψ uniquely deter- mine v in L ∞ (R d ),
• the phaseless scattering data

|R + | 2 on a fixed open subset Ω ′ R of Ω R uniquely determine v in L ∞ (R d ).
Remark 2.1. In all aforementioned results of this section the assumption that v is real-valued can be replaced by the assumption that v is complex-valued and equation (1.5) is uniquely solvable for

ψ + ∈ L ∞ (D) at fixed k ∈ R d , k 2 = E > 0.
Remark 2.2. In Problems 1.1, 1.2 the assumption that v of (1.1) is compactly supported (supported in D) can be replaced by the assumption that v has sufficient decay at infinity; see e.g. [ChS], [E], [F1], [G], [HN], [N1]- [N6], [R], [VW], [W], [WY]. In this case, especially in Problems 1.1b, 1.1c, 1.2b, 1.2c, it is natural to assume that v is a priori known on R d \D. Theorem 2.1, Corollary 2.1 and formula (2.11) remain valid for v with sufficient decay at infinity. Remark 2.3. Theorem 2.1, Corollary 2.1 and formula (2.11) have analogs in dimension d = 1. Using these results, in a subsequent work we plan to consider phaseless inverse scattering in dimension d = 1 when frequency is not fixed.

Remark 2.4. The approach of the present work can be also used for phaseless inverse scattering for obstacles.

Proofs of estimates (2.5), (2.6) and Theorem 2.1

3.1. Proofs of estimates (2.5), (2.6). Note that the asymptotic formula (1.4) holds uniformly in x/|x| ∈ S d-1 . Therefore, we have that

ψ + (x, k) = ψ + 1 (x, k) + δψ + (x, k), (3.1) ψ + 1 (x, k) = e ikx + c(d, |k|) e i|k||x| |x| (d-1)/2 f (k, |k| x |x| ), (3.2
)

|δψ + (x, k)| ≤ δ 1 (|x|, |k|), (3.3) δ 1 (r, ρ) = O ( r -(d+1)/2 ) for r → +∞ at fixed ρ > 0, (3.4)
where

x ∈ R d , k ∈ R d \{0}.
Using (2.3), (3.1), (3.2) we obtain that

|ψ + (x, k)| 2 = ψ + (x, k)ψ + (x, k) = 1 + |x| -(d-1)/2 a 0 (x, k) + |x| -(d-1) |c(d, |k|)| 2 |f (k, |k|x/|x|)| 2 + 2Re (δψ + (x, k)ψ + 1 (x, k)) + |δψ + (x, k)| 2 .
(3.5)

Due to (2.2)-(2.4), (3.5), we have that

δa(x, k) = |x| -(d-1)/2 |c(d, |k|)| 2 |f (k, |k|x/|x|)| 2 + 2|x| (d-1)/2 Re (δψ + (x, k)ψ + 1 (x, k)) + |x| (d-1)/2 |δψ + (x, k)| 2 . (3.6) Note that |f (k, l)| ≤ C f ( √ E), (k, l) ∈ Ω f , (3.7) for some positive C f ( √ E),
where Ω f is defined by (1.12a). Using (3.3), (3.6), (3.7) we obtain (2.5) with δ 0 (r, ρ) given by

δ 0 (r, ρ) = δ 0,1 (r, ρ) + δ 0,2 (r, ρ) + δ 0,3 (r, ρ), δ 0,1 (r, ρ) = r -(d-1)/2 |c(d, ρ)| 2 (C f (ρ)) 2 , δ 0,2 (r, ρ) = 2r (d-1)/2 δ 1 (r, ρ)(1 + c(d, ρ)r -(d-1)/2 C f (ρ), δ 0,3 (r, ρ) = r (d-1)/2 (δ 1 (r, ρ)) 2 , (3.8)
where r > 0, ρ > 0.

Formulas (3.4), (3.8) imply (2.6). This completes the proof of estimates (2.5), (2.6).

3.2. Proof of Theorem 2.1. Due to (2.1), (2.3), (2.7) we have that

a 0 (sl/|l|, k) = 2|c(d, |k|)||f (k, l)| cos(2πT -1 s + α(k, l) + β(d, |k|)), (k, l) ∈ Ω f , k ̸ = l, s > 0, (3.9)
where Ω f is defined by (1.12a), T is defined by (2.7). In addition, one can see that

T > 0 for (k, l) ∈ Ω f , k ̸ = l. (3.10)
Due to (2.4), (3.9), we have that

|f | ( cos(2πT -1 s + β) cos α -sin(2πT -1 s + β) sin α ) = (2|c|) -1 (a(sl/|l|, k) -δa(sl/|l|, k)), (k, l) ∈ Ω f , k ̸ = l, s > 0, (3.11) α = α(k, l), |f | = |f (k, l)|, |c| = |c(d, |k|)|, β = β(d, |k|). Using (3.11) for s = s 1 + nT and s = s 2 + nT , where s 1 , s 2 ∈ [0, T ], n ∈ N, we obtain the system ( cos(2πT -1 s 1 + β) -sin(2πT -1 s 1 + β) cos(2πT -1 s 2 + β) -sin(2πT -1 s 2 + β) ) |f | ( cos α sin α ) = (2|c|) -1 ( a((s 1 + nT )l/|l|, k) -δa((s 1 + nT )l/|l|, k) a((s 2 + nT )l/|l|, k) -δa((s 2 + nT )l/|l|, k)
) .

(3.12) Formula (2.8) follows from (3.12). Formula (2.9) follows from (2.5), (2.6). Theorem 2.1 is proved.

Proof of Theorem 2.2

Note that

|ψ + | 2 = ψ + ψ + , (4.1) |R + | 2 = R + R + . (4.2) Note also that R d \(B r ∪ ∂B r ) ⊆ R d \(D ∪ ∂D) (4.3)
under assumption (1.15).

Proof of item (1)

. Due to equation (1.1) for ψ + and assumptions (1.2), we have that

-∆ψ + (x, k) = Eψ(x, k), x ∈ R d \(D ∪ ∂D), for each k ∈ S d-1 √ E . (4.4) Therefore, ψ + (•, k) is (complex -valued) real -analytic on R d \(D ∪ ∂D) at fixed k ∈ S d-1 √ E .
(4.5)

Using (4.1), (4.4) we obtain that

|ψ + (•, k)| 2 is real -analytic on R d \(D ∪ ∂D) at fixed k ∈ S d-1 √ E . (4.6)
As a corollary of (4.3), (4.6), we have that

|ψ + (sl/|l|, k)| 2 is real -analytic in s ∈]r, +∞[ at fixed (k, l) ∈ Ω f . (4.7)
Therefore, at fixed (k, l) ∈ Ω f , the function |ψ + (sl/|l|, k)| 2 given for all s ∈ [r 1 , r 2 ], where r ≤ r 1 < r 2 , uniquely determines this function for all s ∈]r, +∞[ via real-analytic continuation. This result and Corollary 2.1 imply item (1) of Theorem 2.2. 2). Note that

Proof of item (

(-∆ x + v(x) -E)R + (x, x ′ , E) = δ(x -x ′ ), (-∆ x ′ + v(x) -E)R + (x, x ′ , E) = δ(x -x ′ ), (4.8) where x ∈ R d , x ′ ∈ R d , E > 0, δ is the Dirac delta-function.
In addition, due to (1.2), (4.8), we have that

(-∆ x -E)R + (x, x ′ , E) = δ(x -x ′ ), x ∈ R d \(D ∪ ∂D), x ′ ∈ R d , (-∆ x ′ -E)R + (x, x ′ , E) = δ(x -x ′ ), x ∈ R d , x ′ ∈ R d \(D ∪ ∂D).
(4.9) Therefore,

R + (•, x ′ , E) is (complex -valued) real -analytic on R d \(D ∪ ∂D ∪ {x ′ }) for fixed x ′ ∈ R d , R + (x, •, E) is (complex -valued) real -analytic on R d \(D ∪ ∂D ∪ {x ′ }) for fixed x ∈ R d .
(4.10) Using (4.2), (4.10) we obtain that

|R + (•, x ′ , E)| 2 is real -analytic on R d \(D ∪ ∂D ∪ {x ′ }) for fixed x ′ ∈ R d , |R + (x, •, E)| 2 is real -analytic on R d \(D ∪ ∂D ∪ {x ′ }) for fixed x ∈ R d . (4.11)
As a corollary of (4.3), (4.11), we have, in particular, that .12) If l ̸ = -k, then using (2.11), (4.12) one can see that the phaseless data of item (2) uniquely determine |ψ

|R + (-sk/|k|, s ′ l/|l|, E)| 2 is real -analytic in s ∈]r, +∞[ if l ̸ = -k and in s ∈]r, +∞[\s ′ if l = -k at fixed s ′ ∈ [0, +∞[. ( 4
+ (s ′ l/|l|, k)| 2 , s ′ ∈ [r ′ 1 , r ′ 2 ]
, via real-analytic continuation in s. In turn, the latter data uniquely determine f (k, l) due to the result of item (1). Therefore, it remains to consider the case when l = -k. In addition, in view of (1.10), we can assume that r ′ 2 ≤ r 2 . Under these conditions, using (2.11), (4.12) one can see that phaseless data of item (2) uniquely determine |ψ

+ (s ′ l/|l|, k)| 2 , s ′ ∈ [r ′ 1 , r ′ 2 [
, via real-analytic continuation in s. In turn, the latter data uniquely determine f (k, l) due to the result of item (1).

This completes the proof of item (2) of Theorem 2.2.

Proof of item (3).

Note that, under our assumptions,

ψ + (x, •) is (complex -valued) real -analytic on S d-1 √ E at fixed x ∈ R d . (4.13)
Actually, this result is well-known and follows from consideration of (1.5) with k ∈ C d , k 2 = k 2 1 + . . . + k 2 d = E. Using (4.1), (4.13) we obtain that

|ψ + (x, •)| 2 is real -analytic on S d-1 √ E at fixed x ∈ R d . (4.14) Let A k ′ ,ε = {k ∈ S d-1 √ E : |k -k ′ | < ε}, k ′ ∈ S d-1 √ E , ε > 0, (4.15) B x ′ ,ε = {x ∈ R d : |x -x ′ | < ε}, x ′ ∈ R d , ε > 0. (4.16)
Since Ω ′ ψ is open in Ω ψ , we can take

x ′ ∈ R d \(D ∪ ∂D), k ′ ∈ S d-1 √ E , ε 1 > 0, ε 2 > 0 such that B x ′ ,ε 1 × A k ′ ,ε 2 ⊆ Ω ′
ψ . In addition, using (1.13), (4.6), (4.14) one can see that already |ψ + | 2 on B x ′ ,ε 1 ×A k ′ ,ε 2 uniquely determines |ψ + | 2 on Ω ψ via sequential real-analytic continuations in x and in k.

In turn, |ψ + | 2 on Ω ψ uniquely determines f for all (k, l) ∈ Ω f , k ̸ = l, due to Corollary 2.1. In addition, under our assumptions, f is continuous on Ω f and therefore, is uniquely determined also for k = l.

This completes the proof of item (3) of Theorem 2.2.

Proof of item (4).

Since

Ω ′ R is open in Ω R , we can take y ∈ R d \(D ∪ ∂D), y ′ ∈ R d \(D ∪ ∂D), ε > 0, ε ′ > 0 such that B y,ε × B y ′ ,ε ′ ⊆ Ω ′ R .
In addition, using (1.13), (4.11) one can see that already |R + | 2 on B y,ε × B y ′ ,ε ′ uniquely determines |R + (x, x ′ , E)| 2 for all (x, x ′ ) ∈ Ω R , x ̸ = x ′ , via sequential real-analytic continuations in x and in x ′ . In turn, the latter data uniquely determine |ψ + | 2 on Ω ψ via (2.11). Finally, due to item (3), |ψ + | 2 on Ω ψ uniquely determines f on Ω f . This completes the proof of item (4) of Theorem 2.2.

Acknowledgment. The author is very grateful to M.V. Klibanov for remarks that have helped to simplify considerably the formulations of Theorem 2.1 and Corollary 2.1 of the first version (see [N8]) of the present article.

  and G + 0 depends also on d;f (k, l) = (2π) -d ∫ D e -ily v(y)ψ + (y, k)dy,(1.7)

Problem 1.1c. Reconstruct

  Reconstruct potential v on R d from its scattering data ψ + on some appropriate Ω ′ ψ ⊆ Ω ψ . potential v on R d from its scattering data R + on some appropriate Ω ′ R ⊆ Ω R .

	E:
	Problem 1.1a. Reconstruct potential v on R d from its scattering amplitude f on
	some appropriate Ω ′ f ⊆ Ω f .
	Problem 1.1b.

Problem 1.2a. Reconstruct

  potential v on R d from its phaseless scattering data |f | 2 on some appropriate Ω ′ f ⊆ Ω f .

Problem 1.2b. Reconstruct

  potential v on R d from its phaseless scattering data |ψ + | 2 on some appropriate Ω ′ ψ ⊆ Ω ψ .

Problem 1.2c. Reconstruct

  potential v on R d from its phaseless scattering data |R + | 2 on some appropriate Ω ′ R ⊆ Ω R .

	Note that in quantum mechanical scattering experiments in the framework of model
	described by equation (1.1) the phaseless scattering data |f | 2 , |ψ + | 2 , |R + | 2 of Problems
	1.2a-1.2c can be measured directly, whereas the complete scattering data f , ψ + , R + of
	Problems 1.1a-1.1c are not accessible for direct measurements. Therefore, Problems 1.2
	are of particular interest in the framework of quantum mechanical inverse scattering.


Rosane Ushirobira
email: rosane.ushirobira@inria.fr

Anja Korporal
email: anja.korporal@inria.fr

Wilfrid Perruquetti
email: wilfrid.perruquetti@inria.fr

On an algebraic method for derivatives estimation and parameter estimation for partial derivatives systems

In this communication, we discuss two estimation problems dealing with partial derivatives systems. Namely, estimating partial derivatives of a multivariate noisy signal and identifying parameters of partial differential equations. The multivariate noisy signal is expressed as a truncated Taylor expression in a small time interval. An algebraic method can be then used to estimate its partial derivatives in the operational domain. The same approach applies for the parameter identification problem. The main objects in the aforementioned algebraic method are particular differential operators called annihilators. Stafford's theorem can be applied to formalize the left ideal in the Weyl algebra formed by the annihilators. In the spatial domain, iterated integrals present in the estimates yield noise filtering.

I. INTRODUCTION

Many problems in Engineering concern the estimation of partial derivatives of a noisy multivariate signal. For instance, they appear in the field of signal processing or control, or in economy issues as well.

Numerous methods were developed to tackle this numerical differentiation problem, the finite differences techniques being the most usually used.

The instability of possible solutions to these problems arise from the presence of noise due to the differentiation.

The use of algebraic tools for numerical differentiation in the multivariate case were addressed by Riachy et al. in [START_REF] Riachy | Multivariate numerical differentiation[END_REF], [START_REF] Riachy | An algebraic method for multi-dimensional derivative estimation[END_REF], [START_REF]Différenciation numérique multivariable I : estimateurs algébriques et structure[END_REF]. Their inspiration comes from the algebraic approach initiated by M. Fliess and H. Sira-Ramírez [START_REF] Fliess | An algebraic framework for linear identification[END_REF] and from the solutions proposed by Mboup et al. in [START_REF] Mboup | Numerical differentiation with annihilators in noisy environement[END_REF].

The algebraic method developed here is inspired by the parameter estimation methods elaborated in [START_REF] Ushirobira | Estimation algébrique des paramètres intrinsèques d'un signal sinusoïdal biaisé en environnement bruité[END_REF], [START_REF]Algebraic parameter estimation of a biased sinusoidal waveform signal from noisy data[END_REF], [START_REF]Algebraic parameter estimation of a multi-sinusoidal waveform signal from noisy data[END_REF]. To illustrate this approach on a multi-dimensional parameter estimation problem, we examine a simple particular case of a partial differential equation. The example of the heat conduction on a thin rod is then discussed and treated by algebraic estimations. An earlier work in this direction, also based on algebraic techniques, can be found in [START_REF] Rudolph | An algebraic approach to parameter identification in linear infinite dimensional systems[END_REF].

In Section II, the two estimation problems are briefly described. The algebraic framework for the proposed solutions is given in Section III. More details on the solutions can be found in Section IV. Computer simulations to illustrate our estimates will be presented in the final version. Let x = (x 1 , . . . , x m ) be an element in K m (m ∈ N). An mtuple N ∈ N m will be usually written as N = (N 1 , . . . , N m ). Let us consider the partial order on N m defined by

N M if N i ≤ M i , for all 1 ≤ i ≤ m.
In this paper, we wish to discuss the construction of closed-forms estimators for partial derivatives in the case of a multivariate signal f (x) with f :

U ⊂ R m → R for U some neighborhood of 0. For a given I = (i 1 , . . . , i m ) ∈ N m , we denote |I| := i 1 + • • • + i m , I! := i 1 ! . . . i m !, x I = x i 1 1 . . . x i m m and ∂ I ∂ x I = ∂ i 1 ∂ x i 1 1 . . . ∂ im ∂ x im m .
In practical problems, the available signal is corrupted by a noise. Let us denote by f n the noisy multivariate signal

f n (x) = f (x) + ϖ(x),
where ϖ(x) denotes an additive noise.

For the time being, we examine the function f only. Assume that f admits a Taylor series expansion at 0 and write:

f (x) = ∑ I∈N m a I I! x I , where a I = ∂ I f ∂ x I (0).
For N = (N 1 , . . . , N m ) ∈ N m , the truncated Taylor series f N at order N is given by:

f N (x) = ∑ I N a I I! x I . (1)
Recall that for a function g : R m → R, its multivariate Laplace transform is given by:

G(s) := L(g)(s) = R m + g(t) e -s T t dt, (2)
where T t denotes the transpose of t ∈ R m . That implies, for instance:

L x I I! = 1 s I+1 .
To realize f N (x) in the operational domain, we apply the Laplace transform (2) on the equation [START_REF] Riachy | Multivariate numerical differentiation[END_REF]. It results:

F N (s) = ∑ I N a I s I+1 , (3)
where s = (s 1 , . . . , s m ) denotes the Laplace variable and s

I = s i 1 1 . . . s i m m .
The idea presented here is to individually estimate each derivative a J for J N. To formalize our procedure, we consider the following sets:

Θ = {a I | I N}, Θ est = {a J } and Θ est = Θ \ Θ est .
The definition of Θ, Θ est and Θ est is clear: Θ contains all the parameters, Θ est contains the parameters to be estimated and Θ est the remaining ones. The relation (R) below follows from (3): R :

P(s)F N (s) + Q(s) + Q(s) = 0 (4)
where

P(s) = s N , Q(s) = -a J s N-J-1 ∈ K Θ est [s, 1 s ] and
Q(s) = -∑ I N,I =J a I s N-I-1 ∈ K Θ est s, 1 s (5)
By K Θ est and K Θ est , we denote respectively the algebraic extensions K Θ est = K(Θ est ) and K Θ est = K(Θ est ). Sharing the relation R into the polynomials P, Q and Q is justified by the fact that Q contains the coefficient to be estimated, while Q is formed by all the remaining terms. To obtain an equation containing only known terms and a J , the polynomial Q must be eliminated. That will provide a formula for the estimate of a J .

To annihilate Q, some particular differential operators must be chosen to act on (R). These operators are called annihilators. Such estimators for a J will be constructed by using structural properties of the Weyl algebra.

Let us stress that if Π is an annihilator estimating a J , the partial derivative of f at any other point p ∈ K m can be obtained by computing Π (L ( f (x + p))).

For x, τ ∈ K m , we use the notation: Recall that for a multivariate function g and its Laplace transform G, the inverse Laplace transform satisfies

L -1 1 s I ∂ J G ∂ s J = 1 (I -1)! x 0 (x -τ) I-1 (-τ) J g(τ)dτ. (6)
Setting v I,J (τ) = (xτ) I (-τ) J , a shorter notation can be used:

L -1 1 s I ∂ J G ∂ s J = 1 (I -1)! x 0 v I-1,J (τ)g(τ)dτ (7)
As we have pointed out, the resulting equation from the elimination of Q in (R) contains the desired estimate of a J and known terms. The return to the spatial domain is done by the action of the inverse Laplace transform on this equation.

B. Parameter estimation

A much alike method can be applied to identify parameters in a partial differential equation. To illustrate this assertion, we present the simple example of the one-dimensional heat equation. A similar algebraic approach was studied, for instance, in [START_REF] Rudolph | An algebraic approach to parameter identification in linear infinite dimensional systems[END_REF] for this problem and in [START_REF] Gehring | Algebraic identification of heavy rope parameters[END_REF] for the parameter identification of a linear model of the planar motion of a heavy rope. In both papers, uni-dimensional Laplace transforms are used and operational functions obtained as the solutions of an initial value problem. Contrariwise, here the equation is converted to the operational domain by using two-dimensional Laplace transform.

Consider the problem of the heat conduction in a thin rod of length 1. Let w : (z,t) → w(z,t) be the function representing the temperature at position z at time t. The partial differential equation describing this problem is given by:

∂ 2 ∂ z 2 w(z,t) -β ∂ ∂t w(z,t) -αw(z,t) = 0 (8)
The initial temperature is assumed to be constant, say w(z, 0) = 30. We denote it by f : z → w(z, 0). The boundaries are assumed to be perfectly isolated, so the heat may not flow at the boundaries. Therefore, ∂ ∂ z w(0,t) = ∂ ∂ z w(1,t) = 0. We assume also a Dirichlet boundary condition by setting v : t → w(0,t) = 30 e -t 4 . The goal is to identify the parameters α and β . The algebraic method used in the previous subsection is applied here. The Laplace transform ( 2) is employed to realize the partial differential equation ( 8) as an algebraic equation in the Laplace variable s = (s 1 , s 2 ):

s 2 1 -β s 2 -α W (s 1 , s 2 ) + β F(s 1 ) -s 1 V (s 2 ) = 0, (9)
where W (s), F(s 1 ) and V (s 2 ) denote the Laplace transforms of w(z,t), f (z) and v(t), respectively. Since

F(s 1 ) = 30 s 1 and V (s 2 ) = L (v) (s 2 ) = 120 4s 2 + 1 ,
the equation ( 9) leads to:

s 2 1 -β s 2 -α W (s 1 , s 2 ) + 30 β s 1 -120 s 1 4 s 2 + 1 = 0 (10)
Here the set of parameters Θ est to be estimated is

Θ est = {α, β },
while Θ est = / 0. Rewriting [START_REF] Gehring | Algebraic identification of heavy rope parameters[END_REF] in the form of a R-relation (see ( 4)) gives:

R : P(s)W (s) + Q(s) + Q(s) = 0 ( 11
)
where

P(s) = (4s 2 + 1) s 2 1 -β s 2 -α and (12)
Q(s) = 30 4 s 2 s 1 + 1 s 1 β ∈ K Θ est s, 1 s and
Q(s) = -120 s 1 ∈ K Θ est [s] (13)
The problem of annihilating Q is simpler than in the previous subsection and it is tackled by finding suitable Qannihilators that will lead to a system in Θ est . It is interesting to notice that in this case, the parameters to be estimated also appear in the coefficients of the polynomial P.

III. ANNIHILATORS VIA THE WEYL ALGEBRA

We have indicated in the previous Section that our aim is to annihilate the polynomials Q, see ( 5) and [START_REF] Mboup | Parameter estimation for signals described by differential equations[END_REF]. That will be done by the action of annihilators, meaning certain differential operators (or polynomials in the variables ∂ I ∂ s I ) with polynomial coefficients (or rational functions) in the variables s 1 , . . . , s m . The Weyl algebra appears naturally in this context and its structural properties will be quite useful in the choice of the annihilators.

This algebraic viewpoint is inspired by the work of M. Fliess et al. [START_REF] Fliess | An algebraic framework for linear identification[END_REF], [START_REF] Fliess | Questioning some paradigms of signal processing via concrete examples[END_REF], [START_REF] Mboup | Parameter estimation for signals described by differential equations[END_REF]. Details about the algebraic notions defined in the sequel can be found in [START_REF] Mcconnell | Noncommutative Noetherian Rings[END_REF].

A. The Weyl Algebra

We recall below some basic definitions and properties of the Weyl algebra.

Definition 1: Let m ∈ N \ {0}. The Weyl algebra A m (K) (or A m ) is the K-algebra with generators p 1 , q 1 , . . . , p m , q m and relations

[p i , q j ] = δ i j , [p i , p j ] = [q i , q j ] = 0, ∀ 1 ≤ i, j ≤ m where [•, •] is the commutator defined by [u, v] := uv -vu, for all u, v ∈ A m (K).
The Weyl algebra A m can be realized as the algebra of polynomial differential operators on the polynomial ring K[s] by setting:

p i = ∂ ∂ s i and q i = s i ו , ∀ 1 ≤ i ≤ m.
That implies that A m can be written as

A m = K[q][p] = K[s] ∂ ∂ s .
The algebra of differential operators B m (K) (or B m ) on K[s] with coefficients in the rational functions field K(s) is naturally related to A m (K). In this case, we can write

B m := K(q)[p] = K(s) ∂ ∂ s .
A basis for A m is given by q p J | I, J ∈ N m where q = q i 1 1 . . . q i m m and p = p i 1 1 . . . p i m m . An operator F ∈ A m can be written in a canonical form,

F = ∑ I,J λ IJ q I p J with λ IJ ∈ K.
Similarly, an element F ∈ B m can be written as

F = ∑ I g I (s) ∂ I ∂ s I , where g I (s) ∈ K(s).
The order of F is defined as ord(F) = max{[I| | g I (s) = 0}. This definition holds for the Weyl algebra A m as well, since A m ⊂ B m . Some useful properties of A m and B m are given by the following propositions: Proposition 2: A m is a domain. Moreover, A m is simple and Noetherian. These properties are shared by B m . In addition, A m is neither a principal right domain, nor a principal left domain. Nevertheless this is true for B 1 :

Proposition 3: B 1 admits a left division algorithm, that is, if F, G ∈ B 1 , then there exists Q, R ∈ B 1 such that F = QG+R and ord(R) < ord(G). Consequently, B 1 is a principal left domain. Alas, this proposition does not hold for B m for m ≥ 2. But a celebrated theorem by T. Stafford (see [START_REF] Stafford | Module structure of Weyl algebras[END_REF]) provides an remarkable property on the number of generators of a left ideal in the Weyl algebra. Namely, Stafford proved that every left ideal of D (D = A m or B m ) can be generated by two elements in D:

Theorem 4 (Stafford): Let a be a left ideal of D generated by three elements F 1 , F 2 and F 3 ∈ D. Then, there exist G 1 and G 2 ∈ D such that

a = D(F 1 + G 1 F 3 ) + D(F 2 + G 2 F 3 ).
An effective implementation in Maple, named Stafford, of this important theorem can be found in the work of A. Quadrat and D. Robertz [START_REF] Quadrat | Computation of bases of free modules over the weyl algebras[END_REF].

Remark 5: It is important to notice that the principality of B 1 was largely used in the initial works on algebraic methods applied to univariate numerical differentiation, such as [START_REF] Mboup | Numerical differentiation with annihilators in noisy environement[END_REF] or parameter estimation in ordinary differential equations, see for instance [START_REF]Algebraic parameter estimation of a multi-sinusoidal waveform signal from noisy data[END_REF]. In the multivariate case, the principality holds no longer, therefore the importance of Stafford's theorem.

To close this part, we remark a useful identity: Remark 6: For arbitrary N, M ∈ N r , we have

∂ N ∂ s N 1 s M = ∑ 0 J N N J (-1) |N-J| M N-J s M+N-J ∂ J ∂ s J ,
where N J = n 1 j 1 . . . n r j r , M N = m n 1 1 . . . m n r r and m n i i denotes the rising factorial (m n i i = m i (m i + 1) . . . (m i + n i -1)).

B. Annihilators

From now on, we consider m ≥ 2. In this subsection, we construct annihilators. These differential operators help to construct algebraic estimators, either of the partial derivatives, or else of parameters.

Definition 7: Let R ∈ K Θ est s, 1 s . A R-annihilator w.r.t. B m is an element of Ann B m (R) = {F ∈ B m | F (R) = 0}.
Let us remark that Ann B m (R) is a left ideal of B m . Therefore, by Stafford's theorem (Theorem 4), Ann B m (R) is generated by two generators Π 1 and Π 2 ∈ B m :

Ann B m (R) = B m Π 1 + B m Π 2 .
We call the annihilators Π 1 and Π 2 minimal R-annihilators w.r.t. B m . Notice that Ann B m (R) contains annihilators in finite integral form, i.e. operators with coefficients in

K 1 s . Lemma 8: Consider R (s) = αs N , N = (n 1 , . . . , n m ) ∈ Z m with α ∈ K Θ est . A minimal R-annihilator is given by s i ∂ ∂ s i -n i , ∀1 ≤ i ≤ m.
Recall that the degree of a monomial s I ∈ K s, 1 s is |I|. The total degree of a polynomial in s is the maximum degree of its monomials.

Remark 9: Consider R ∈ K Θ est s, 1 s with a monomial s I of maximal degree. So R has total degree |I|.

Let i k = max{i j | j = 1, . . . , m}. If |I| > 0, then ∂ i k +1 ∂ s i k +1 k is clearly an R-annihilator.
Now, recall that the polynomial to be annihilated is:

Q(s) = -∑ I N,I =J a I s N-I-1 ∈ K Θ est s, 1 s .
By the previous Remark, it results immediately: Lemma 10: The differential operators

∂ n k ∂ s n k k are Q- annihilators, for all 1 ≤ k ≤ m.
In order to construct an alternative annihilator, an algorithm is sketched below:

Algorithm 11:

Input: a polynomial R = ∑ I∈Z m finite b I s I in K Θ est s, 1 s of total degree d ∈ N Output: an R-annihilator 1) set Π = 1 ∈ D.
2) choose a monomial of degree d, say b J s J with J = ( j 1 , . . . , j m ) (so

|J| = d). 3) choose j k = min{ j > 0 | = 1, . . . , m}. 4) apply π = s k ∂ ∂ s k -j k (see Lemma 8) on R. 5) (a) if π (R) = 0, then return π and stop the algorithm. (b) if π (R) = 0, then set Π = π • Π and return to step (2) with R ← π(R). Example 12: Consider m = 2 and R(s 1 , s 2 ) = a 00 s 2 1 s 2 + a 01 s 2 1 + a 10 s 1 s 2 + a 20 s 2 ∈ K Θ est [s 1 , s 2 ]. A R-annihilator constructed with the above algorithm is s 1 ∂ ∂ s 1 -2 • s 2 ∂ ∂ s 2 -1 .
The concept of an estimator must be defined in order to take into account the remaining terms in the relation (R), after the action of an annihilator. Notice that the parameters to be estimated might appear in the set of coefficients of both polynomials Q and P, but they might be present exclusively in one of the two. Therefore a Q-annihilator must not eliminate all terms in Θ est , as formalized in the definition below:

Definition 13: An estimator π ∈ B is a Q-annihilator satisfying coeffs (π(R)) ∩ K Θ est = / 0,
where coeffs (R) denotes the set of coefficients of a polynomial R ∈ K Θ s, 1 s . It is implied by this definition that the criteria on the coefficients must be considered in the choice of annihilators in algorithm 11.

IV. DERIVATIVE ESTIMATION

For the sake of simplicity, the two-dimensional case only is presented in the sequel. In addition, to lighten the notation, the signal is written as

f (x, y) = ∑ (i, j)∈N 2 a i j i! j! x i y j , where a i j = ∂ i+ j f ∂ x i ∂ y j (0). We write also B instead of B 2 .
Assume that the parameter to be estimated is a k . For a truncated Taylor series at N = (N 1 , N 2 ), the polynomials appearing in the relation (R) are (see [START_REF] Mboup | Numerical differentiation with annihilators in noisy environement[END_REF]):

P(s 1 , s 2 ) = s N 1 1 s N 2 2 , Q(s 1 , s 2 ) = -a k s N 1 -k-1 1 s N 2 --1 2 ∈ K Θ est s, 1 s and Q(s 1 , s 2 ) = -∑ (i, j) N (i, j) =(k, ) a i j s N 1 -i-1 1 s N 2 -j-1 2 ∈ K Θ est s, 1 s .
In what follows, we provide estimates for all derivatives a i j , 0 ≤ i, j ≤ 3. Due to clear reasons of symmetry, it is enough to present the estimations of a 00 , a 10 , a 20 , a 30 , a 11 and a 21 .

For each a k to be estimated, we use a truncation of the Taylor series at N = (k, ). Notice that in this case, from [START_REF] Stafford | Module structure of Weyl algebras[END_REF] it follows:

Q(s 1 , s 2 ) = - a k s 1 s 2 .
The first step of the estimation is to determine minimal Q-annihilators. To begin, Lemma 8 helps to find two Qannihilators s 1

∂ ∂ s 1 + 1 and s 2 ∂ ∂ s 2
+ 1. However, they are not estimators since they clearly eliminate Q as well.

We follow the algorithm 11 to determine a Q-annihilator that is also a estimator.

The next table summarizes the various data about the coefficients a k :

coefficient P Q Q -annihilator a 00 1 0 1 s 1 s 2 a 10 s 1 -a 00 s 2 1 s 2 1 s 2 ∂ ∂ s 1 a 20 s 2 1 -a 00 s 1 s 2 -a 10 s 2 1 s 3 1 s 2 ∂ 2 ∂ s 2 1 a 30 s 3 1 - a 00 s 2 1 s 2 -a 10 s 1 s 2 -a 20 s 2 1 s 4 1 s 2 ∂ 3 ∂ s 3 1 a 11 s 1 s 2 -a 00 -a 10 s 1 -a 01 s 2 1 s 2 1 s 2 2 ∂ 2 ∂ s 1 s 2 a 21 s 2 1 s 2 -a 00 s 1 -a 10 -a 01 s 1 s 1 -a 11 s 2 -a 20 s 1 1 s 3 1 s 2 2 ∂ 2 ∂ s 1 s 2 s 1 ∂ ∂ s 1 -1
Let us remark that some annihilators are proposed in [START_REF] Riachy | Multivariate numerical differentiation[END_REF]. Using remark 6, it can be shown that the annihilators given above coincide with their minimal annihilators.

Remark 14: Another way to proceed would be to construct other two annihilators for each Q. Then, to apply Stafford's theorem to obtain two minimal generators by using the package Stafford. The final step is to observe the criteria in Definition 13 in these generators to obtain an estimator. Applying the inverse Laplace transform (6) provides the consequent estimates. Some summarized computations can be found below:

1) Estimating a 00 . The action of the estimator provides the resulting expression in the operational domain:

F N (s 1 , s 2 ) s 1 s 2 = a 00 s 2 1 s 2 2 .
After applying the inverse Laplace transform (6), we obtain the estimate:

a 00 = 1 xy x 0 f (τ)dτ
2) Estimating a 21 . The action of the estimator provides the resulting expression in the operational domain: After applying the inverse Laplace transform (6), we obtain the estimate:

a 21 = - 360 x 5 y 3 x 0 (v 0,0,2,1 + v 0,1,2,0 + 4v 1,0,1,1 + 4v 1,1,1,0 +v 2,1,0,1 + v 2,1,0,0 ) (τ) f (τ)dτ V. PARAMETER ESTIMATION The polynomial Q to be annihilated in the equation (12) is Q(s) = -s 1 ∈ K Θ est [s]. Some natural Q-annihilators are: π 1 = ∂ 2 ∂ s 2 1 , π 2 = ∂ ∂ s 2 and π 3 = ∂ ∂ s 2 ∂ ∂ s 1 . (14)
Although being unnecessary, the following Lemma can be checked by using Stafford:

Lemma 15: The two minimal Q-annihilators are π 1 and π 2 . Applying π 1 on the relation [START_REF] Fliess | Questioning some paradigms of signal processing via concrete examples[END_REF] gives:

2W + 4s 1 ∂ ∂ s 1 W + s 2 1 -β s 2 -α ∂ 2 ∂ s 2 1 W + 60 s 3 1 β = 0
where W = W (s 1 , s 2 ). In order to apply the inverse Laplace transform (6), we divide the above equation by s 3 1 s 2 2 . Using the notation v I,J = v I,J (τ, η) = (zτ) i 1 (tη) i 2 (-τ) j 1 (-η) j 2 , for all I = (i 1 , i 2 ), J = ( j 1 , j 2 ) ∈ N 2 , we obtain in the spatial domain: (v 2,1,0,0 -4v 1,1,1,0 + v 0,1,2,0 ) w(τ, η)dτdη.

In the same manner, an expression resulting from the action of π 2 can be found: 4s 4v 0,1,0,1 + 2v 0,2,0,0 + 1 2 v 0,2,0,1 w(τ, η)dτdη.

Solving the system provides the estimates of α and β .

VI. CONCLUSION

The goal of this communication was to examine two estimation problems under an algebraic viewpoint. The proposed algebraic method is inspired by the work of Fliess et al.

We have discussed the numerical derivation problem in the multivariate case. Basic properties of the Weyl algebra were applied to obtain the desired estimates.

A simple example of parameter estimation for the onedimensional heat equation was also presented in order to demonstrate the algebraic method.

x m 0 g

 0 (τ 1 , . . . , τ m ) dτ 1 . . . dτ m .

A 1 whereA 11

 111 11 α + A 12 β = B 0,2,0 w(τ, η)dτdη and B 1 = -(z,t) 0

  2 1 -4α -(8s 2 + 1) β W + (4s 2 + 1) s 2 1β s 2αAfter dividing the above equation by a suitable monomial in s, we obtain A 21 α + A 22 β = B 2 in the spatial domain. A system on Θ est results from the aforementioned actions of Q-annihilators:

							∂ ∂ s 2	W
							+	120 s 1	β = 0.
		A 11 A 12 A 21 A 22	α β	=	B 1 B 2
	with A 21 = -
	0	(z,t)	1 4	v 2,2,0,0 +	1 2	v 2,1,0,1 + 2v 2,0,0,1

(z,t)

0 1 4 v 2,2,0,1 + v 2,2,0,0 + 2v 2,1,0,1 w(τ, η)dτdη, A 22 = -+4v 2,1,0,0 ) w(τ, η)dτdη + 10t 2 z 3 and B 2 = -(z,t) 0


Antoine Liutkus

Zafar Rafii

Bryan Pardo

Derry Fitzgerald

Laurent Daudet

Laurent Daudet Kernel Spectrogram

Kernel Spectrogram models for source separation

Keywords: audio source separation, spatial filtering, spectrogram models

come

I. INTRODUCTION

Source separation is a field of research that gathered much attention during the last 20 years [START_REF] Comon | Handbook of Blind Source Separation: Independent Component Analysis and Blind Deconvolution[END_REF]. Its objective is to recover several unknown signals called sources that were mixed together into observable mixtures. In audio signal processing, the sources are typically understood as different auditory streams [START_REF] Bregman | Auditory Scene Analysis, The perceptual Organization of Sound[END_REF] that make sense perceptually. In music processing for instance, they correspond to different instruments playing in a song. In spoken speech enhancement, one source may be the target voice whereas others correspond to background noise to filter out.

One of the dominating paradigm today for the separation of audio waveforms is the use of generalized Wiener filtering [START_REF] Benaroya | Audio source separation with a single sensor[END_REF], [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF] under Gaussian assumptions. In practice, this approach requires good models of the spectrograms of This work is supported by LABEX WIFI (Laboratory of Excellence within the French Program "Investments for the Future") under references ANR-10-LABX-24 and ANR-10-IDEX-0001-02 PSL. Correspondance should be addressed to antoine@liutkus.net each source along with its spatial characteristics and permits very good separation provided these parameters are well estimated. The main challenge in achieving good separation then mainly becomes the devising of good spectrogram models that catch the main features of the sources to separate while requiring few parameters. Techniques such as Nonnegative Tensor Factorizations (NTF, see [START_REF] Cichocki | Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multiway Data Analysis and Blind Source Separation[END_REF]) are often used to this purpose. Their principle is to assume that the spectrogram of each source may be decomposed as the sum of only a few spectral templates activated over time. In spite of their appealing tractability, NTF models often come with some limitations. First, they often fail at efficiently decomposing in a concise way many sources such as voice that exhibit a great variety of spectra. Second, they typically assume that different sources are characterized by different sets of spectra, which may not be realistic, e.g. for mixtures of speech.

In this study, instead of decomposing the spectrograms of the sources as a combination of fixed patterns, we rather focus on their regularities to identify them from the mixtures. Auditory Scene Analysis [START_REF] Bregman | Auditory Scene Analysis, The perceptual Organization of Sound[END_REF] indeed demonstrated on perceptual grounds that apart from the commonly used harmonic property that refers to an absolute feature of many auditory sources, local features such as repetitivity, continuity or common fate are fundamental in our ability to discriminate them within a mixture. These dynamic features can be seen not to depend on any particular spectral absolute template modelled by NTF, but rather on local regularities concerning their evolution over time, frequency and space.

In order to model dependencies within the spectrograms of the sources, we use kernel local parametric models, that are deeply rooted in the local regression approach [START_REF] Cleveland | Locally weighted regression: An approach to regression analysis by local fitting[END_REF]. Basically, the value of the spectrogram of a source at some time-frequency (TF) bin is supposed to be close to its values nearby. There are often sophisticated ways to decide whether two TF bins have similar values. In full generality, proximity kernels are introduced which give the proximity of two TF points from the perspective of a source. There are several ways of building such kernels, including direct analytical expressions or through the use of feature spaces. In any case, the value of a source spectrogram is supposed to be correctly estimated using its values at locations whose proximity is high. Separation of additive sources in this context can be performed using a variant of the backfitting algorithm [START_REF] Hastie | Generalized additive models[END_REF]. Different sources are modelled through different proximity kernels. The approach, coined in as Kernel Additive Modelling (KAM), is flexible enough to permit taking prior knowledge about the dynamics of many kinds of signals into account. The proposed methodology encompasses many popular methods for audio source separation, such as DUET [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF], ADRESS [START_REF] Barry | Real-time sound source separation using azimuth discrimination and resynthesis[END_REF], REPET and REPET-SIM [START_REF] Rafii | Repeating pattern extraction technique (REPET): A simple method for music/voice separation[END_REF], [START_REF] Fitzgerald | Vocal separation using nearest neighbours and median filtering[END_REF], [START_REF] Rafii | Music/voice separation using the similarity matrix[END_REF], median filtering for drums removal [START_REF] Fitzgerald | Harmonic/percussive separation using median filtering[END_REF], etc. Moreover, it provides an efficient way to devise new specific separation algorithms for sources that are characterized by local features, rather than by a global additive model such as NTF. We show its performance on music/voice separation and provide a complete MATLAB implementation.

II. MODEL AND METHOD

II-A. Notations and model

Let the mixture x be a set of I time series, where x (n, i) denotes the value of the i th channel of the mixture at sample n. In music processing, we often have I = 2 in the stereo case. We assume that the mixture is the sum of J sources sj :

x (n, i) = J j=1 sj (n, i).
Let {s j } j=1•••J and x be the Short Term Fourier Transforms (STFTs) of the J sources and of the mixture, respectively. They are all N f × N t × I tensors, where N f is the number of frequency bands and N t the number of frames. s j (f, t) is the I × 1 vector that gives the value of the STFT s j for all channels (e.g. left and right) at TF bin (f, t).

Under the Local Gaussian Model [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF], the vectors s j (f, t) for all TF bins of a multichannel audio signal are assumed to be independent, each one of them being distributed with respect to a multivariate centered complex Gaussian distribution:

∀ (f, t) , s j (f, t) ∼ N c (0, s j (f, t) R j (f )) . (1)
In expression (1), boldfaced s j (f, t) ≥ 0 indicates the spectrogram of source j at TF bin (f, t). It is a nonnegative scalar that basically accounts for the energy of that source at TF bin (f, t). R j (f ) is a I × I positive semidefinite matrix that is called the spatial covariance matrix of source j at frequency band f . It encodes the covariance between the different channels of s j at that frequency 1 . Such a model notably encompasses the popular linear instantaneous and convolutive cases [START_REF] Comon | Handbook of Blind Source Separation: Independent Component Analysis and Blind Deconvolution[END_REF], that correspond to a rank-1 R j (f ) [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF].

Since the mixture x (f, t) is the sum of J independent random Gaussian vectors s j (f, t), it also has a Gaussian distribution. If the parameters s j and R j are known or estimated as ŝj and Rj , it can be shown that the Minimum Mean-Squared Error (MMSE) estimates ŝj of the STFTs of the sources are readily obtained through generalized spatial Wiener filtering [START_REF] Benaroya | Audio source separation with a single sensor[END_REF], [START_REF] Cemgil | Prior structures for Time-Frequency energy distributions[END_REF], [START_REF] Liutkus | Gaussian processes for underdetermined source separation[END_REF], [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF], using:

ŝj (f, t) = ŝj (f, t) Rj (f )   J j =1 ŝj (f, t) Rj (f )   -1
x (f, t) .

(2) The waveforms of the sources in the time domain are then easily obtained through inverse STFT.

II-B. Kernel constant models for spectrograms

In many source separation studies, the spectrograms s j of the sources are taken as the activation over-time of a few K spectral templates W j (f, k):

s j (f, t) = K k=1 W j (f, k) H j (k, t) , (3)
where H j gives the activation gains of these templates over time. This approach leads to the popular Nonnegative Matrix Factorization (NMF) framework [START_REF] Fitzgerald | On the use of the beta divergence for musical source separation[END_REF], [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF], [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF], [START_REF] Févotte | Notes on nonnegative tensor factorization of the spectrogram for audio source separation : statistical insights and towards self-clustering of the spatial cues[END_REF] for audio source separation.

Here, we do not assume that the spectrogram s j of a source is properly described using a parametric model such as [START_REF] Benaroya | Audio source separation with a single sensor[END_REF]. Instead, we will draw from the ideas of local regression [START_REF] Cleveland | Locally weighted regression: An approach to regression analysis by local fitting[END_REF] to model spectrograms only locally. More specifically, prior knowledge about the source comes as neighbourhoods I j (f, t), called the proximity kernel of source j, which indicates the TF points where the spectrogram has a value equal to s j (f, t):

∀ (f , t ) ∈ I j (f, t) , s j (f , t ) ≈ s j (f, t) .
In musical signals for instance, percussive elements are known to be self-similar along the frequency axis, while harmonic stable sounds are self-similar along time [START_REF] Fitzgerald | Harmonic/percussive separation using median filtering[END_REF], leading to proximity kernels that are either vertical or horizontal, as depicted in figure 1(a) and 1(b) respectively. Alternatively, source j may be known to be repetitive at period T j , so that I j (f, t) includes {(f, t + kT j )} k∈Z , as in figure 1(c).

This approach is flexible enough to take prior knowledge about the dynamics of many kinds of signals into account. It encompasses a large number of recently proposed methods for source separation [START_REF] Fitzgerald | Harmonic/percussive separation using median filtering[END_REF], [START_REF] Rafii | A simple music/voice separation method based on the extraction of the repeating musical structure[END_REF], [START_REF] Liutkus | Adaptive filtering for music/voice separation exploiting the repeating musical structure[END_REF], [START_REF] Rafii | Repeating pattern extraction technique (REPET): A simple method for music/voice separation[END_REF], [START_REF] Rafii | Music/voice separation using the similarity matrix[END_REF], [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF] and provides an efficient way to devise new specific separation algorithms for sources that are characterized by local features, instantiated by the definition of the neighbours I j (f, t) of any point (f, t), rather than by a global model such as NTF.

II-C. The kernel backfitting algorithm

Assume that s j is not observed exactly but only through a noisy observation z j whose likelihood p (z j (f, t) | s j (f, t)) is known. This likelihood accounts for the fact that even if Fig. 1. Examples of proximity kernels to account for prior knowledge about sources. (a) vertical, for percussive elements, (b) horizontal, for stable harmonic elements, (c) periodic, for repetitive elements, (d) cross-like, for smoothly varying spectrograms such as vocals. z j (f, t) is likely to be close to s j (f, t), important discrepancies may occur during iterations. This can for example be taken into account by choosing a Laplacian likelihood, leading to:

-log p (z j (f, t) | s j (f, t)) = |z j (f, t) -s j (f, t)| .
With such a likelihood and supposing that z j (f, t) are all independent, s j (f, t) may be estimated through maximum likelihood as:

ŝj (f, t) = argmin sj (f,t) (f ,t )∈Ij (f,t) |z j (f , t ) -s j (f, t)| ,
which is readily shown to be equivalent to:

ŝj (f, t) = median {z j (f , t ) | (f , t ) ∈ I j (f, t)} , (4)
so that ŝj is readily estimated through a median filtering of z j , which can be achieved in linear complexity thanks to efficient implementations found in most numerical computing libraries.

The kernel backfitting algorithm we propose for estimation of the sources spectrograms s j is strongly inspired by the original backfitting procedure proposed in the context of nonparametric additive modelling [START_REF] Hastie | Generalized additive models[END_REF], [START_REF] Hastie | Generalized additive models[END_REF]. This algorithm proceeds in an iterative fashion, where separation and reestimation of the parameters are performed alternatively. In this procedure, the spectrograms z j (f, t) of the current estimates ŝj of the sources STFT are used as noisy observations of their true value, and re-estimation of ŝj from z j is achieved through median filtering (4).

The whole procedure is summarized in algorithm 1, where • denotes conjugate transpose and tr (•) stands for the trace of a square matrix. For more details about re-estimation of the spatial covariance matrices R j (f ), the reader is referred to [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF], [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF]. Remarkably, all sources can be handled in parallel during step 3, leading to a computationally efficient technique for source separation. Typical computing time is about 5 times slower than real time on a modern desktop computer. Computational complexity furthermore scales linearly with the duration of the audio to process and the number of iterations (typically 5).

III. EVALUATION

Separating vocals from the musical background in popular music is a very challenging task that has many applications Algorithm 1 Kernel backfitting for multichannel audio source separation with locally constant spectrogram models and binary proximity kernels.

1) Input:

• Mixture STFT x (f, t) • Neighbourhoods I j (f, t) as in figure 1. • Number L of iterations 2) Initialization • l ← 1 • ∀j, ŝj (f, t) ← x (f, t) x (f, t) /IJ • R j (f ) ← I × I identity matrix 3)
Compute estimates ŝj of all sources using (2) 4) For each source j:

a) C j (f, t) ← ŝj (f, t) ŝj (f, t) b) Rj (f ) ← I T t Cj (f,t) tr(Cj(f,t)) c) z j (f, t) ← 1 I tr Rj (f ) -1 C j (f, t) d) ŝj (f, t) ← median {z j (f , t ) | (f , t ) ∈ I j (f, t)} 5)
If l < L then set l ← l + 1 and go to step 3 6) Output:

sources spectrograms ŝj and spatial covariance matrices Rj (f ) to use for filtering [START_REF] Bregman | Auditory Scene Analysis, The perceptual Organization of Sound[END_REF].

in the entertainment industry and in the automatic indexing and querying of musical databases [START_REF] Gómez | Predominant fundamental frequency estimation vs singing voice separation for the automatic transcription of accompanied flamenco singing[END_REF]. In the recent years, it has been the topic of numerous research studies and many different techniques were devised for this purpose [START_REF] Durrieu | A musically motivated midlevel representation for pitch estimation and musical audio source separation[END_REF], [START_REF] Durrieu | Musical audio source separation based on user-selected F0 track[END_REF], [START_REF] Liutkus | Adaptive filtering for music/voice separation exploiting the repeating musical structure[END_REF], [START_REF] Huang | Singing-voice separation from monaural recordings using robust principal component analysis[END_REF], [START_REF] Rafii | Music/voice separation using the similarity matrix[END_REF], [START_REF] Fitzgerald | Vocal separation using nearest neighbours and median filtering[END_REF], [START_REF] Rafii | Repeating pattern extraction technique (REPET): A simple method for music/voice separation[END_REF]. In the following, we detail and evaluate a voice/music separation procedure based on KAM.

III-A. Data and metrics

In our experiments, the processed data consists of 10 complete stereo tracks from the album The Pet Sounds by the popular band THE BEACH BOYS. This band published an extensive set of studio recordings for this album in 1997 as a commercial release2 , which includes separated vocals and background as stereo tracks. After some manual synchronization, they were mixed down so as to produce the full-length stereo mixtures to separate.

For the purpose of evaluation, all separated full-length vocals and backgrounds tracks from each technique are segmented into 10s excerpts, yielding 168 such excerpts, for which separation performance is evaluated on both the separated vocals and background music. The metrics considered are the classical Source to Distorsion Ratio (SDR) from the BSSEVAL toolkit [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF]. It is given in dB and is higher for better separations.

Since different excerpts may yield very different separation difficulties, it is known that directly averaging BSSEVAL metrics is not meaningful [START_REF] Liutkus | Informed source separation through spectrogram coding and data embedding[END_REF]. For this reason, the deltametric ∆ SDR is considered instead. It gives the difference of the performance with those obtained through oracle Wiener separation [START_REF] Vincent | Oracle estimators for the benchmarking of source separation algorithms[END_REF], which uses the true spectrograms of the sources in [START_REF] Bregman | Auditory Scene Analysis, The perceptual Organization of Sound[END_REF]. This delta-metric were shown to be more reliable for averaging over a corpus [START_REF] Liutkus | Informed source separation through spectrogram coding and data embedding[END_REF].

III-B. Techniques and parameters

For performance comparison, each full track of the corpus was separated using the techniques IMM [START_REF] Durrieu | A musically motivated midlevel representation for pitch estimation and musical audio source separation[END_REF], RPCA [START_REF] Huang | Singing-voice separation from monaural recordings using robust principal component analysis[END_REF], REPET-SIM [START_REF] Fitzgerald | Vocal separation using nearest neighbours and median filtering[END_REF], [START_REF] Rafii | Music/voice separation using the similarity matrix[END_REF], adaptive REPET [START_REF] Liutkus | Adaptive filtering for music/voice separation exploiting the repeating musical structure[END_REF] and adaptive REPET with a further DUET processing [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF].

For KAM separation, background and vocals were both modelled as having locally constant spectrograms as described in section II-B. In one so called KAM multirepet setting, the musical accompaniment is modelled as the sum of 5 repeating patterns as in figure 1 (c). In another KAM multirepet+harm setting, a further stable harmonic source is included in the background model as in figure 1 (b) and its length corresponds to 2s. In all cases, the vocal part was modelled using the cross-like kernel of figure 1 (d), whose height and length were respectively set to 50Hz and 0.4s. Frames of 90ms with an overlap of 85% were considered for the computation of STFTs and the periods of the patterns were estimated by a peak-picking of the beat-spectrum [START_REF] Rafii | A simple music/voice separation method based on the extraction of the repeating musical structure[END_REF].

Running time is approximately 5 times slower than real time and varies linearly with the number of iterations (typically 5) and the number of sources. A full MATLAB implementation of the proposed method is made publicly available on the companion webpage of this paper 3 , along with audio examples.

III-C. Results

Considering the results given on figure 2, we note that objective performance of the proposed KAM setup is at the state of the art level. More precisely, the background is consistently shown to be better estimated than with other competitive methods, which is very interesting for karaoke applications and is an encouraging result. However, the scores tend to show that other techniques such as adaptive REPET [START_REF] Liutkus | Adaptive filtering for music/voice separation exploiting the repeating musical structure[END_REF] give better estimates for the vocals. Several remarks may be done considering this evaluation.

First, the lack of a large full-tracks audio corpus prevents giving separation performance for different kinds of music. Higher is better. The background is shown to be well separated.

The Pet Sounds album is indeed characterized by mostly center-panned vocals, which fits well the assumptions of a DUET approach [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF] and explains the very high scores of the aREPET+DUET method. Extensive informal perceptual testing has shown that KAM is very robust to different kinds of music, from black metal to jazz to electro-pop music.

The reader is strongly encouraged to listen to the separated signals on the companion webpage of this paper and to run the provided MATLAB script on his own sound examples. Second, the scores given here only hold for the particular choice of proximity kernels we made in this voice/music separation task. KAM may be used in many other settings or yield improved performance with more adequate proximity kernels depending on the track considered. Remarkably, 5 out of the 7 techniques evaluated here can be understood as particular instances of KAM.

IV. CONCLUSION

In this paper, we have proposed a new framework for audio source separation, where each source is modelled through the local regularities of its spectrogram. The spectrogram taken at some time-frequency bin is supposed to be close to its values nearby, where nearness is defined through a source-specific proximity kernel. Separation is performed using a variant of the backfitting algorithm, coined in as kernel backfitting. The proposed method comes as a unifying framework for many state-of-the-art techniques for source separation and yields an easy and principled way to combine local models in order to build sophisticated mixture models. The corresponding algorithms are easy to implement and provide good performance.

Fig. 2 .

 2 Fig.2. Distribution of the ∆ SDR score over 10s excerpts. Higher is better. The background is shown to be well separated.

Thus, x (f, t) and s j (f, t) are I × 1 vectors, boldfaced s j (f, t) is a scalar and R j (f ) is a I × I matrix. Estimates are denoted ŝj , ŝj and Rj (f ).

THE BEACH BOYS, The Pet Sounds Sessions, Capitol rec. 1997.