A DISCRETE LIOUVILLE IDENTITY FOR NUMERICAL RECONSTRUCTION OF SCHR¨ODINGER POTENTIALS

,


Introduction
Consider the boundary value problem (1) L σ,q v ≡ −∇ • [σ∇v] + qv = 0 in Ω, v = f on ∂Ω, in a simply connected, bounded domain Ω ⊂ R 2 with C 2 boundary ∂Ω, and f ∈ H 1/2 (∂Ω).The positive and bounded coefficient σ(x) is called the conductivity and q(x) is the Schrödinger potential.They are the unknowns in inversion, to be determined from the Dirichlet to Neumann (DtN) map Λ σ,q : H 1/2 (∂Ω) → H −1/2 (∂Ω) defined by where n is the unit outer normal at ∂Ω and v solves (1).The case q = 0 is known as the inverse conductivity or electrical impedance tomography problem.The case σ = 1 is the inverse Schrödinger problem.Our goal in this paper is to introduce a novel method, based on parametric model reduction, for the numerical reconstruction of the solution of the inverse Schrödinger problem (σ = 1) in the absorptive case q ≥ 0. Parametric model reduction is mainly used for approximating efficiently the response of dynamical systems for design, optimal control and uncertainty quantification [3].We are interested in parametric model reduction for improving the inversion process.This is a largely unexplored direction, but recent progress has been made in [9,21] for parabolic equations, in [28] for the wave equation, in [5,8,7,6] for the inverse conductivity problem and in [4] for a related inverse spectral problem.The construction of the reduced models varies between problems because it must respect the underlying physics.For example, the projection-based reduced models in [28] approximate the wave propagator and use causality, whereas the projection-based models in [9,21] for parabolic equations are obtained by rational approximation of the transfer function, the Laplace transform with respect to time of the measurement map.The parametric reduced models for the inverse conductivity problem are not projectionbased.They are resistor networks constructed from very accurate approximations of the DtN map and the resistances in the network play the role of the parameters in the parametric model reduction.These resistor network reduced models rely on the graph theory developed in [12,14,15,16,17,18,19,20,27].
The advantage of using parametric reduced models for inversion is that they can lead to iterative algorithms that perform much better than the usual least squares data fit methods.They converge in one or two iterations and give quantitatively superior reconstructions.The disadvantage is that it is difficult to construct good parametric reduced models.These must retain the structure of the governing partial differential equation so we can extract the unknown parameters from them, and capture correctly important phenomena such as the decreased sensitivity of the measurements to changes in the parameters away from the boundary in inverse elliptic and parabolic problems.
In this paper we show how to use the resistor networks with circular planar graphs, the reduced models for the inverse conductivity problem in [5,7,6,8], to solve the inverse Schrödinger problem with absorptive potentials in two dimensions.The inverse conductivity and Schrödinger problems are connected in the continuum by a well known Liouville identity [29].Here we show how to connect them in the discrete (network) setting.This is difficult because conductivities are defined on edges of the graph of the network, and the potential is associated with the nodes.The parametric reduced models (networks) obtained as in [8], which can be interpreted as five point stencil difference schemes for Schrödinger's equation, encode information about the unknown potential q, but this is not restricted to the diagonal of the finite difference operator, as expected.Thus, it is not straightforward to obtain a reconstruction of q from the network.
Discrete Liouville identities have been proposed before in [26] and in [1,2] for the inverse Schrödinger problem on networks.However, these studies are not concerned with the connection with the continuum problem, and in fact it is not clear if the discrete Schrödinger problem considered in [1,2] is consistent with measurements of the DtN map in the continuum.
Here we derive a discrete generalized Liouville identity for solving the continuum inverse Schrödinger problem with networks.The networks, with graph G, are parametric reduced models that approximate the DtN map Λ 1,q .The Liouville identity is defined on the line graph G of G, which establishes an isomorphism between the edges of G where the conductivities lie and the nodes of G that support the Schrödinger potential.We use the Liouville identity to formulate a preconditioned Gauss-Newton inversion algorithm.The preconditioner is obtained from the parametric reduced model and leads to an efficient method, as demonstrated with numerical simulations.
1.1.Contents.We start in §2 by recalling how the Liouville transform relates the inverse problem for the absorptive Schrödinger equation to the inverse conductivity problem.The inverse problem for resistor networks stated in §3 is the discrete analogue of the inverse conductivity problem.We can transform it to a discrete analogue of the inverse Schrödinger problem using the generalized Liouville identity derived in §3.4.The connection between the continuum and discrete inverse problems is in §4.We show in section §4.1 how to relate measurements of the DtN map of the Schrödinger equation to the discrete DtN map of a unique resistor network.In §4.2 we use the generalized Liouville identity to obtain a discrete Schrödinger potential from this resistor network.The discrete potential is used in §5 to obtain a reconstruction of the continuum Schrödinger potential.The performance of the inversion method is assessed with numerical simulations in §6.We conclude with a summary in §7.

Continuum inverse conductivity and Schrödinger problems
A well-known relation between the Schrödinger L 1,q and conductivity L σ,0 differential operators is through the Liouville identity (see e.g.[29]) where σ > 0, σ ∈ C 2 (Ω) and ( 4) In the left hand side of (3) we have the composition of three linear operators.Two of them are the operator v → σ −1/2 v denoted, in an abuse of notation, by σ −1/2 .In linear algebra two matrices A and B are congruent if there is an invertible matrix S such that A = SBS T (see e.g.[25]).Borrowing the terminology, we say that two linear differential operators A and B are Liouville congruent if there is a positive function s for which A = s • B • s.This allows us to restate the Liouville identity as follows: The operators L σ,0 and L 1,q are Liouville congruent when q is given by (4).
Note that the Dirichlet to Neumann maps of the conductivity and Schrödinger problems satisfy [29] (5) Thus, when σ equals a constant near ∂Ω, the DtN maps are also Liouville congruent.When σ has a nonzero normal derivative at ∂Ω, the DtN maps are congruent up to a diagonal (or multiplication by a function) operator.We assume henceforth that σ| ∂Ω = 1.Equations (3)- (5) show that in the continuum setting any inverse conductivity problem for sufficiently smooth σ can be formulated as a Schrödinger problem.The converse is true only for potentials q that give a positive solution σ of (4).This is the case for absorptive potentials q ≥ 0. Indeed, the strong maximum principle (see e.g.[22,§6.4,Theorem 4]) guarantees that when q ≥ 0, the solution of (6) −∆s + q s = 0, in Ω, s = 1, on ∂Ω, satisfies 0 < s ≤ 1.Hence, σ = s 2 satisfies (4) and the operators L σ,0 and L 1,q are Liouville congruent.We can take the notion of Liouville congruence further than the classical relations (3)-( 4).Instead of relating operators L σ,0 and L 1,q , we can also consider operators L σ1,0 and L σ0,q , for two positive conductivities σ 0 and σ 1 and a potential q.It follows by straightforward calculations that (7) ( with q given by ( 8) The importance of this generalized Liouville identity becomes clear in the next section, where we derive a discrete analogue of ( 7)-(8).

Discrete inverse conductivity and Schrödinger problems
As stated in the introduction, it is useful to view our study in the context of parametric model reduction for inversion.We need reduced models that keep the underlying structure of the differential operators L σ,q so that we can obtain reconstructions of σ and q from them and, in addition, give better conditioned optimization algorithms than the usual output least squares.Such parametric reduced models are known for the inverse conductivity problem in two dimensions.They are based on the circular planar resistor networks studied in [20,27] and are described briefly below.We refer to [5,7,6] and the review [8] for details on how to use the parametric reduced models to determine a discrete Laplacian which is consistent with the measurements of the DtN map in the continuum setting, and to recover the unknown conductivity.
In this section we describe the basic tools for extending the inversion approach to the Schrödinger problem.We begin in §3.1 with the formulation of the discrete analogues of the inverse conductivity and Schrödinger problems.Then we review briefly in §3.2 the relevant facts about the resistor networks.The line graph introduced in §3.3 and the discrete Liouville identity defined in §3.4 allow us to use the networks for solving the inverse Schrödinger problem.
3.1.The discrete conductivity and Schrödinger operators.The discrete structure of a resistor network is an undirected graph G = (V, E), where V is a finite set of vertices (nodes) and the edge set E is a subset of {{i, j} | i, j ∈ V, i = j}.We denote the set of functions from a finite set X to R by R X , and write f (x) for f ∈ R X and x ∈ X.All functions and operators related to finite sets are written henceforth in bold to distinguish them from their continuum counterparts.
The discrete gradient on a network is the linear operator A sign needs to be specified for each edge, but as long as this sign convention is fixed, it does not change the subsequent definitions.
A resistor network is defined by its graph G = (V, E) and a positive discrete conductivity function γ ∈ (0, ∞) E .We may also define a discrete Schrödinger potential q ∈ R V on the vertices of the graph, and introduce the discrete Schrödinger operator L γ,q : R V → R V , which maps potentials Here D * : R E → R V is the adjoint of the discrete gradient D and the product f g for functions f , g ∈ R X (where The special case L γ,0 is the weighted graph Laplacian with weight γ (see e.g.[11,18]).This is a resistor network.The edge e of the graph is a resistor with conductance γ(e) and the nodes represent electrical connections.The absorptive case L γ,q with q > 0 is also a resistor network, with each node i connected to the ground (zero potential) by a resistor with conductance q(i).
To define the DtN map of a resistor network, we collect the nodes V in two disjoint sets: the boundary nodes in the set B ⊂ V , and the interior nodes I = V \B.The Dirichlet boundary value problem for the network is to find the potential u ∈ R V such that (12) (L γ,q u) I = (L γ,q ) II u I + (L γ,q ) IB u B = 0, and Here u I is the restriction of u to I, and the linear operator (L γ,q ) BI : R I → R B is defined by (L γ,q ) BI v = (L γ,q u) B where u I = v and u B = 0.The other operators are defined similarly.
Lemma 3.1.When the subgraph of G induced by I is connected, the Dirichlet problem for q ≥ 0 has a unique solution for any f ∈ R B .
Proof.The case q = 0 is proven in [19,13] using a discrete analogue to the maximum principle.Here we proceed by first showing that (L γ,0 ) II is positive definite.
If the subgraph of G induced by I is connected, we have that the weighted graph Laplacian on this subgraph, which is denoted by L γ I ,0 , is a symmetric positive semidefinite matrix with a one dimensional nullspace spanned by the constant vector 1 ∈ R I .The sub-matrix (L γ,0 ) II can be written where E is a diagonal matrix, whose only non-zero entries are positive and correspond to the interior nodes that have an edge in common with a boundary node.
and therefore (L γ,0 ) II must be positive definite.When q ≥ 0, notice that (L γ,q ) II = (L γ,0 ) II + diag(q I ) is the sum of a positive definite and a positive semidefinite matrix, and must then be positive definite.This proves the solvability result for q ≥ 0.
The Dirichlet to Neumann (DtN) map Λ γ,q : R B → R B associates to the boundary potentials f ∈ R B the boundary currents (13) (L γ,q u) B = (L γ,q where u solves (12).The linear map Λ γ,q is symmetric, and in the case q = 0, it has a one dimensional null space spanned by the vector of all ones 1 ∈ R B .To write it more explicitly, let us simplify notation as L BB := (L γ,q ) BB , and similar for the BI, IB and II blocks of L γ,q .By solving for u I in (12) and substituting in (13), we can write Λ γ,q as a Schur complement of L γ,q , (14)

3.2.
The inverse problem for resistor networks.The discrete analogue of the inverse conductivity problem is: Find the conductivity γ ∈ (0, ∞) E from the DtN map Λ γ,0 , given the underlying graph G = (V, E).In what follows we refer to the solution as γ(Λ γ,0 ).The discrete inverse problem has been solved for circular planar graphs [20,16,12] i.e., when G can be embedded in the plane with no edge crossings so that the boundary nodes are on a circle and the interior nodes inside the circle.The recoverability result in [20,16,12] can be summarized as follows.
Theorem 3.1.The DtN map for a critical circular planar resistor network with positive conductivity determines uniquely the conductivity.
A circular planar graph is said to be critical if the two following conditions hold: (i) any two sets of boundary nodes with the same number p of elements and lying on disjoint segments of the boundary circle can be linked with p disjoint paths.(ii) the deletion of any edge in the graph breaks condition (i).
For example, among the family of graphs C(m, n) illustrated in Figure 1, with n odd, the graphs are critical when m = (n − 1)/2 i.e., when the number |E| of edges is equal to the number n(n − 1)/2 of entries in the strictly upper triangular part 1 of Λ γ,0 .
1 Since Λ γ,0 is symmetric with rows summing to zero, its independent entries lie in its strictly upper triangular part.
Remark 3.1.There are other families of critical circular planar graphs that allow the number of boundary nodes to be even, such as the pyramidal networks in [7,20].In this paper we use graphs C n(n − 1)/2, n with n odd.
There exist direct (optimization-free) reconstruction methods for networks with graphs C((n − 1)/2, n).They solve the discrete inverse conductivity problem in a finite number of algebraic operations by peeling off sequentially the layers of the network [16].Layer peeling and any other methods of reconstructing γ(Λ γ,0 ) become unstable as the number of edges in the network grows.This is similar to the instability of the continuum inverse conductivity problem.Thus, in practice it is impossible to approximate the solution of the continuum problem by the discrete conductivity of larger and larger networks.However, we can relate the continuum and discrete problems by using networks of modest size, as explained in §4.
In §2 we discussed a well-known connection between the inverse conductivity and Schrödinger problems in the continuum via the Liouville transform and its generalization ( 7)- (8).To use the existing results on the discrete inverse conductivity problem, we seek to relate the discrete conductivity γ to a discrete Schrödinger potential q.This is not straightforward because the discrete conductivity γ ∈ R E is defined on the edges of the graph, while q ∈ R V is defined at the nodes.To resolve this difficulty, and to derive the discrete analogue of the generalized Liouville transform, we introduce below the line graph and a discrete Schrödinger operator associated with it.

The line graph.
To a resistor network with graph G = (V, E) we associate a line graph G = ( V , E).All quantities and operators associated with the line graph appear henceforth with tilde.The vertices (nodes) V and edges E of the line graph are defined in terms of V and E of the original graph so that: (i) there is one vertex of G per edge of G, and so V is isomorphic to E. (ii) there is an edge between two vertices of G if and only if the corresponding edges of G share a node.
An example of a graph and its line graph is given in Figure 2.
If γ ∈ (0, ∞) E is a conductivity on the graph G, we define the conductivity γ ∈ (0, ∞) E on the associated line graph by geometric averages of γ.To be more precise, if e, e ∈ E are distinct edges in G that share a vertex, then {e, e } ∈ E is an edge of the line graph G and (15) γ({e, e }) = γ(e)γ(e ).
The line graph together with the conductivity γ is itself a resistor network.As in section 3 we can define a discrete gradient D : R V → R E , and for some q ∈ R V , we can define a discrete Schrödinger operator L γ, q : Since V is isomorphic to E and R V is isomorphic to R E , we can think of L γ, q as a discrete Schrödinger operator acting on the edges E of G, with a Schrödinger potential q defined on E as well (i.e., a function in [0, ∞) E ). 3.4.The discrete generalized Liouville identity.With the line graph defined above, we derive here a discrete analogue to the generalized Liouville identity ( 7)- (8).The need for this identity can be explained as follows.Unlike the continuum case, where σ ≡ 1 corresponds to L 1,0 = ∆, the case of constant discrete conductivity γ = 1 has no special significance in the discrete setting.However, the reduced model (resistor network) for the reference conductivity σ o ≡ 1 (i.e., q ≡ 0) is a good approximation of the Laplacian on the graph G, as shown in [5,7,8].The geometric averages of its edge conductivities γ 0 define the line graph Laplacian.We explain in the next section how to obtain the reduced model (resistor network) for the unknown parameter σ, related to the unknown q.The discrete Liouville identity stated in the next theorem allows us to estimate q from this reduced model and the line graph Laplacian.
Theorem 3.2 (Discrete Generalized Liouville Identity).Let γ 0 and γ 1 be two conductivities in (0, ∞) E .Then their associated weighted graph Laplacians on the line graph satisfy where division and power of vectors is understood componentwise and We say that L γ1,0 and L γ0, q are Liouville congruent, when q is given by (18).
Proof.We obtain from (16) and the definition of the discrete gradient D that for any γ ∈ (0, ∞) E and e, f ∈ E, e = f , Hence for any γ ∈ (0, ∞) E , and distinct e, f ∈ E, we have Then the off-diagonal entries of (17) follow from writing (19) for γ equal to γ 0 and γ 1 , and equating.It remains to calculate q, so that the diagonal entries of ( 17) are equal.Because the off-diagonal entries are equal, this is the same as equating the sum of the rows Since D1 = 0, we must also have L γ,0 1 = 0, so Solving for q in this equation gives (18).
Remark 3.2.Theorem 3.2 shows that two weighted graph Laplacians on the line graph, L γ1,0 and L γ0,0 , are congruent up to a diagonal term which involves the potential q.Moreover, the matrix that carries the congruence relation is diagonal, given by diag(γ 1 /γ 0 ) 1/2 .This is the discrete analogue of the congruence relation (7), and the definition of q is the analogue of (8).

Connection between the continuum and discrete inverse problems
To solve the continuum inverse problem, we obtain in §4.1 the networks as parametric reduced models for lumped measurements of the DtN map Λ 1,q .These measurements define the DtN map Λ γ,0 of the network with discrete conductivity γ, related to the unknown potential q via the discrete generalized Liouville identity, as described in §4.2.

4.1.
Resistor networks as parametric reduced models.The exponential instability of the inverse Schrödinger and conductivity problems limits the resolution of the reconstructions.In our context this means that the size of the reduced models (the networks) cannot be too large, specially for noisy data.For the networks with graphs C n(n − 1)/2, n considered here, the number n of boundary nodes can be chosen based on the noise level, as explained in Appendix A. We show here how to define from the measurements of Λ 1,q a discrete DtN map of the network with graph C n(n − 1)/2, n , our reduced model for the unknown potential q and therefore conductivity.
In principle, we could just take point measurements of Λ 1,q at the n boundary nodes.However, studies such as [23] tell us that smooth boundary currents (i.e., with Fourier transform supported at small frequency) are better for sensing inhomogeneities inside the domain.Thus, we smooth the measurements of Λ 1,q by lumping them at the n boundary points of the network.The lumping is achieved with n compactly supported smooth functions φ 1 , . . ., φ n , whose disjoint supports are numbered consecutively in ∂Ω, in counter-clockwise order.We normalize them by ∂Ω φ j = 1.
From Λ 1,q we obtain the n × n data matrix φ i , Λ 1,q φ j for i, j = 1, . . ., n, where •, • denotes the H 1/2 (∂Ω), H −1/2 (∂Ω) duality pair.It follows from the relation (5) between the DtN maps that for absorptive potentials q ≥ 0 there is a unique conductivity σ satisfying (4) and σ| ∂Ω = 1, so that (20) for i, j = 1, . . ., n.Since σ| ∂Ω = 1, we can simplify the first term in (20), and the second term vanishes for i = j because the φ i have disjoint supports.We then get (21) φ i , Λ 1,q φ j = φ i , Λ σ,0 φ j , for i, j = 1, . . ., n, i = j, so the off-diagonal entries of the data matrix are lumped measurements of the DtN map for the associated conductivity problem.We denote this map by M (q), and 0.4 0.7 define its entries M ij as ( 22) This definition satisfies the conservation of currents relation M 1 = 0, for 1 ∈ R B , the vector of all ones at the boundary nodes.Since the measurement matrix M corresponds to taking measurements for a conductivity problem, we can use Theorem 1 in [5] to get the following result.Theorem 4.1.Let q ≥ 0 and M (q) be the n × n measurement matrix of Λ 1,q obtained as in (22).There is a unique reduced model, which is the resistor network with DtN map equal to M (q), graph C n(n − 1)/2, n , and discrete conductivity γ M (q) .

4.2.
From resistor networks to Schrödinger potentials.We have now shown that there is a unique parametric reduced model (resistor network) with DtN map M (q) and graph C(n(n − 1)/2, n).It remains to determine the Schrödinger potential from its discrete conductivity γ M (q) .
A naive approach to reconstructing the potential would be to first obtain a continuum conductivity σ from γ(M (q)) e.g., using the method based on resistor networks from [5], and then use the continuum Liouville identity (4) to get q from σ.This does not work because σ is never exact and errors are amplified when taking the Laplacian in (4).To illustrate this, we show in figure 3 that the conductivities corresponding to two very different Schrödinger potentials are hard to tell apart.
Our inversion algorithm takes a different route.Its outline is given below, and all the steps are described in detail in the next section.
(1) Choose the number n of boundary nodes of the graph C(n(n − 1)/2, n).
(3) Find discrete conductivities γ 1 ≡ γ(M (q)) and γ 0 ≡ γ(M (0)).( 4) Estimate the average q avg of the unknown potential q from M (q).( 5) On the line graph associated to C(n(n − 1)/2, n) compute the geometric averages γ 0 from γ(M (0)) and form the discrete Laplacian L γ0,0 .( 6) Use the discrete generalized Liouville identity identity ( 17)-( 18) with γ 1 , γ 0 from step (3) and L γ0,0 from step ( 5) to obtain the discrete potential q. (7) Use the discrete potential q and the estimated average q avg to reconstruct the unknown Schrödinger potential q.The use of γ 0 at step (3) and of q avg at step ( 7) is a calibration so that the reconstruction is exact for q ≡ 0 and q ≡ q avg .In the context of parametric model reduction, we work with three reduced models: one for the unknown q which we wish to determine, one for the zero potential and one for the constant q avg potential, the average of the unknown q.We ask that the reduced models share the same discrete Laplacian, calculated for q = 0.This discrete Laplacian is a finite difference scheme with steps defined by γ 0 .The average potential at step ( 7) is used as a scaling.We motivate this scaling by noticing that when Ω is the unit disk and q is radially symmetric (i.e.q ≡ q(r)), the problem of finding q from M (q) is related to the inverse spectral problem of finding q from the spectrum of the operator −∆ + q with Dirichlet boundary conditions [8].To solve the inverse spectral Schrödinger problem it is essential to estimate q avg accurately as it corresponds to a shift in the spectrum (see e.g.[10,4]).In our case it is difficult to relate M (q) to the spectrum of −∆ + q, so we use step (7) to ensure that the reconstruction is correct for the constant q avg .If the parametric reduced models share the discrete Laplacian i.e., they give accurate approximations of M (q) for a set of functions q on the grid defined by γ 0 , we expect that they account for the asymptotes of the spectrum.

From discrete Schrödinger potentials to continuum ones
A numerical estimate of the solution of the inverse Schrödinger problem is conventionally obtained with the optimization (output least-squares) formulation (23) min where • F is the Frobenius norm and K(q) denote measurements of Λ 1,q .This formulation requires regularization, based on prior information about q that may not be available.The method is widely used but has a high computational cost, tends to get stuck in local minima and gives solutions that are biased by the priors.Instead of minimizing the misfit in the data, we apply a non-linear mapping Q : R n×n → R n(n−1)/2 that acts as an approximate inverse of M (q), our lumped measurements of Λ 1,q .Thus, we solve instead the optimization problem (24) min with the Gauss-Newton iteration (25) for k = 0, 1, . . .
Here † denotes the Moore-Penrose pseudoinverse [24] and DM [q] : L 2 (Ω) → R n×n and DQ[M ] : R n×n → R n(n−1)/2 are the Jacobians of the mappings M (q) and Q(M ).Following the approach in [5] for the inverse conductivity problem, we define the nonlinear preconditioner Q(M ) by solving a discrete Schrödinger inverse problem with data M ( §5.1).The initial guess q 0 is obtained by a linear interpolation of the entries of Q(M (q true )) on a carefully chosen grid §5.3.This ensures that the Gauss-Newton method converges quickly, usually in one step.
5.1.The nonlinear preconditioner.The mapping Q(M ) is defined by solving a discrete inverse Schrödinger problem with data M , and ensuring that with q avg estimated as in §6.1.The computation of Q involves the following steps: (1) Find the parametric reduced models for the zero potential, the constant potential q avg and the unknown q.That is, determine the discrete conductivities γ 0 ≡ γ(M (0)), γ avg ≡ γ(M (q avg )), and γ(M ).( 2) Use the discrete Liouville identity (18) with reference conductivity γ 0 to find the discrete potentials q(γ(M )) and q(γ avg ) from γ(M ) and γ avg .(3) The map Q is defined by ( 26) where the division is understood componentwise.
The scaling at step (3) ensures that the map Q(M ) gives exact reconstruction of the constant potential q = q avg , under the following assumption: Assumption 5.1.All entries of q(γ avg ) are nonzero.
It is not clear how to prove that this holds, but the assumption has been verified numerically for many graphs and values of q avg on the unit disk.
5.2.The sensitivity functions.A good preconditioning map Q means that Q(M ) is some approximation of the identity.We illustrate this numerically in figure 4 where we show the sensitivity or Fréchet derivative of some of the entries of the vector Q(M (q)) to changes in q.These are the "rows" of the Jacobian of the preconditioned mapping, DQ[M (q)]DM [q].The sensitivity functions are evidently well localized, so to first order approximation the entries Q(M (q)) are local averages of q.However DQ[M (q)] is not a preconditioner in the usual sense because it has a non-trivial nullspace, as stated in the next proposition.This is handled by taking its pseudoinverse in the Gauss-Newton iteration (25).Proof.We first characterize the left and right null space of the linearization of the mapping q(γ) defined by the discrete Liouville identity (18) with discrete reference In particular, if we let δγ = γ we obtain from (18) that D q[γ]γ = 0, so γ is a right null vector of the Jacobian.This may be seen directly from (18), by realizing that the function γ → q(γ) is homogeneous of degree zero2 .A similar calculation gives that γ/γ 0 is a left null-vector for D q[γ], Going back to the definition of Q(M ), we recall that for critical circular planar graphs, the Jacobian Dγ[M ] is invertible [20, §12].Thus, step (1) in the definition of Q gives an invertible linearization map.
Step (3) is also invertible because it involves the multiplication with a diagonal matrix with nonzero entries, so the left and right null-vectors of DQ[M ] can be found from those of D q[γ(M )].
The left null vector follows from z/ q avg = γ(M )/γ 0 being in the left nullspace of D q[γ(M )], as shown above, where the division is understood componentwise.For the right null vector, consider a scalar h > −1.Then by the homogeneity of order 1 between the DtN map and the conductors in the network, we must have This and definition (26)  The "x" are for q = 1 and the "•" for q = 3.

5.3.
The initial guess.The localization of the sensitivity functions displayed in Figure 4 allows us to view Q(M (q true )) as an approximation of q true at the maxima of the sensitivity functions.These maxima define the nodes of the sensitivity grids used in [7] for the inverse conductivity problem.The sensitivity grids depend weakly on the Schrödinger potential, as illustrated numerically in figure 5 for two constant Schrödinger potentials.Thus, we can compute them ahead of time.
To obtain a good initial guess q 0 of the iteration (25), we linearly interpolate the values of Q(M (q true )) on a Delaunay triangulation of the sensitivity grid nodes.Some of these initial guesses are illustrated in figures 6 and 7, with the same color scales as the corresponding q true .We note that q 0 is already close to the actual Schrödinger potential, capturing not only its geometrical features but also its magnitudes.Moreover the computation of q 0 is inexpensive, because all the operations involved in the calculation of Q(M ), with given M , are carried on a relatively small network.

Numerical results
We show here numerical reconstructions obtained from data M (q true ) approximated by solving (1) numerically with a finite volume method on a fine grid, and then using the lumped measurement formulas (22).With the same finite volume method, fine grid and measurement formulas, we also compute M (0) so that we can obtain γ 0 ≡ γ(M (0)) needed for the computation of Q[M ].This is to avoid discretization discrepancies in the calculation of the preconditioner mapping.In practice, one may use any fine grid that gives a good approximation of M (0).6.1.Estimating the average Schrödinger potential.To compute Q[M ] we need to estimate the average q avg of the Schrödinger potential from the data M (q true ).We do this with a direct search.Explicitly, given some trial values q 1 , . . ., q m for the average, we estimate it as q avg = argmin q∈{q1,...,qm} M (q) − M (q true ) F , where • F is the Frobenius norm.To avoid the inverse crime, in this estimation we compute {M (q j )} 1≤j≤m on a different fine grid than the one used for computing the synthetic data M (q true ).6.2.Calculating the non-linear preconditioner mapping.The computation of Q[M ] involves γ avg , the discrete conductivity corresponding to the constant Schrödinger potential q ≡ q avg .We obtain it from M (q avg ), computed as described above for M (q true ) and M (0).The fine grid in the computation of M (q avg ) is the same as that in the calculation of M (0), to avoid discretization discrepancies.6.3.The Gauss-Newton iterations.The reconstructions obtained with Gauss-Newton iteration (25) for a smooth and a piecewise constant potential are shown in figures 6 and 7, respectively.We include only the initial guess q 0 (which is obtained as explained in 5.3) and the first iterate q 1 .Subsequent iterates q k , k ≥ 2 are not shown because they are indistinguishable from q 1 .Thus, for all practical purposes the Gauss-Newton iteration converges in one iteration.We believe this is because the initial guess q 0 is already close to q true , and the Jacobian of Q • M is in some sense close to the identity, as explained in §5.2.We observe that q 1 seems to be a better reconstruction than q 0 .For example in figure 7, q 0 has artifacts close to the center due to the Delaunay triangulation used for the linear interpolation.These artifacts are diminished in q 1 .
We include a typical convergence history in figure 8. Since the Gauss-Newton iterations are designed to minimize the residual (24), we expect that the norm of the preconditioned residual Q(M (q k )) − Q(M (q true )) 2  2 decreases with k.This is true for the first iteration, but then the residual stagnates.To explain this, consider a vector z k = 0 spanning the null space of the Jacobian of Q(M (q)) evaluated at q = q k .An explicit formula for such a vector is given in proposition 5.1.Clearly the Gauss-Newton update q k+1 − q k defined in (25) is independent of the component of the preconditioned residual Q(M (q k )) − Q(M (q true )) in the direction z k .In other words, the update q k+1 − q k = −DM [q k ] † DQ[M (q k )] † (Q(M (q k )) − Q(M (q true )) + αz k ), is independent of α.We believe the stagnation of the preconditioned residual is because the component of the preconditioned residual along z k is not zero, and the component orthogonal to z k is still reduced by the iterations.To test this hypothesis, we include in figure 8 the norm of the projected preconditioned residual P k (Q(M (q k )) − Q(M (q true ))) 2  2 , where P k = I − z k z T k /(z T k z k ) is the projector on the subspace orthogonal to z k .We observe that this quantity does decrease with k, at least until reaching machine precision (at k = 2).For reference we also include the norm of the unpreconditioned residual, i.e.M (q k ) − M (q true ) 2 F .This quantity remains basically unchanged with the iterations.

Summary
We introduced and analyzed a novel inversion algorithm for determining absorptive Schrödinger potentials from measurements of the Dirichlet to Neumann map.The method can be viewed in the context of parametric model reduction, with reduced models having the physical meaning of resistor networks.Networks have been used successfully for the inverse conductivity problem.Here we show how to use them for Schrödinger's equation.The connection is made by a new discrete q true q 0 q 1 n = 17  generalized Liouville identity, defined on the line graph of the network.The set of nodes of the line graph, where the potential is discretized, is isomorphic to the set of edges of the network, where the conductivities are defined.

Iteration
Figure 8.A typical convergence history for the preconditioned Gauss-Newton method.We show convergence in terms of the unpreconditioned residual (green), the preconditioned residual (red) and the projected preconditioned residual (blue).
We use the discrete Liouville identity to obtain an approximate inverse of the forward map, and formulate a preconditioned Gauss-Newton iteration for reconstructing the continuum Schrödinger potential.The method is superior to traditional output least squares because it converges quickly, in one step, it is computationally inexpensive and gives good quantitative results.off-diagonal ones by enforcing conservation of currents, so we can guarantee [5] (at least in the noiseless case) that there is a unique resistor network with graph C n(n − 1)/2, n and DtN map M .
The lumping has two regularizing effects.The first is that by summing more measurements of Λ 1,q we get noise cancellation.The second is that the estimation of the resistors in a network from its DtN map is more stable if the network is smaller.

Figure 1 .
Figure 1.Graphs of class C(m, n), where n the number of boundary nodes and m is the number of radial and angular layers.The boundary nodes are shown as circles and the interior nodes as filled (black) circles.

Figure 2 .
Figure 2. A graph G in black and its line graph G in red.