D I Florent

Bartolo

Connected Creation: The Art of Sharing

Keywords: Peer-to-peer network, shared creation, poetics of data

come

Introduction

Art does not evolve in a separate sphere, preserved from reality. Art is transitive and interacts with social space. It establishes relations, and configures spaces. In this regard, over time, artists have used several techniques such as perspective to structure spaces of representation. For example, during the Renaissance, painters used perspective to organize "the prose of the world" [START_REF] Foucault | Les mots et les choses : Une archéologie des sciences humaines[END_REF] by placing human figures inside a unified space. At its origins (1430-1460), perspective had a political dimension and was used as a tool to build the scene of history and civic freedom (Arasse, 2006, p. 133). Perspective represents a system that differs significantly from the system of places and images that was widely used in the late Middle Ages to confer order to memory. Coined as "art of memory", this system was used from its introduction by Simonides de Ceos around 477 BC up until the 17 th century to structure memory impressions, improve recall, and contribute to the combination of ideas [START_REF] Yates | L'Art de la mémoire[END_REF]. However, we might ask what type of system do artists use to structure and arrange elements today? They use database systems within which they configure spaces and establish connections between different data sets. Online databases make it possible to retrieve data using advanced search criteria, in order to use their content to create dynamic compositions. Gregory Chatonsky, for example, has created over the last ten years several web applications such as The revolution took place in New York (2002) that use Flickr to dynamically retrieve vernacular pictures -taken by multiple individuals all over the world -in order to associate them with other images and textual statements (in a slideshow format). Database systems are also used to configure virtual worlds such as the one in which the art installation called The tunnel under the Atlantic (Maurice Benayoun, 1995) allows to dig. A database system is used to manage the visibility of a collection of pictures in a virtual world based on the time that each user spends to contemplate them: The more time a user spends in front of a picture, the more chances he/she will have of finding pictures belonging to a same category. However, database systems do not only make it possible to build dynamic compositions and virtual worlds. They also help to organize social networks: both centralized and decentralized.

In centralized networks, databases tend to be used to collect large amounts of data and to secure their access. Thereby, they allow organizations to control and monetize access to their content. Yet, databases can also serve to collect metadata in order to facilitate the exchange of data hosted by peers. Databases are often used in this manner to collect metadata associated with files that are available in peer-to-peer networks: in terms of connection and communication between peers, peer-to-peer networks offer different possibilities from centralized networks. In a peer-to-peer network, users interact directly with each other to exchange streams of data without having to first obtain permission (from a central entity). Many researchers working in different fields have used this technology over the last thirty years to communicate and collect data. For example, in 1997, the Space Sciences Laboratory at the University of California used peer-to-peer technology to launch a distributed computing system aimed at analyzing radio signals, searching for signs of extra terrestrial intelligence with the help of a large Internet-based community. This project, called SETI@home, "captured, as stated, the public's imagination and achieved tremendous popularity" (SETI@home Classic: In Memoriam, 2005). However, artists have also used peer-topeer systems to collect and exchange data. They may even have used this technology to create original networks: networks of "shared creation" that allow new types of interactions between themselves and their audience.

Launched in 2006 by the CiTu-Paragraphe, the In/Out project is one of these attempts to create a network dedicated to artistic creation using peerto-peer technology and a database. The In/Out network aims indeed to create an ecosystem of artworks by making the exchange of data between artists possible in real-time. The artistic potential of this type of network is addressed hereafter. First, the types of connections that the network allows to establish between artists are defined. Second, several artworks that have been designed especially for the In/Out network -today called Art Collider -are described. Finally, the relations that exist between a connected artwork and the network itself -considered as whole entity -are examined in order to better understand how they work together.

Establishing a connection

The In/Out network is the result of a collaborative work between different artists and programmers. Several teams were set up which were responsible for creating art installations that shared data streams in real-time using a same system of communication.

However, the different installations do not share the same type of streams. The terms of the exchanges are based on the individual artists' personal interests. Each installation effectively reflects a different artistic position and proposes a unique experience of the network thanks to different types of data streams. In this regard, the In/Out network allows not only to share video streams, but also to share audio and raw data streams. In addition, artists can use these streams of data in many ways.

In order to establish meaningful connections between the different installations, over the last years the CiTu has organized several meetings between artists and programmers. For example, ahead of the second exhibition of the network that took place in 2008, the CiTu organized a oneweek workshop to discuss the projects and their associated streams. The workshop enabled each artist to adapt their project to the needs of the others participants. The workshop also made clear that the network needed to be organized. It was therefore decided to use a server to keep track of each stream generated in order to make these accessible to everyone as soon as they would be available.

Being a peer-to-peer network, the In/Out network allows direct connection between peers, but in order to provide access to all the streams generated, the network also keeps track of all the individual streams. Each data stream of data is indeed indexed in an online database. When queried, the database releases information about all data streams currently available. None of these are actually stored. They are only indexed by the database that collects and shares metadata such as their name and their URL, that is to say metadata that can be used to identify them. However, the database is also used to store metadata other than name and URL so as to make it possible for users to search for streams matching specific criteria. Thus, for each stream, the database records a description, a screen capture, and a list of tags. The database also provides the name of the person at the origin of each stream, allowing, for example, to retrieve dynamically only the data streams produced by a specific person or tagged as "bright".

It is not necessary to query the database in order to establish a connection with another artist. Anyone can directly ask for the data stream produced by another art installation provided they know its URL. However, the database plays a key role by enabling artists to search in real-time for streams of data that match advanced search criteria, so that they do not have to rely on chance to obtain a stream that corresponds to their needs. In this respect, the database offers the possibility to create connections that are both dynamic and meaningful. Moreover, the qualities of the connections that can be established also depend on the number of streams available at the time that a request is initiated. Indeed, the network needs to reach a critical mass, that is to say a sufficient number of peers who generate their own stream of data. To reach this stage, the CiTu has launched a wiki platform that provides several software programs such as the Flomixer: The Flomixer enables anyone to join the network by providing a video stream that can be the result of a compositing process (Getting Started, 2011).

The In/Out network uses the database in a manner similar to the way that Internet services such as Napster have utilized their own database. Developed by Shawn Fanning in 1999 to offer the possibility for all to share their MP3 collections online, the Napster platform also used a database to keep track of each file available on the hard drives of its users in order to facilitate their access. By making copyrighted materials available for free, Napster defied -until its closure in 2001the rules of the market economy. In the same time, the Napster platform contributed to the creation of online communities based on peers' interests, allowing a new relational sensibility to develop among young people. The In/Out network shares these aspects with the Napster platform. However, unlike Napster, the In/Out network is not used to share copyrighted materials. It serves to coordinate the circulation of data streams generated in-real time by artistic installations.

Connected artworks

To further examine potential relations between different art installations belonging to the In/Out network during an exhibition, two specific examples are presented. The first one is an installation called Driving that we created ourselves in 2008. Driving is a video installation that shows a pony car driving through the night. The car (a 1970 Dodge Challenger R/T) was modellized and rendered using 3D software. The frames displayed by the installation are completely computer generated and are projected onto one of the four walls of a small and dark projection room. Yet, unlike frames projected in a classical movie theatre, the frames projected by Driving respond to an external signal. This signal is a data stream generated by another art installation (currently connected to the In/Out network). Indeed, the Driving installation uses an external audio stream to control the speed of the car displayed on the screen: This higher the sound level, the faster the car moves and conversely the lower the sound level, the slower the speed of the car will be. The car may even come to a stop if the stream connection is lost and no other audio stream is available. The installation does not show a movie with a beginning and an end. Instead, it plays the same dozen of video sequences in a loop, showing a car travelling on an endless road. The audio stream retrieved by the installation is also used to produce car engine sounds. Created mainly by altering the pitch of the audio stream, the sound of the car engine confers to the Dodge Challenger a real presence inside the projection room.

The Driving installation only shows outside views of the pony car. Visitors can not see the inside: they are invited to sit on a bench in front of the screen to watch the car drive by closely at different speeds. Moreover, in order to establish a more direct contact between the visitor and the network, the audio stream retrieved by the installation is also used to flash the car lights. By doing so, the installation aims to reveal the presence of the network behind the image, enabling visitors to imagine all the streams of data that are being generated and retrieved simultaneously by the different installations that are part of the network.

In addition, in order to retrieve an audio stream, the installation generates its own stream of data. The installation releases the video frames that are being projected inside the projection room in realtime, thereby offering others artists the possibility to use them as raw material. As an example, the video frames made available by Driving could be used by the computer program Crystallizer. Born from the collaboration between Robin Gareus and Joseph Nechvatal in 2009, Crystallizer is a program that attacks video frames by using colonies of viruses. The program is inspired by the work of mathematician John Horton Conway on artificial life (particularly Conway's Game of Life). It uses video frames generated by installations within the network as a host environment for viruses that feed on their colors, and multiply until they have completely altered their code and original appearance. Crystallizer is one of the programs that Joseph Nechvatal has used throughout his entire career to influence the outcome of an image thanks to autonomous agents [START_REF] Sikora | Le projet Computer Virus[END_REF]). Yet, unlike other computer virus programs that Joseph Nechvatal has used previously in his viral works, the Crystallizer infects an entire network. The attacks are not limited to a collection of pictures previously selected by the artist. They contaminate an entire living network by producing and proposing to the other peers, transformed and degraded streams of their own data.

The In/Out network is characterized by a continuous displacement of data from one peer to another, allowing artworks to invest a space without barriers. Thanks to the network artworks are not confined to an enclosed space, they can be distributed. This changes most of the rules that apply when organizing an exhibition. As Geert Lovink has stated, new media are increasingly distributed and require a new aesthetics that moves beyond the concepts of form and medium. They require namely "a distributed aesthetics" that takes into account simultaneously "the dispersed and the situated, the asynchronous production and multi-user access to artefacts on the one hand, and the highly individuated and dispensed allotment of information, on the other" (Munster & Lovink, 2005).

Distributed creation

In recent years, the In/Out network has been presented on multiples occasions. In 2008 and 2010 at La Bellevilloise in Paris for example. Yet, the exhibition of the network does not generally occur at one place. The exhibition typically takes place simultaneously at several locations in order to reflect the distributed architecture of the network. For instance, the artworks selected for the In/Out x.0 exhibition (that took place in 2008) were distributed between several exhibition venues located in Paris (like La Bellevilloise) but also in Montreuil at the premises of PétaHertz studio and at the Maison Populaire. An interactive installation connected to the network and created by Yves-Marie L'Hour and Benoit Meudic was even exhibited at the Cube (Issy-lès-Moulineaux).

The distance that separates each exhibition space from the others, confers a physical dimension to the network. Each exhibition space helps to define the presence of the network in the city by covering different territories. The exhibition spaces represent points of connection. They enable visitors to enter the network, to follow the path taken by data to travel from one place to another. Yet, inside each exhibition space, how are the connections that are being made between the different installations that belong to the network made visible?

In order to make the network visible to each visitor, the CiTu has developed a monitoring system. This system uses the information collected by the database described previously, to show in real-time all the connections that are being established or closed between peers. The monitoring system has a structuring effect on how visitors understand the structure and dynamics of the network. It shows what kinds of streams of data are being shared and where each installation is located, enabling viewers to choose between the different exhibition spaces. Presented in all exhibition spaces, the monitoring system proposes several views of the network called "Connections", "Physical and Metaphysical Locations", "Sociability and Nervousness", "Treeview and Influences", "Affinity Groups" and "Flux and Reflux". Each view offers a different representation of the network making it possible for the visitors to observe many of its characteristics. The "Physical and Metaphysical Locations" view for example places all the selected installations inside circles that represent all the exhibition spaces. This view thus allows visitors to observe the data exchanges that are being made between the different exhibition spaces. On the contrary, the "Connections" view does not take into account the location of the installations. By representing all the streams of data that are being exchanged between the different art installations it draws the viewers' attention to the connections established directly between artists. The view even specifies the type of streams of data that are being transmitted (raw data, audio or video) by using different colors to display them: Each stream is represented by a colored line that connects two installations (themselves represented as cells).

The monitoring system helps to map the network by drawing a landscape of connections. However, the map is not the network. A network is indeed opaque, ubiquitous, and non-formal. It cannot be reduced to its representation. It must be apprehended within the complex ecology in which it is developing, that is to say not merely with the monitoring system that only makes certain types of connections visible. The smooth and seamless lines that the monitoring system displays cannot account for the human efforts required to create links between each physical installation or to maintain them. The views of the monitoring system work only in relation to the installations that are being presented in each exhibition space: The monitoring system is in charge of the global representation of the network, whereas the individual installations allow to interact directly with them, to observe how data is being moved and transformed. The poetics of data that the network carries is not fully displayed by the monitoring system. It has to be experienced by entering the network, that is to say by interacting with the exhibited installations, and visiting several exhibition spaces. Unlike the monitoring system, the installations do not offer access to a global representation of the network. They allow to capture what is not yet described, what is not yet visualized, as illustrated with the flashing lights in the Driving installation. The installations also enable viewers to observe how the data streams are progressively transformed as they move from one installation to the next, in an endless relay.

Conclusion

The In/Out network is not only a tool. At a time where artistic presence over the Internet tends to be standardized, it is a collective experiment that proposes an original way of being present as a connected artist on Internet but also in an exhibition space. The In/Out network makes it possible to imagine new types of collaborations and external participations in favor of a poetics of collaboration. The In/Out network is an open space with no barriers that invites one to be part of its making, part of its development and change. The network does more than simply manage and map data. Potentially it enables anyone to play a role in its organization (structure) by proposing a stream of data or by developing a computer program that uses and transforms data streams shared by other peers. However, like any network, In/Out does not constantly grow. The life of the In/Out network is organized around exhibition projects that bring together groups of collaborating artists. The In/out network is a continuous emergent project characterized by a distributed aesthetics that invites one to embrace the surrounding flow of data. For all these reasons, the In/Out network can be considered as an artistic proposition in itself, a collective artwork. Indeed, the In/Out network does not only connect artworks, it also embodies an artwork that has to be entered and experienced through multiple views, as a visitor or as a participant.

Acknowledgements

We wish to extend our thanks to all the members of CiTu-Paragraphe who have participated in the activities of the In/Out network. We are especially grateful for the support of Julien de la Fuente and Robin Gareus successively in charge of the network development under the artistic direction of Maurice Benayoun. Our thanks go also to all the artists who have collaborated with us.


Marcus J Grote

Marie Kray

Uri Nahum

Adaptive Eigenspace Method for Inverse Scattering Problems in the Frequency Domain

Keywords: Inverse medium problem, Helmholtz equation, full waveform inversion, PDE constrained optimization, frequency stepping, parameter estimation, adaptive eigenspace inversion, regularization, sample averaging approximation. Ns

A nonlinear optimization method is proposed for the solution of inverse scattering problems in the frequency domain, when the scattered field is governed by the Helmholtz equation. The time-harmonic inverse medium problem is formulated as a PDE-constrained optimization problem and solved by an inexact truncated Newtontype iteration. Instead of a grid-based discrete representation, the unknown wave speed is projected to a particular finite-dimensional basis of eigenfunctions, which is iteratively adapted during the optimization. Truncating the adaptive eigenspace (AE) basis at a (small and slowly increasing) finite number of eigenfunctions effectively introduces regularization into the inversion and thus avoids the need for standard Tikhonov-type regularization. Both analytical and numerical evidence underpins the accuracy of the AE representation. Numerical experiments demonstrate the efficiency and robustness to missing or noisy data of the resulting adaptive eigenspace inversion method.

Introduction

Inverse scattering problems occur in a wide range of applications such as radar and sonar technology, non-destructive testing, geophysical exploration or medical imaging. By illuminating an unknown body, the scatterer, with waves of various directions or wavelengths, one attempts to obtain information about that body from the scattered waves recorded at some distance. In the inverse medium problem, the scatterer is a penetrable, bounded inhomohogeneity inside the medium characterized by one or several varying physical parameters and the inverse problem consists in estimating these parameters from scattering data. Typical inverse medium problems include oil and gas exploration [START_REF] Van Doren | Parameter identification in large-scale models for oil and gas production[END_REF] in geophysics or breast tumor detection [START_REF] Mojabi | Microwave biomedical imaging using the multiplicative regularized Gauss-Newton inversion[END_REF] in medical imaging. Numerical methods for the solution of inverse scattering problems essentially fall into either of two classes: qualitative and quantitative methods.

Qualitative methods [START_REF] Cakoni | Qualitative methods in inverse scattering theory: an introduction[END_REF] generally require little a priori knowledge about the scatterer and permit to estimate its location and shape quite efficiently. Examples are the MU-SIC (MUltiple SIgnal Classification) algorithm [START_REF] Cheney | The linear sampling method and the MUSIC algorithm[END_REF][START_REF] Kirsch | The MUSIC-algorithm and the factorization method in inverse scattering theory for inhomogeneous media[END_REF], the DORT (Decomposition of the Time Reversal Operator) method [START_REF] Prada | Decomposition of the Time Reversal Operator: Application to detection and selective focusing on two scatterers[END_REF][START_REF] Cassier | Multiple scattering of acoustic waves by small soundsoft obstacles in two dimensions: mathematical justification of the Foldy-Lax model[END_REF], the linear sampling method [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF][START_REF] Haddar | An improved time domain linear sampling method for Robin and Neumann obstacles[END_REF], the probe method [START_REF] Potthast | A survey on sampling and probe methods for inverse problems[END_REF], and the factorization method [START_REF] Kirsch | The MUSIC-algorithm and the factorization method in inverse scattering theory for inhomogeneous media[END_REF][START_REF] Lechleiter | The factorization method is independent of transmission eigenvalues[END_REF]. These methods provide a criterion for deciding whether any point lies inside or outside the scatterer only on the basis of farfield measurements of the scattered field. Hence, they effectively determine the support of the scatterer and thus permit to quickly detect material defects or obstacles, but do not provide any further quantitative information about its physical characteristics such as the local sound speed.

Quantitative methods typically reformulate the inverse scattering problem as a PDEconstrained optimization problem, where the unknown physical parameters, the control u, are determined by minimizing an appropriate objective functional [START_REF] Tarantola | Inverse Problem Theory And Methods For Model Parameter Estimation[END_REF]. The objective functional, L[y, u], measures the misfit between the simulated wave field, the state y, and the true scattered data, the observations y obs . Both optimize-then-discretize or discretizethen-optimize strategies are common [START_REF] Haber | On optimization techniques for solving nonlinear inverse problems[END_REF][START_REF] Hinze | Optimization with PDE Constraints[END_REF]. To tackle the ill-posedness of the inverse problem, a penalization term is usually included in the objective functional, such as standard Tikhonov regularization [START_REF] Vogel | Computational Methods for Inverse Problems[END_REF][START_REF] Beilina | Approximate Global Convergence and Adaptivity for Coefficient Inverse Problems[END_REF].

At least two strategies are then available: the full-space approach, where u and y are sought simultaneously, and the reduced-space approach, where y is eliminated from the objective functional as L[u] = L[y(u), u], taking advantage of the linearity of the underlying wave equation [START_REF] Haber | On optimization techniques for solving nonlinear inverse problems[END_REF][START_REF] Métivier | Full waveform inversion and the truncated Newton method[END_REF][START_REF] Grote | Inexact interior-point method for PDE-constrained nonlinear optimization[END_REF][START_REF] Van Leeuwen | 3D frequency-domain seismic inversion with controlled sloppiness[END_REF]. Then, standard Newton or quasi-Newton methods from nonlinear optimization can be applied [START_REF] Nocedal | Numerical Optimization[END_REF][START_REF] Herzog | Algorithms for PDE-constrained optimization[END_REF]. For inverse medium problems, however, the exact solution of the Newton equations at each iteration may be prohibitively expensive due to the very large number of (unknown) parameters. In recent years, inexact truncated Newton methods [START_REF] Eisenstat | Choosing the forcing terms in an inexact Newton method[END_REF][START_REF] Dembo | Truncated-Newton algorithms for large-scale unconstrained optimization[END_REF][START_REF] Nash | A survey of truncated-Newton methods[END_REF], where at each (outer) iteration the (quasi-)Newton equations are solved by using only a few (inner) Krylov subspace iterations, have proved particularly effective for large-scale inverse medium problems [START_REF] Métivier | Full waveform inversion and the truncated Newton method[END_REF][START_REF] Rieder | Inexact Newton regularization using conjugate gradients as inner iteration[END_REF].

The spatial discretization of the parameter u on a standard finite difference (FD) or finite element (FE) mesh typically results in a huge number of (unknown) nodal values, in fact increasingly so, as the mesh is refined. For time-dependent scattering problems, the adaptive eigenspace inversion (AEI) method [START_REF] De Buhan | Logarithmic stability in determination of a 3D viscoelastic coefficient and a numerical example[END_REF][START_REF] De Buhan | A new approach to solve the inverse scattering problem for waves: combining the TRAC and the Adaptive Inversion methods[END_REF] recently achieved a significant reduction in the number of parameters by projecting u to a finite-dimensional basis of eigenfunctions of a particular elliptic differential operator, which is iteratively adapted during the nonlinear optimization.

Here, we propose the AEI method for the inverse medium problem in the frequency domain, where the scattered field is governed by the Helmholtz equation. In Section 2, we describe our AEI approach, which combines state-of-the-art techniques from large-scale nonlinear optimization, such as inexact truncated Newton-like methods and frequency stepping [START_REF] Chen | Inverse scattering via Heisenberg's uncertainty principle[END_REF][START_REF] Bao | Imaging of local surface displacement on an infinite ground plane: The multiple frequency case[END_REF][START_REF] Métivier | Full waveform inversion and the truncated Newton method[END_REF], with an adaptive eigenspace representation of u for regularization. Next, in Section 3, we present both analytical and numerical evidence which underpins the remarkable accuracy of our particular choice of basis functions. In particular, we show how adapting the dimension of the eigenspace basis effectively builds regularization into the inversion. Finally, we subject our AEI method in Section 4 to a series of numerical tests that demonstrate not only its accuracy and robustness with respect to missing or noisy data, but also its versatility by combining it with a sample averaging approximation [START_REF] Haber | An effective method for parameter estimation with PDE constraints with multiple right-hand sides[END_REF].

Adaptive Eigenspace Inversion (AEI)

In this section, we present the Adaptive Eigenspace Inversion (AEI) method for inverse medium problems in the frequency domain, when the scattered wave field is governed by the Helmholtz equation. First, we formulate the inverse medium problem for the squared wave speed u(x) as a PDE-constrained optimization problem. Next, we introduce the adaptive eigenspace representation for u, the distinguishing feature of the AEI method. Finally we describe in detail the full AEI algorithm where both the dimension and the basis functions themselves are adapted during the optimization process.

Inverse medium problem

We consider a time-harmonic scattering problem in unbounded space from a penetrable inhomohogeneity located inside a bounded convex domain Ω ⊂ R d , d = 1, 2, 3. Outside Ω, the wave speed is known and may vary. Inside Ω, the scattered field y satisfies the Helmholtz equation together with a Sommerfeld-type radiation condition at the boundary Γ = ∂Ω, for simplicity:

   -ω 2 y -∇ • (c 2 (x)∇y) = f (x),
in Ω, ∂y ∂n -ik(x)y = 0, on Γ.

(

) 1
Here ω > 0 denotes the constant frequency, c(x) > 0 the wave speed and k(x) = ω/c(x) the wave number. Next, we perform N s illuminations of the medium inside Ω with source terms f = f , = 1, ..., N s , and denote by y the corresponding (unique) solutions of (1). Given the measurements y obs on Γ, or part of it, we seek to reconstruct the (unknown) squared wave speed u = c 2 inside Ω such that every solution y of (1) with f = f coincides at Γ with the measurements y obs , = 1, ..., N s . In doing so, we assume that the wave speed c is known on the boundary Γ.

To solve the inverse medium problem, we now formulate it as a PDE-constrained optimization problem and thus seek a minimizer u of the standard data misfit functional

L[u] = 1 2
Then, the inverse problem reads:

Find u ∈ V = L ∞ (Ω) such that u = argmin v∈V L[v],
with y (u) satisfying (1) for f = f , = 1, ..., N s .

(

) 3
Since the inverse problem (3) is generally ill-posed, Tikhonov regularization term is typically added for stability. Instead, we shall incorporate regularization by restricting u to a finite-dimensional subspace V K , adaptively determined during the optimization.

To solve (3), we shall consider standard Newton or Quasi-Newton methods, which require the gradient of L[u]. For a direction p, it is given by

∇L[u], p = Ns =1 Ω (∇y (u) • ∇z (u)) p , (4)
where each Lagrange multiplier z , = 1, ..., N s , solves the adjoint problem

   -ω 2 z -∇ • (u(x)∇z ) = 0, in Ω, ∂z ∂n + ik(x)z = (y -y obs ), on Γ. (5)
For the solution of (3), we can choose between the optimize-then-discretize or the discretize-then-optimize approaches. In the former case, we solve the inverse problem (2)-(4), whereas in the latter case, we consider the discrete functional

L h [u] = 1 2 Ns =1 P y (u) -y obs 2 2 , (6)
together with its corresponding discrete gradient

∇L h [u] = Ns =1 ∂y ∂u P P y (u) -y obs , (7)
where P denotes a projection matrix from the underlying discrete computational subspace to the observations. We will consider either approach indifferently and in doing so demonstrate that our approach does not rely on any particular choice of discretization.

Adaptive eigenspace basis

Instead of a standard nodal basis (FD, FE grid-based discretization), we shall use a basis of (global) eigenfunctions {φ m } m≥1 to represent the parameter u(x) as

u(x) = u 0 (x) + m≥1 β m φ m (x). ( 8
)
Here the "background" u 0 (x) ∈ H 1 (Ω) solves the elliptic problem:

-∇ • (µ(x)∇u 0 (x)) = 0, ∀ x ∈ Ω, u 0 (x) = c 2 (x), ∀ x ∈ Γ, (9)
where µ(x) is defined by

µ(x) = 1 max{|∇u(x)|, ε} , ∀ x ∈ Ω, ε > 0. ( 10
)
The parameter ε > 0, which ensures that the denominator of µ does not vanish, is typically set to a very small value such as ε = 10 -6 . While the primary role of u 0 is to accommodate the (known) inhomogeneous boundary values of u, we shall demonstrate in Section 3, that u 0 also captures much of the behavior of u in the interior.

Following [START_REF] De Buhan | A new approach to solve the inverse scattering problem for waves: combining the TRAC and the Adaptive Inversion methods[END_REF], we choose for the functions φ m the orthonormal basis of eigenfunctions φ m ∈ H 1 0 (Ω) of the elliptic operator,

-∇ • (µ(x)∇φ m (x)) = λ m φ m (x), ∀ x ∈ Ω, φ m (x) = 0, ∀ x ∈ Γ, (11)
with corresponding eigenvalues 0 < λ m ≤ λ m+1 , m ≥ 1. Clearly, at higher λ m , the eigenfunctions φ m in (11) will be increasingly oscillatory.

In Section 3, we shall provide analytical and numerical evidence which underpins the remarkable accuracy of this basis for representing any given u(x). In our AEI approach, the eigenfunctions {φ m } m and the background u 0 are repeatedly recomputed as the underlying control u(x) varies during the optimization. Hence, we call {u 0 } ∪ {φ m } m≥1 an adapted eigenspace (AE) basis.

Since u(x) is precisely the quantity we seek, and thus unknown, we always use in [START_REF] Dembo | Truncated-Newton algorithms for large-scale unconstrained optimization[END_REF] the value from the previous optimization step. At the first step, when no information about u inside Ω is available yet, we simply set µ(x) ≡ 1. Then, u 0 is a harmonic prolongation of c 2 from Γ into Ω while the basis {φ m } m≥1 simply corresponds to the eigenfunctions of the Laplacian operator in Ω.

Remark 1. The elliptic operator in [START_REF] De Buhan | Logarithmic stability in determination of a 3D viscoelastic coefficient and a numerical example[END_REF] and [START_REF] Eisenstat | Choosing the forcing terms in an inexact Newton method[END_REF] essentially coincides with the gradient of the penalized total variation (TV) regularization term [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Vogel | Iterative methods for total variation denoising[END_REF],

R T V (u) = 1 2 Ω |∇u| 2 + ε 2 dx,
given by

∇R T V (u) = -∇ • 1 |∇u| 2 + ε 2 ∇u -∇ • (µ(x)∇u) .
Penalized TV-regularization is well-known in image processing for noise removal, while preserving sharp interfaces. Instead of adding a Tikhonov regularization term to the objective functional, the AEI approach projects u to the basis of eigenfunctions of the gradient of the penalized TV-regularization functional; hence, the AE basis inherits similar properties.

AEI algorithm

In practice, we truncate the infinite sum in (8) at a finite number K ≥ 1:

u(x) = u 0 (x) + K m=1 β m φ m (x). ( 12
)
To keep both the memory requirements and the computational effort low, it is imperative to keep the number K of eigenfunctions minimal. The truncation of the eigenfunction expansion is also crucial for numerical stability, as it builds regularization into the AEI approach -see Remark 1 above but also Remark 3 below. At higher frequencies, waves detect and carry more detailed information about the scatterer, yet the number of local minima of L[u] may also increase. To minimize the chance of landing in a (false) local minimum, we also apply a standard frequency continuation procedure [START_REF] Chen | Inverse scattering via Heisenberg's uncertainty principle[END_REF][START_REF] Bao | Imaging of local surface displacement on an infinite ground plane: The multiple frequency case[END_REF]. First, we solve the inverse problem (3) at the lowest frequency ω 1 . Then we progressively increase ω = ω 2 , . . . , ω n while re-initializing the optimization at every ω j from the previous lower frequency ω j-1 . In doing so, we assume that the measurements are available through a range of frequencies, for instance via Fourier transform of a timedependent signal.

AEI Algorithm.

Input: initial guess u = 1, observations y obs . Output: u * .

1. Choose K ≥ 1 and compute {φ m } K m=1 from [START_REF] Eisenstat | Choosing the forcing terms in an inexact Newton method[END_REF] and u 0 from (9) with µ ≡ 1

2. Expand u(x) = u 0 (x) + K m=1 β m φ m (x) 3. For ω = ω 1 , . . . , ω n (a) Compute L[u] and ∇L[u], set H (b) STOP: if ∇L[u] ≤ Tol i. Solve Hp = -∇L[u] ii. Determine step size α and set u := u + αp iii. Update L[u], ∇L[u] and H (c) Set µ from (10) with ∇u (d) Update K, compute {φ m } K m=1 from (11) and u 0 from (9) (e) Expand u(x) = u 0 (x) + K m=1 β m φ m (x) 4. u * = u
The AEI approach applies regardless of the underlying optimization method used. Here we consider truncated Newton-like methods [START_REF] Dembo | Truncated-Newton algorithms for large-scale unconstrained optimization[END_REF][START_REF] Nash | A survey of truncated-Newton methods[END_REF] and denote by H either the true Hessian or some approximation of it, depending on (Newton, BFGS, or Gauss-Newton) method used. In all cases the linear system in 3(b) is solved by a truncated CG-iteration with the Eisenstat-Walker criterion [START_REF] Eisenstat | Choosing the forcing terms in an inexact Newton method[END_REF]. In 3(b)ii, the step size α of the search direction p is determined either by Armijo (Newton, Gauss-Newton) or Wolfe-Powell (BFGS) stepsize control, depending on the underlying method [START_REF] Nocedal | Numerical Optimization[END_REF]. In steps 1 and 3(d) of the AEI Algorithm, we compute the first K eigenfunctions in (11) by using a standard restarted Lanczos iteration [START_REF] Lehoucq | Deflation techniques for an implicitly re-started arnoldi iteration[END_REF].

Approximation properties of the AEI expansion

In the AEI method presented in Section 2, the (unknown) parameter u is expanded as in [START_REF] Grote | Inexact interior-point method for PDE-constrained nonlinear optimization[END_REF] in the L 2 -orthogonal basis of eigenfunctions {φ m } m=1,...,K defined by [START_REF] Eisenstat | Choosing the forcing terms in an inexact Newton method[END_REF] together with u 0 defined by [START_REF] De Buhan | Logarithmic stability in determination of a 3D viscoelastic coefficient and a numerical example[END_REF]. In this entire section, we shall assume that u is known and shall now provide some analytical and some numerical evidence which underpins the remarkable accuracy provided by our particular choice of u 0 and the AE basis.

One-dimensional case

In one space dimension, (9) reduces to

   - d dx µ(x) d dx u 0 (x) = 0 ∀ x ∈ (a, b), u 0 (a) = c 2 (a) , u 0 (b) = c 2 (b), (13)
where

µ(x) = 1 max{|u (x)|, ε} , ∀x ∈ (a, b). ( 14
)
The following result is immediate.

Proposition 1. Let u : Ω → R, Ω = (a, b), a < b, u continuous, piecewise differentiable and with C ≥ u (x) ≥ ε > 0 or -C ≤ u (x) ≤ -ε < 0, at every differentiability point x ∈ Ω, u(a) = c 2 (a) and u(b) = c 2 (b). If u 0 solves (13) then u 0 (x) = u(x), ∀x ∈ Ω.
Proof. We first assume that u is strictly increasing and therefore u (x) ≥ ε > 0. Since µ(x) = 1/ max{|u (x)|, ε} = 1/u (x), u itself also satisfies [START_REF] Haber | On optimization techniques for solving nonlinear inverse problems[END_REF] with u 0 (a) = c 2 (a) and u 0 (b) = c 2 (b). As µ is strictly positive through Ω, the (weak) solution of ( 13) is unique and therefore u 0 (x) = u(x), ∀x ∈ Ω. If u is strictly decreasing, the proof is similar for µ(x) = -1/u (x) > 0.

Remark 2. From Proposition 1 we conclude that our choice for u 0 is in fact optimal in one space dimension when u is strictly positive (or negative) throughout Ω, since it automatically yields u itself. The same conclusion immediately holds in higher dimensions, if the background medium is layered, that is if u(x 1 , x 2 , . . . , x n ) depends on a single variable.

In general, however, u 0 will not equal u, in particular when u is not monotonic. First, we illustrate the usefulness of our particular choice for u 0 in ( 13)-( 14) by comparing it to a straightforward choice where u 0 is harmonic and satisfies (13) with µ(x) ≡ 1. If u is strictly increasing (or decreasing) over Ω = (0, 1), u 0 defined by ( 13)-( 14) automatically coincides with u -see Proposition 1. Clearly, if u 0 is merely harmonic over Ω, it does not coincide with u, as shown in Fig. 1. If u both increases and decreases over Ω, neither the harmonic nor the adaptive u 0 will coincide with u. However, the adaptively computed u 0 will in general better approximate u inside Ω, as shown in Fig. 1; in fact, the adaptive u 0 coincides with u over [0.6, 1].

Next, we approximate the difference u -u 0 in span{φ 1 , φ 2 , . . . , φ K }, where all φ m satisfy (11) in one space dimension, that is

   - d dx µ(x) d dx φ m (x) = λ m φ m (x) ∀ x ∈ (a, b), φ m (a) = 0 , φ m (b) = 0, (15)
with µ defined in [START_REF] Haber | An effective method for parameter estimation with PDE constraints with multiple right-hand sides[END_REF].

Clearly, the local behavior of φ m in the neighborhood of any fixed x 0 ∈ Ω is essentially determined by the magnitude of

C = |u (x 0 )|. If |u (x)| C > ε in a neighborhood of x 0 , φ m essentially behaves like φ m (x) A m sin( Cλ m x) + B m cos( Cλ m x)
near x 0 and hence is strongly oscillatory at higher m. However, if u is essentially constant near x 0 , |u (x)| 0, then µ = 1/ε there and φ m essentially behaves like

φ m (x) A m sin( ελ m x) + B m cos( ελ m x).
Since ε is very small, φ m will also remain essentially constant there for moderate values of m. To illustrate this behavior, we now consider the piecewise linear profile u(x) shown in Fig. 2, where u (x) = 3, 0 and -9 on the interval (0, 0.3), (0.3, 0.9) and (0.9, 1), respectively. Since u is piecewise constant, we can immediately determine all eigenfunctions as:

φ m (x) =      A m,1 sin( √ 3λ m x) + B m,1 cos( √ 3λ m x), x ∈ [0, 0.3), A m,2 sin( √ ελ m x) + B m,2 cos( √ ελ m x), x ∈ [0.3, 0.9), A m,3 sin( √ 9λ m x) + B m,3 cos( √ 9λ m x), x ∈ [0.9, 1]. (16)
In Fig. 2, we show u together with some of the eigenfunctions from [START_REF] Haddar | An improved time domain linear sampling method for Robin and Neumann obstacles[END_REF]. On every subinterval, φ m has a different frequency determined by the local value of √ Cλ m . In [0.3, 0.9), the frequency ελ m is very small, as ε = 10 -4 , and φ m appears essentially constant. As λ m further increases, the frequency √ ελ m increases as well and oscillations appear. Clearly, the smaller ε, the more eigenfunctions φ m essentially behave as constants wherever u is essentially constant. We remark that φ 1 nearly coincides with u up to a scaling factor. 15) either with µ ≡ 1 or µ as in [START_REF] Haber | An effective method for parameter estimation with PDE constraints with multiple right-hand sides[END_REF]. Again, we observe that the adaptively determined eigenfunction φ 1 captures well the main features of u up to a scaling factor, unlike the first eigenfunction of the Laplacian. To reach a relative error below 1% for the best approximation with respect to the L 2 -norm, only the first six eigenfunctions φ 1 , . . . , φ 6 are needed in the adaptive case. In contrast, the first six eigenfunctions of the Laplacian yield a seven times larger relative L 2 -error, which drops below 1% only once thirteen eigenfunctions are included in the approximation; hence, the AE basis better captures the essential information about u.

Two-dimensional case

To illustrate the remarkable approximation properties of the AE basis in two space dimensions, we now consider the profile u(x), x = (x 1 , x 2 ), shown in Fig. 4. Next, we compute u 0 from ( 9) with ε = 10 -6 and µ as in [START_REF] Dembo | Truncated-Newton algorithms for large-scale unconstrained optimization[END_REF]. In Fig. 5, we observe that u 0 matches with remarkable accuracy the background medium but misses the embedded kite-shaped obstacle. The first eigenfunction φ 1 from (11), however, ignores the background and capture precisely the remaining obstacle. Using u 0 and φ 1 , we expand u as in [START_REF] Grote | Inexact interior-point method for PDE-constrained nonlinear optimization[END_REF] with K = 1 and compute its best L 2 -approximation. Shown in Fig. 5, it is hardly distinguishable from the true u with well-defined sharp contours and a relative L 2 -error below 2%.

In contrast, if we repeat the same experiment with µ ≡ 1, we observe in Fig. 5 how u 0 indeed matches the boundary values of u but fails to capture any additional features inside Ω. Similarly, the first eigenfunction of the Laplacian is independent of u(x) and thus, as expected, carries no information about it. Again, we expand u as in [START_REF] Grote | Inexact interior-point method for PDE-constrained nonlinear optimization[END_REF] but now use the first 1000 eigenfunctions of the Laplacian instead. Although the L 2 best approximation with 1,000 Laplacian eigenfunctions now yields a reasonable approximation of u with 6% relative L 2 -error, the contours are blurred while small high-frequency ripples appear due to the well-known Gibbs-phenomenon. These results illustrate the remarkable accuracy even of but a few eigenfunctions of the AE basis.

Remark 3. For a given profile u, the corresponding AE basis {u 0 , φ 1 , . . . , φ K } of relatively small dimension usually yields a remarkably accurate representation of u. At higher eigenvalues λ m , however, the (mutually orthogonal) AE basis functions φ m become increasingly oscillatory and no longer carry useful information about u -see Fig. 6. Truncating the expansion in (12) at a finite value, K, thus effectively builds regularization into the AEI approach. 9), the first AE eigenfunction φ 1 from (11), the L 2 best approximation of u using {u 0 , φ 1 }. Bottom, from left to right: the harmonic u 0 , the first Laplace eigenfunction φ 1 , the L 2 best approximation of u using the first 1,000 Laplace eigenfunctions {u 0 , φ 1 , . . . , φ 1000 }.

Numerical Results

We shall now illustrate the usefulness and versatility of the AEI method through a series of numerical experiments. Clearly, the squared velocity u of the medium is now unknown and we shall attempt to recover it from boundary measurements by solving [START_REF] Beilina | Approximate Global Convergence and Adaptivity for Coefficient Inverse Problems[END_REF].

First, we compare the adaptive eigenspace representation to a standard grid-based nodal representation of the control u to demonstrate the resulting significant reduction in degrees of freedom. Next we add yet another level of adaptivity by solving the auxil- iary elliptic eigenvalue problem with an adaptive finite element method. We also present numerical experiments that underpin our choices for the background profile u 0 and the varying dimension of V K with increasing frequency. Then, we show that the AE basis functions φ m are highly localized, and hence easily sparsified; thus, the memory requirement of the adaptive eigenspace basis is kept to a minimum. Finally, we demonstrate the robustness of the AEI approach in the presence of noise or missing data and combine it with the sample average approximation (SAA) approach [START_REF] Haber | An effective method for parameter estimation with PDE constraints with multiple right-hand sides[END_REF].

In the entire section, we consider the profile u displayed in Fig. 7, which mimics a layered material with regions of different wave speed. Unless specified otherwise, the typical parameter settings in the numerical experiments are the following: nine equispaced Gaussian sources are located along the upper boundary at (0.1, 0.8), . . . , (0.9, 0.8), whereas the receivers are located on the four lateral boundaries of Ω = (0, 1)×(0, 1). We use secondorder staggered finite differences on a 500×500 Cartesian mesh for the discretization of ( 2), ( 9)- [START_REF] Eisenstat | Choosing the forcing terms in an inexact Newton method[END_REF]. To avoid any inverse crime, the reference solution is computed on a separate finer mesh, which does not contain the coarser computational mesh.

In the AEI algorithm described in Section 2.3, we always set the initial guess to u(x) ≡

1. Starting at the lowest frequency ω = 8, we progressively increase the frequency ω = 10, 12, 14, . . . , 90. The number of eigenfunctions K starts at K = 16 or 32 and increases linearly with the frequency ω. For the optimization, we use a standard truncated Gauss-Newton method [START_REF] Haber | On optimization techniques for solving nonlinear inverse problems[END_REF][START_REF] Nocedal | Numerical Optimization[END_REF] without extra regularization term. The search direction is by a truncated Conjugate-Gradient iteration with the Eisenstat-Walker stopping criterion [START_REF] Eisenstat | Choosing the forcing terms in an inexact Newton method[END_REF][START_REF] Métivier | Full waveform inversion and the truncated Newton method[END_REF][START_REF] Nash | A survey of truncated-Newton methods[END_REF][START_REF] Dembo | Truncated-Newton algorithms for large-scale unconstrained optimization[END_REF] and the step-size by a standard Armijo rule. In the definition of µ in [START_REF] Dembo | Truncated-Newton algorithms for large-scale unconstrained optimization[END_REF], we always set ε = 10 -6 .

In the following numerical experiments, we either use finite differences for the discretizethen-optimize approach, or finite elements for the optimize-then-discretize approach. We shall not emphasize any particular choice for the discretization as it did not affect the results.

Adaptive eigenspace vs. nodal basis

The use of an adaptive eigenspace (AE) basis for the control variable u instead of a standard grid-based nodal basis is the distinguishing feature of the AEI method. Thus, we now compare the AE vs. a standard nodal representation for the reconstruction of u. In both cases, we omit extra Tikhonov-type regularization.

First, we include an additional tenth source located at (0.15, 0.15). The two corresponding reconstructed profiles are shown in Fig. 8. Although both methods recover the essential features of the medium, the AEI method clearly yields much crisper boundaries but also higher accuracy inside the various subregions. Moreover, the AEI method achieves the higher accuracy with fewer than N u = 360 degrees of freedom vs. N u = 501, 000 for the standard nodal representation -see Table 1. Clearly, adding regularization would certainly remove some of the artifacts in the grid-based approach and thus yield a smoother, but not necessarily more accurate, reconstruction. Next, we repeat the previous experiment but now omit all the receivers at the lower boundary of Ω together with the tenth source located at (0. information about the lower part of the medium is available in the data. Nevertheless, as shown in Fig. 9, the AEI method is still able to recover u everywhere inside Ω, unlike the standard grid-based approach. Indeed, as shown in Table 1, the relative L 2 -error for the nodal approach has now almost doubled whereas the error for the AEI method has hardly changed. Although the AEI method uses much fewer control variables than the grid-based approach, i.e. N u = K, the reconstructions appear remarkably accurate and tolerant to missing data.

Adaptive finite element discretization

The AEI method uses two separate computational meshes, the first for the forward and adjoint problems (1), ( 5) and the second for the auxiliary elliptic eigenvalue problems ( 9), [START_REF] Eisenstat | Choosing the forcing terms in an inexact Newton method[END_REF]. So far both meshes were spatially uniform. Here we include yet another level of adaptivity by adapting the finite element (FE) mesh in the solution of ( 9) and ( 11) to better capture small-scale features of u.

In the AEI Algorithm, mesh adaptation is performed after the while loop in step 3b and before the update in step 3c. We use a standard a priori FE adaptive mesh strategy based on the Hessian of the current u [START_REF] Quarteroni | Numerical Models for Differential Problems[END_REF], which is available in the open source software FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF]. Again, the mesh to solve the forward problem (1) and the adjoint problem (5) remains fixed and uniform (200 × 200 Cartesian mesh) throughout the entire computation. For the optimization, we now use a truncated Quasi-Newton (BFGS) method with Wolfe-Powell step-size control.

In Fig. 10, we present the numerical results for the AEI method with and without mesh adaptation. The reconstruction on a uniform triangular mesh yields 4.86% relative L 2error, whereas the reconstruction with adaptive FE strategy yields 4.17% relative L 2 -error, yet with a ten times smaller number of vertices. Hence, we have not only reduced the error in the reconstruction with even crisper and smoother edges, but also greatly reduced the number of degrees of freedom in the control thereby dividing the overall execution time by 2.5. As shown in Fig. 10, the adapted mesh automatically concentrates the degrees of freedom along variations of u and is refined only where it is needed.

Adaptive vs. harmonic background u 0

As shown in Section 3.2, the background state u 0 defined through [START_REF] De Buhan | Logarithmic stability in determination of a 3D viscoelastic coefficient and a numerical example[END_REF] allows the AEI method to accomodate varying boundary data. Here we compare the accuracy in the background state u 0 harmonic adaptive without mesh adaptation 10.73% 4.86% with mesh adaptation 10.90% 4.17% Table 2: Adaptive vs. harmonic background. Relative L 2 -error in u.

Constant vs. adaptive dimension of V K

Our AEI method does not require extra regularization term, such as standard Tikhonov regularization. Instead it restricts the optimization to the AE subspace V K = span{u 0 , φ 1 , . . . , φ K }, which effectively acts as inherent regularization of the inverse problem. Here we focus on the choice of K to understand the regularization effect of the AEI method. To prevent any intrinsic regularization from a very fine mesh or over-abundant data, we omit the receivers at the lower boundary of Ω (missing data) and use a coarser 200 × 200 Cartesian mesh.

Figure 12: Constant vs. adaptive dimension of V K . Top: for constant K = 100, truncated full-Newton method with relative L 2 -error = 14.79% (left) and truncated Gauss-Newton method with relative L 2 -error = 7.95% (right). Bottom: for linearly varying K, truncated full-Newton method with relative L 2 -error = 6.50% (left) and truncated Gauss-Newton method with relative L 2 -error = 6.72% (right).

First, we keep the number of eigenfunctions constant through the entire frequency stepping process. In Fig. 12, we present the numerical results for constant K = 100 using either a truncated full-Newton or Gauss-Newton method. For the full-Newton method, the regularization is not sufficient to reduce the perturbation at the top of the computational domain, although the reconstruction is more accurate than with the nodal basis. In contrast, for the Gauss-Newton method, the regularization is too strong and the method has difficulty reconstructing the kite. Similar conclusions hold for other constant values of K, which underlines the need for K to vary with the frequency ω.

Next, we let K vary linearly with ω, starting at a small number of eigenfunctions (K = 32) to reduce the risk of landing in a false local minimum. As ω increases during frequency continuation, we slowly increase K to capture smaller details of the scatterer. In Fig. 12, both reconstructions are now quite accurate with an L 2 -error of 6.50% for the Newton and 6.72% for the Gauss-Newton method, respectively. Remark 4. The number of eigenfunctions K controls the regularization, similarly to the parameter α in a standard Tikhonov regularization term αR[u]. For Tikhonov regularization, α is initially large but then gradually decreases to zero as the nonlinear iteration approaches the desired minimum. In contrast for our AEI approach, K is initially small but then gradually increases during optimization so that the AE basis includes a higher number of eigenfunctions for better accuracy.

Sparse AEI method

The AEI method greatly reduces the number of parameter values by restricting the optimization to the subspace V K = span{u 0 , φ 1 , . . . , φ K } of much smaller dimension. Since the discretized version of the eigenvalue problem [START_REF] Eisenstat | Choosing the forcing terms in an inexact Newton method[END_REF] leads to a sparse, symmetric and positive definite matrix, the first K eigenfunctions can be efficiently computed via a standard Lanczos iteration [START_REF] Lehoucq | Deflation techniques for an implicitly re-started arnoldi iteration[END_REF]. Still, the storage of the first K eigenfunctions, which are global functions in Ω, may at first appear quite large, especially in three space dimensions.

In fact, much of the information contained in those eigenfunctions is highly localized in space and essentially negligible in most of Ω. Again we consider the example described at the beginning of Section 4. Next, for each eigenfunctions φ m , we monitor all entries smaller than η • φ m ∞ in magnitude. In Fig. 13, we display the percentage of small entries in φ m for each frequency averaged over all m ≤ K; recall that K increases linearly with ω.

As ω increases, the percentage of small entries also increases and quickly saturates above 90% at higher frequencies; hence, more than 90% of all entries are in fact negligible. To save memory space, we can therefore set to zero all those small values and simply replace the eigenfunctions φ m by their sparse approximations.

In Fig. 14, we show the resulting reconstruction with the sparse AEI approach, where all small entries in the eigenfunctions below η φ m ∞ are set to zero for η = 0.1, 0.05, or 0.01. Remarkably, the reconstruction hardly changes and, in fact, is even slightly better (L 2errors below 4%), while saving more than 90% of memory.

Noisy data

To illustrate the robustness of the AEI method with respect to noise, we now add multiplicative noise to the observations for each observation point x i and frequency ω j :

y obs,δ (x i , ω j ) = y obs (x i , ω j )(1 + δ • ξ i,j ),
where ξ i,j are i.i.d. Gaussian random variables with mean zero and variance equal to one. The level of noise is denoted by δ.

Again, we consider the parameter settings described at the beginning of Section 4, but use the adaptive finite element strategy from Section 4.2. In Fig. 15 we present the reconstruction results for two different levels of noise: for δ = 10% (left) we obtain a relative L 2 -error of 4.01% and for δ = 20% (right) we obtain a relative L 2 -error of 5.22%. As shown in Fig. 15, the AEI method is still able to reconstruct the profile without any added regularization and without artifacts due to noise. In fact, at the smaller noise level δ = 10%, the relative L 2 -error is even slightly better than that without noise, 4.01% vs. 4.17%; each case, however, leads to (slightly) different meshes and eigenfunctions due to the adaptive finite element strategy.

Sample Average Approximation

As the number of sources N s in (2) or (6) increases, the cost of computing N s forward and adjoint solutions in (1), ( 5) may become prohibitive. To limit the computational cost without ignoring any of the available data, we consider the sample average approximation (SAA) approach from [START_REF] Haber | An effective method for parameter estimation with PDE constraints with multiple right-hand sides[END_REF], which replaces the sources f by N r "super-shots", N r N s ,

F j = Ns =1 ξ j f , j = 1, . . . , N r , (17)
where the ξ j are i.i.d. random variables with zero mean and unit variance and corresponding observations

Y obs j = Ns =1 ξ j y obs , j = 1, . . . , N r . (18)
During frequency stepping, we choose for each ω a different realization of ξ j -here ξ j = ±1 with probability 0.5. Again, we consider the parameter settings described at the beginning of Section 4, but now with N s = 201 Gaussian sources located at (0.1, 0.8), (0.11, 0.8), . . . , (0.89, 0.9), (0.9, 0.8).

The SAA approach [START_REF] Haber | An effective method for parameter estimation with PDE constraints with multiple right-hand sides[END_REF] with only a single "super-shot", N r = 1, yields the reconstruction shown in Fig. 16. For comparison, we also display the reconstruction without SAA with a single source located at (0.5, 0.8). Although the computational cost of both approaches is identical, the SAA approach yields better accuracy, remarkably so, given that only a single ("super-shot") source is used for the approximation.

Figure 16: SAA approach. Reconstruction of u with a single source. Left: without SAA and with N s = 1 (10.05% relative error). Right: with a single SAA "super-shot", i.e. N r = 1, from N s = 201 sources (5.79% relative error). Note that the computational effort is identical.

Concluding remarks

We have presented a nonlinear optimization method for the solution of inverse scattering problems in the frequency domain, when the scattered field is governed by the Helmholtz equation. Instead of a standard (FD or FE) grid-based representation, the unknown (squared) sound speed u is projected to the finite-dimensional subspace V K = span{u 0 , φ 1 , . . . , φ K } of much smaller dimension. The "background" u 0 is determined by solving (9) whereas the remaining orthonormal basis functions φ m are determined by computing the first K eigenfunctions in [START_REF] Eisenstat | Choosing the forcing terms in an inexact Newton method[END_REF]. The time-harmonic inverse medium problem is formulated as a PDE-constrained optimization problem and solved by an inexact truncated Newton or quasi-Newton iteration. During the optimization process, which includes frequency continuation, both the basis and the dimension of V K are repeatedly adapted to the current iterate. The full Adaptive Eigenspace Inversion (AEI) Algorithm is given in Section 2.3. For monotonic one-dimensional or layered media, we have proved that our choice for u 0 is in fact optimal. For arbitrary media, our numerical results suggest that it is clearly superior to a straightforward harmonic extension from the known boundary values. Together with but a few eigenfunctions, the adaptive eigenspace basis yields a remarkably accurate representation of u. At higher eigenvalues, the eigenfunctions become increasingly oscillatory while no longer carrying useful information about u. Hence adapting the dimension of the eigenspace basis in [START_REF] Grote | Inexact interior-point method for PDE-constrained nonlinear optimization[END_REF], effectively builds regularization into the inversion, so that no additional Tikhonov regularization is needed -see Remark 4. As the frequency ω increases, smaller scale features of the scatterer become visible and the dimension of V K thus ought to increase accordingly.

In contrast to a standard grid-based nodal representation, the AEI method uses much fewer control variables for u. Still the reconstructions are remarkably accurate, display less artifacts and prove more tolerant to partial or missing data. Moreover, our AEI method leads to a significant reduction in execution time and proves robust with respect to added noise. By combining it with the Sample Average Approximation (SAA) approach from [START_REF] Haber | An effective method for parameter estimation with PDE constraints with multiple right-hand sides[END_REF], it also efficiently handles large numbers of sources.

Since the discrete version of the eigenvalue problem [START_REF] Eisenstat | Choosing the forcing terms in an inexact Newton method[END_REF] leads to a sparse, symmetric and positive definite matrix, the first K eigenfunctions can be efficiently computed via a standard Lanczos iteration. If finite element mesh adaptation is used for the numerical solution of [START_REF] Eisenstat | Choosing the forcing terms in an inexact Newton method[END_REF], small-scale features and interfaces are captured with even greater accuracy in the reconstruction, without increasing the computational effort. Although the eigenfunctions are global, their information content is highly localized in space so that most entries are in fact negligible.

Our current work involves the extension of the AEI method to multi-parameter inverse scattering problems [START_REF] Prieux | Multiparameter full waveform inversion of multicomponent Ocean-Bottom-Cable data from Valhall field. Part 1: imaging compressional wave speed, density and attenuation[END_REF][START_REF] Prieux | Multiparameter full waveform inversion of multicomponent Ocean-Bottom-Cable data from Valhall field. Part 2: imaging compressional and shear-wave velocities[END_REF][START_REF] Tarantola | Inversion of seismic reflection data in the acoustic approximation[END_REF]. We shall report on those results elsewhere in the near future.

Figure 1 :

 1 Figure 1: One-dimensional case. Comparison of adaptive and harmonic choices for u 0 . Left: u monotonically increasing; right: u arbitrary.

25 λ 25 =

 2525 1793.3732

Figure 2 :

 2 Figure 2: One-dimensional case. The true profile u (top left), together with the eigenfunctions φ 1 , φ 2 , φ 3 , φ 24 and φ 25 from (15) with ε = 10 -4 .

1 AEI: µ from ( 14 )Figure 3 :

 1143 Figure 3: One-dimensional case. The true smooth profile u (top) and the first eigenfunction φ 1 from (15): µ ≡ 1 (left), and µ as in (14) (right).

Figure 4 :

 4 Figure 4: Two-dimensional case. Left: the exact profile u. Right: cross-cut at x 2 = 0.4.

Figure 5 :

 5 Figure 5: Two-dimensional case. Top, from left to right: the background u 0 from (9), the first AE eigenfunction φ 1 from (11), the L 2 best approximation of u using {u 0 , φ 1 }. Bottom, from left to right: the harmonic u 0 , the first Laplace eigenfunction φ 1 , the L 2 best approximation of u using the first 1,000 Laplace eigenfunctions {u 0 , φ 1 , . . . , φ 1000 }.

Figure 6 :

 6 Figure 6: First AEI eigenfunction φ 1 (left) and last AEI eigenfunction φ N (right) on a 201 × 201 grid.

Figure 7 :

 7 Figure 7: Left: true profile u; Right: real part of y with f at (0.1, 0.8) and ω = 90

Figure 8 :

 8 Figure 8: Adaptive eigenspace vs. nodal basis. Full boundary data: reconstruction with a nodal basis (left) or an AE basis (right).

Figure 9 :

 9 Figure 9: Adaptive eigenspace vs. nodal basis. Missing boundary data: reconstruction with a nodal basis (left) or an AE basis (right).

Figure 10 :

 10 Figure 10: Adaptive FE discretization. Top: uniform triangular mesh with 30,534 vertices (left) and recovered u (right). Bottom: with mesh adaptation, final mesh for ω = 90 with 2,783 vertices (left) and recovered u (right).

Figure 13 :

 13 Figure 13: Sparse AEI method. Percentage of entries smaller than η • φ m ∞ in magnitude averaged over all m ≤ K at each frequency.

Figure 14 :

 14 Figure 14: Sparse AEI method. Reconstruction of u with sparsified eigenfunctions. Left: with η = 0.1 (error = 3.83%). Center: η = 0.05 (error = 3.81%). Right: η = 0.01 (error = 3.90%)

Figure 15 :

 15 Figure 15: Noisy data for varying noise level δ. Left: δ = 10% (L 2 -error = 4.01%). Right: δ = 20% (L 2 -error = 5.22%).

Table 1 :

 1 15, 0.15). Hence, much less Adaptive eigenspace vs. nodal basis. Relative L 2 -error and number of degrees of freedom for u.

		Full data	Partial data
	Basis L 2 -Error	N u	L 2 -Error	N u
	Nodal 15.91% 501,000 30.24% 501,000
	AEI	4.65%	≤ 360	4.80%	≤ 360

reconstruction for two different choices for u 0 :

-harmonic: u 0 is computed only once at the beginning of the AEI Algorithm by solving ( 9) with µ ≡ 1 -adaptive: u 0 is recomputed at each frequency step by solving [START_REF] De Buhan | Logarithmic stability in determination of a 3D viscoelastic coefficient and a numerical example[END_REF] with µ as in [START_REF] Dembo | Truncated-Newton algorithms for large-scale unconstrained optimization[END_REF] Figure 11: Adaptive vs. harmonic background. Top: without mesh adaptation, harmonic u 0 (left) and adaptive u 0 (right). Bottom: with mesh adaptation, harmonic u 0 (left) and adaptive u 0 (right).

In Fig. 11, we compare the recovered u for the above two different background states u 0 , both either with or without the adaptive mesh strategy from Section 4.2. Clearly, the adaptive background u 0 improves the accuracy of the reconstruction as the relative L 2errors are approximately halved -see Table 2. Moreover, the artifacts near the exterior boundary, visible in the harmonic background approach, are now completely absent.

Hence the more accurate background state u 0 enables the AEI method to better approximate the remainder u -u 0 with the AE basis. Still, the adaptive strategy is also slightly more expensive since it adaptively recomputes u 0 for every frequency. These conclusions also appear to hold if an adaptive FE strategy is included in the solution of ( 9), [START_REF] Eisenstat | Choosing the forcing terms in an inexact Newton method[END_REF], which in fact seems to have little impact on the overall accuracy here.