User:Daviddaved/Rational difference operators and applications

From Wikibooks, open books for an open world
Jump to: navigation, search

Rational difference operators and applications

by D.V. Ingerman


  • an algorithm for solving inverse conductivity (Calderon's) problem for planar electrical networks based on representation of DN operators as quotients of small bandwidth/sparse matrices; the algorithm suggests modifying DN operator measurements to special Cauchy data measurements;
  • relationship b/w Radon/X-ray transform and DN map of a planar network, including a boundary rigidity theorem (discrete analog of the result by L. Pestov & G. Uhlmann);
  • a connection b/w DN monodromy matrices of a planar electrical network and Jones polynomial of a corresponding knot.


Statement of the problem[edit]

The following inverse problem, stated by Calderon, has potential practical applications to the medical and geophysical imaging and had received a lot of attention in the past decades, see []. This is the problem of recovering the conductivity of a body from the corresponding Dirichlet-to-Neumann operator. Given a domain with positive measurable function on it, the Dirichlet-to-Neumann operator connects the Dirichlet and Neumann boundary values of \gamma-harmonic functions

 \nabla\cdot(\gamma\nabla u) = 0, u\in H^1(\Omega),

defined on the domain. It is a pseudo-differential operator of order 1.

\Lambda_\gamma:H^{1/2}(\partial\Omega)\rightarrow H^{-1/2}(\partial\Omega)

The inverse conductivity problem is the problem of inverting the map

 \gamma \rightarrow \Lambda_{\gamma}.

In the 2D case one can make sense of the operator for measurable conductivities (vs. differentiable), (see [1]), using the Hilbert transform b/w \gamma-harmonic conjugate functions.

Dirtichlet-to-Neumann operator[edit]

The Dirichlet-to-Neumann operator is a special type of Poincaré–Steklov operator. On a manifold w/boundary it is the pseudo-differential operator from the Dirichlet boundary values (potential) to the Neumann boundary values (current) of harmonic functions. It is well-defined because of uniqueness and existence of the solution of the Dirichlet problem.

Definition of the Dirichlet-to-Neumann operator for a domain

The inverse conductivity problem is the problem of inverting the map

 \gamma \rightarrow \Lambda_{\gamma}.

In the 2D case one can make sense of the operator for measurable conductivities (vs. differentiable), (see [1]), using the Hilbert transform b/w \gamma-harmonic conjugate functions. The dimension of the unknown parameter fit the dimension of the measurement data in 2D. The ideas of the solutions of the inverse problems are presented in a way bridging the discrete and continuous domain settings.

Discretization technique

The techniques of solutions are motivated by making the following diagram commutative:

Discrete and Contiuous Inverse problems.jpg

Electrical networks[edit]

A planar electrical network can be defined as a weighted planar graph with natural boundary

G(V, \partial G, E, \gamma),\partial G\subset V.

The weight function defined on the edges of the network is called conductivity. An important object containing information about a network G(V,E,w) is its Laplacian matrix/Kirchhoff matrix. The boundary measurements for many inverse problems can be expressed conveniently in terms of the matrix blocks. Given a numbering of the vertices of the graph, it's an n by n square matrix LG, where n is the number of vertices in the graph, with the entries:

\sum_{v_k \rightarrow v_l} w_{kl} \mbox{  if}\ k = l, \\
-w_{kl} \mbox{  if}\ v_k \rightarrow v_l, \\
0,  \mbox{  otherwise,}

where vkvl means that there is a directed edge from vertex vk to the vertex vl, and where w is the weight function.

Given a network with boundary it is often convenient to number its boundary vertices first and to write its Laplacian matrix in the block form.

L_G = 
A & B \\
C & D

Exercise (*). Prove that a function/vector u is harmonic at the interior nodes of a graph G if

u|_{int G} = -D^{-1}Cu|_{\partial G}.


M=\begin{pmatrix} A & B \\ C & D \end{pmatrix}

be a block matrix with an invertible square block D.

Then the Schur complement of the matrix M w/respect of the non-singular block D is the matrix

 M/D = A-BD^{-1}C.

Exercise (*). Prove the following determinant identity for a square block matrix M:

\det M = \det(M/D) \det D.

The following matrix W(G) consisting of random walk exiting probabilities (sums over weighted paths in a graph) plays an important role as boundary measurement for inverse problems. Suppose a weighted graph G has N boundary nodes, then the kl 'th entry of the N by N matrix equals to the probability that the first boundary vertex, that a particle starting its random walk at the boundary vertex vk occupies, is the boundary vertex vl. For a finite connected graph the columns of the matrix W(G) add up to 1.

Exercise (**). Derive an explicit formula for the matrix W(G) in terms of the blocks of Laplace matrix L(G) of the graph G.

W_G = I - D_A^{-1}(A - BD^{-1}C).

Exercise (***). Prove the following expansion formulas for entries and blocks of the matrix W(G ),

  • for two boundary vertices pk and pl of a graph G
w_{kl} = \sum_{v_k\xrightarrow[]{path} v_l}\prod_{e=(p,q)\in path}w(e)/l_{pp},

Exercise (*). Prove that the Dirichlet-to-Neumann operator of an electrical network is symmetric and equal to the Schur complement of its Kirchhoff matrix.

\Lambda_G = 
 A & B \\
 B^T  & C 
/ (C) =  A-BC^{-1}B^T.

The Laplace equation gives the direct connection between the exiting probability of the random particle started at the boundary and the value of a harmonic function at a vertex/point, see [10]. The connection can be expressed using the sum of the geometric series identity applied to the blocks of the Kirchhoff matrix of the network/graph.

 \Lambda_G= D_A(I-W_G) = A-B^T D_C \sum_k(I-D_C^{-1}C)^k B.

This is a special case of the convergent geometric Neumann series applied to the diagonally dominated matrix.

Effective conductivities[edit]

The effective conductivity b/w two distinct nodes of an electrical network G is equal to the total current flowing b/w the nodes due to unit difference of potentials applied @ the two nodes, if there are no sources or sinks of the electrical current @ the other nodes of the network. The effective conductivities can be calculated in terms of the Schur complement of the Kirchhoff matrix of the network. The map from conductivities to effective conductivities of a complete network is non-linear and essentially an involution! The effective conductivities between all pairs of nodes of a subset ∂G of a network G determine the DN operator of {G,∂G} and vice versa.

Rotation invariant networks[edit]

If a network is rotation (sphere) or shift (half-plane) invariant then its Dirichlet-to-Neumann operator is diagonal in Fourier coordinates (since solutions of Laplace equation commute w/rotations and shifts) and is given by Stieltjets continued fraction.

\Lambda = z\beta(z)(\sqrt{-\Delta})

Continued fractions[edit]

Let \{a_k\} be a set of n positive numbers. The Stieltjes continued fraction is an expression of the form, see also [],

\beta_a(z) = a_nz + \cfrac{1}{a_{n-1}z + \cfrac{1}{ \ddots + \cfrac{1}{a_1 z} }}

or its reciprocal

\beta^*_a(z) = \frac{1}{\beta_a(z)}.

The function defines a rational n-to-1 map of the right half of the complex plane onto itself,

\beta_a,1/\beta_a:\mathbb{C^+}\xrightarrow[]{n\leftrightarrow 1}\mathbb{C^+}.

The function \beta_a is determined by its symmetric, interlacing zeros and poles. It is also determined by the pre-image set \{\mu_k\} of the point {z = 1}, since

\beta_a(z) = \frac{p(z^2)}{zq(z^2)}=1 \iff p(z^2)-zq(z^2) = 0,

and a complex polynomial is determined by its roots up to a multiplicative constant by the fundamental theorem of algebra.

If a network is rotation (sphere) or shift (half-plane) invariant then its Dirichlet-to-Neumann operator is diagonal in Fourier coordinates (since solutions of Laplace equation commute w/rotations and shifts) and is given by Stieltjets continued fraction.

\Lambda_{\gamma} = \sqrt{-\Delta}\beta_{\gamma}(\sqrt{-\Delta})

Therefore, one can reduce the layered inverse problem to a Pick-Nevanlinna interpolation problem, since given the Dirichlet-to-Neumann operator of a layered network, one can find its eigenvalues and the interpolating rational function. The conductivities of the edges of the network are then given by the coefficients of continued fraction and their reciprocals and can be found by the Eucledian algorithm.

Monodromy of difference operator[edit]

An alternative solution of the discrete inverse problem w/out use of the Pick-Nevanlinna interpolation can be based on the generalization of the calculation of the zeros and poles of the Stieltjes continued fraction for the example of the following network.

Monodromy matrix for a rotation invariant network w/4 circles and 9 radii

Since continued fraction is a rational function it can be written as a ratio of two polynomials

\beta(z) = \frac{zp(z^2)}{q(z^2)},

therefore, from the above formulas, the Dirichlet-to-Neumann operator can be written as a quotient of two high order difference operators. The zeros and poles of the continued fraction can be then interpreted as modified eigenvalues of the monodromy matrices of the difference operators.

The four eigenvalues of the monodromy matrix

M_G = \Lambda_{G}((2,3,4,5);(6,7,8,9))\Lambda_{G}((1,2,3,4);(6,7,8,9))^{-1},

are positive, come in reciprocal pairs and

q(-(\sqrt{\lambda}-\frac{1}{\sqrt{\lambda}})^2) = 0 \iff \lambda\in\sigma(M_G),


\beta_{\gamma}(z) = \cfrac{1}{\cfrac{1}{\gamma_4}z + \cfrac{1}{\gamma_{3}z+\cfrac{1}{\cfrac{1}{\gamma_2}z +\cfrac{1}{\gamma_1 z}}}}.

Exercise(**). Find the similar formula for the zeros of the continued fraction.

Exercise(*). Prove that

\det M_G = 1.

Exercise(**). Prove the following trace formula:

\frac{1}{2}Tr((\sqrt{M_G}-\sqrt{M_G}^{-1})^2)= \frac{1}{2}\sum_{k=1}^4 (\sqrt{\lambda_k}-\frac{1}{\sqrt{\lambda_k}})^2 = \frac{\gamma_2}{\gamma_1}+\frac{\gamma_2}{\gamma_3}+\frac{\gamma_4}{\gamma_3}.

Exercise(**). Prove the following determinant formula:

\det(\sqrt{M_G}-\sqrt{M_G}^{-1})= \prod_{k=1}^4 (\sqrt{\lambda_k}-\frac{1}{\sqrt{\lambda_k}}) = \frac{\gamma_2 \gamma_4}{\gamma_1 \gamma_3}.

Generalization to non-layered domain[edit]

General 2D inverse problem[edit]

The contimued fraction approach can be applied to domains w/inhomogenuous isotropic conductivity case in 2D and higher dimensions using layers discretization. The resulting pseudo-differential Dirichlet-to-Neumann operator can be written as a ratio of two high order differential operators that satisfy three-term recurrence (similar to the numerators and denominators of functional continued fraction). The fundumental solutions of the differential operators can be directly read from the kernel of the Dirichlet-to-Neumann operator. The conductivity can be then found by a Eucledian type algorithm, reversing the three-term recurrence, which is equivalent to solving inverse conductivity (Calderon's) problem for a planar network.

More formally: Let f be the potential and g be the current on a layered domain. Then,

g_0 = D\gamma_0 D f_0, \\
f_{k+1}= f_k + \alpha_k g_k,\\
g_{k+1} = g_k + D\gamma_k Df_k.

\mbox{ In the dual domain, }

g_0 = \gamma_0 (f_0 - \frac{\int\gamma_0 f_0}{\int\gamma_0}), \\
f_{k+1}= f_k + \alpha_k g_k,\\
g_{k+1} = g_k + D\gamma_k Df_k.

Therefore, for ordinary differential operators F and G,

f_k=F_k f_0, \\
g_k=G_k f_0, \\
\Lambda_k = G_k F_k^{-1} = (F_k^*)^{-1}G_k^*,

since, the Dirichlet-to-Neumann operator is self-adjoint. And

F_k^* G_k = G_k^* F_k, \\
G_k^*f_k = F_k^*g_k, \\
G_k = \Lambda_k F_k.

A non-obvious fact: the determinants of the monodromy matrices of the operators F and G are equal to 1. The eigenvalues of the matrices are positive and interlace. They come in the reciprocal pairs in the rotational invariant case.

Since, the operator G is differential, the support of the function Gf belongs to the closure of the support of the potential f, which allows one to read the solutions of the equation F*g=0 from the Dirichlet-to-Neumann operator and G*f=0 from the dual problem. (Need to find a family of fundamental solutions). The coefficients of the differential operators F, G, F* and G* can be then obtained from the Wronskians of the fundamental solutions.

One also gets from the systems above that,

F_{k+1} = \alpha_0\gamma_0\ldots,\alpha_k\gamma_kD^{k+1}+\ldots \mbox{ and } \\
G_{k+1} = \alpha_0\gamma_0\ldots,\alpha_{k-1}\gamma_{k}D^{k+2}+\ldots.

Therefore, one can find the conductivity of the outmost layer from as ratio of the leading coefficients of the corresponding operators F an G, which allows one to reverse the three-term recurrence and find all the conductivities by induction.

Boundary rigidity[edit]

The following construction plays an important role in studying the graphs properly embedded into surfaces. One considers an embedded into a surface (2D manifold) finite graph w/vertices each w/exactly four neighbors. The union of edges of such graph is equal to a union of closed curves on the surface, called geodesics. The faces of the graph can be two-colored.

A medial graph (morphed chess board)

The corresponding graph G and its dual G* can be obtained by putting vertices in all faces of the same color and connecting two vertices by an edge if the corresponding faces of the medial graph M(G) have a common corner. Note that M(G)=M(G*).

One cam modify the construction accordingly for the case of an embedded graph w/boundary by allowing vertices of the medial graph to have one neighbor.

The following transformations of the geodesics of the medial graph M(G) are similar to the Reidemeister moves of the knot theory. They correspond to the Y-Δ transformations of the graph G and G*.

The three moves of geodesics curves

The following matrix consists of the distances (the smallest number of jumps across the geodesics or the edge distance in M*(G)) between boundary faces of the medial graph below.

 0 & 1  & 2 & 3 & 2 & 3 & 2 & 1 \\
 1 & 0  & 1 & 2 & 3 & 4 & 3 & 2 \\
 2 & 1  & 0 & 1 & 2 & 3 & 4 & 3 \\
 3 & 2  & 1 & 0 & 1 & 2 & 3 & 4 \\
 2 & 3  & 2 & 1 & 0 & 1 & 2 & 3 \\
 3 & 4  & 3 & 2 & 1 & 0 & 1 & 2 \\
 2 & 3  & 4 & 3 & 2 & 1 & 0 & 1 \\
 1 & 2  & 3 & 4 & 3 & 2 & 1 & 0
Two examples of boundary isometric medial graphs with 4 and 5 geodesics

The distances between boundary faces of an embedded medial graph are obviously invariant under the geodesics moves.

The geodesics of the medial graph form an arrangement of pseudo-lines if

  • none of them intersect itself
  • none of them is a closed curve
  • no two of them intersect more than once

Looking for the pattern

k+1 & k \\
k & k+1

in the boundary distances matrix one can prove the following boundary rigidity theorem:

(a) a planar medial graph M(G) is determined, up to the moves, by the distances between its boundary faces; (b) every planar medial graph M(G) can be moved into an arrangement of pseudo-lines without changing the distance matrix.

The following min-cut max-flow result was formulated w/help of Derek Jerina, a student of REU summer school on inverse problems at the University of Washington.

Exercise (**). Prove that the ranks of submatrices of the Dirichlet-to-Neumann operator of a planar electrical network w/natural boundary determine the shortest distances between the boundary faces of its medial graph, and therefore determine the planar graph up to a finite sequence of Y-Δ transformations.

(Hint.) It follows from the determinant identity in the section Special matrices that the ranks of the submatrices count the numbers of disjoint paths connecting boundary nodes of the planar graph.

Numerical example[edit]

Knot of a planar electrical network[edit]


The author would like to thank V. Druskin for fruitful discussions.