A discrete Gelfand - Levitan theory

F. Natterer Institut für Numerische und instrumentelle Mathematik Westf. Wilhelms-Universität Münster Einsteinstrasse 62 D-4400 Münster, West-Germany

September 1989

1 Introduction

In the pioneering paper [5], Gelfand and Levitan introduced an elegant and explicit method for computing the potential of a Sturm-Liouville Operator from its eigenvalues and certain values of its eigenfunctions. It has been noted in Burridge [3] that the Gelfand-Levitan method is closely related to inversion methods for recovering the potential in a hyperbolic equation with focused initial state, thus giving a unified theory for the methods of Gelfand-Levitan, of Gopinath and Sondhi [8] and of Parijskij [11] and Blagoveshchenskij [1]; (see also Romanov [12, p. 39]), the latter ones going unnoticed in [3].

In the present note we analyse discrete analogues of these problems. The spacial differential operator of the Gelfand-Levitan theory is replaced by a symmetric tridiagonal matrix. It will turn out that in this setting the Gelfand-Levitan method reduces essentially to a Cholesky decomposition. The hyperbolic equation in the other theories is replaced by a recursion relation involving an arbitrary tridiagonal matrix, the essential step for the inversion being an LU-decomposition. In our approach, (approximate) transmutation operators as originally used in [5] (see Levitan [10] for a systematic study) play a paramount role.

We are aware of the fact that much work has been done on discrete inverse problems. We mention in particular Burridge [3, p. 514-537], Case and Kack [4], Bruckstein and Kailath [2], Landau [9], Gladwell and Willms [7]. Our approach differs from others in that it gives a unified treatment in terms of transmutation operators.

2 Transmutation matrices

Let T, T_{0} be tridiagonal (n, n)-matrices, i.e.

T = ( \begin{matrix}α_{1} & β_{1} \\ γ_{1} & α_{2} & β_{2} \\ . & . & . \\ & & & & β_{n−1} \\ & & & γ_{n−1} & α_{n} \end{matrix} ) , T_{0} = ( \begin{matrix}α_{1}^{0} & β_{1}^{0} \\ γ_{1}^{0} & α_{2}^{0} & β_{2}^{0} \\ . & . & . \\ & & & & β_{n−1}^{0} \\ & & & γ_{n−1}^{0} & α_{n}^{0} \end{matrix} ) .

By a transmutation matrix for T, T_{0} one usually means a (n, n)-matrix L with T_{0}L = LT. We shall use (approximate) transmutation matrices for which this holds only in rows 1 through n − 1. More precisely, with E the projection

E = ( \begin{matrix}1 \\ & 1 \\ & & . & . & . \\ & & & & & 1 \\ & & & & & & 0 \end{matrix} )

we define a transmutation L by

E T_{0} L = E L T .  (2.1)

Theorem 2.1 Let β_{i}^{0} ≠ 0, i = 1, …, n − 1. Then there is a unique lower triangular transmutation L = (ℓ_{ij}) whose (1, 1)-entry is 1.

Proof: With ℓ_{1}, …, ℓ_{n} the rows of L, (2.1) reads

\begin{aligned}α_{1}^{0}ℓ_{1} + β_{1}^{0}ℓ_{2} & = & ℓ_{1}T , \\ γ_{i−1}^{0}ℓ_{i−1} + α_{i}^{0}ℓ_{i} + β_{i}^{0}ℓ_{i+1} & = & ℓ_{i}T , i = 2, …, n − 1 . \end{aligned}

Since β_{i}^{0} ≠ 0 this recursion determines the ℓ_{i} uniquely once ℓ_{1} is given. For ℓ_{1} = (1, 0, …, 0), ℓ_{i} is nonzero only in its first i components since T is tridiagonal.

□

We remark that Theorem 2.1 holds also for lower Hessenberg matrices T_{0}, T.

3 The inverse eigenvalue problem

Now assume T to be symmetric with non-zero off-diagonal elements, i.e. γ_{i} = β_{i} > 0, i = 1, …, n − 1. Let λ_{1}, …, λ_{n} be the eigenvalues and x(λ_{1}), …, x(λ_{n}) the normalized eigenvectors of T, i.e.

T x(λ_{k}) = λ_{k}x(λ_{k}) , k = 1, …, n .

We consider the inverse problem: Determine T from λ_{1}, …, λ_{n} and x_{1}(λ_{1}), …, x_{1}(λ_{n}). It is well known that this problem can easily be solved by the Lanczos method.

The Lanczos algorithm computes for a symmetric matrix A a symmetric tridiagonal matrix T and a unitary matrix U such that

T U = UA  (3.1)

in the following way. With u_{1}, …, u_{n} the rows of U, (3.1) reads

α_{1}u_{1} + β_{1}u_{2}  = u_{1}A
β_{1}u_{1} + α_{2}u_{2} +β_{2}u_{3}  = u_{2}A
...
β_{n−1}u_{n−1} + α_{n}u_{n} = u_{n}A

Now let u_{1} be arbitrary. Then, from equation 1,

α_{1} = u_{1}Au_{1}^{∗} ,
ũ_{2} = u_{1}A − α_{1}u_{1} , β_{1} = ∥ũ_{2}∥ (euclidean norm)
u_{2} = ũ_{2}/β_{1} .

Thus α_{1}, β_{1}, u_{1}, u_{2} are determined. Likewise, from equation 2,

α_{2} = u_{2}Au_{2}^{∗}
ũ_{3} = u_{2}A − β_{1}u_{1} − α_{2}u_{2} , β_{2} = ∥ũ_{3}∥
u_{3} = ũ_{3}/β_{2} ,

yielding β_{2}, u_{3}. Proceeding in this fashion in equations up to and including n − 1 we obtain α_{1}, …, α_{n−1}, β_{1}, …, β_{n−1} and u_{1} … u_{n}. α_{n} is obtained from equation n.

In order to solve our inverse problem we simply apply the Lanczos method to the matrix A = diag(λ_{1}, …, λ_{n}). The same problem can be solved by the Gelfand-Levitan method in the following way.

Let T_{0} be an arbitrary known symmetric tridiagonal matrix, subject only to the condition β_{i}^{0} ≠ 0, i = 1, …, n − 1. We introduce solutions φ(λ), φ^{0}(λ) ∈ R^{n} of

\begin{aligned}φ_{1}(λ) & = & 1 , \\ α_{1}φ_{1}(λ) + β_{1}φ_{2}(λ) & = & λφ_{1}(λ) , \\ β_{i−1}φ_{i−1}(λ) + α_{i}φ_{i}(λ) + β_{i}φ_{i+1}(λ) & = & λφ_{i}(λ) , i = 2, …, n − 1 \end{aligned}

and correspondingly for φ^{0} using T_{0} instead of T. In other words, φ(λ), φ^{0}(λ) satisfy

ET φ(λ) = λEφ(λ) , φ_{1}(λ) = 1  (3.2)
ET_{0}φ^{0}(λ) = λEφ^{0}(λ) , φ_{1}^{0}(λ) = 1 .  (3.3)

Note that for λ = λ_{k}, φ(λ_{k}) is an eigenvector of T, hence

T φ(λ_{k}) = λ_{k}φ(λ_{k}) , k = 1, …, n ,  (3.4)

i.e. the projection operator E can be dropped in that case.

Theorem 3.1 Let β_{i}^{0} ≠ 0, i = 1, …, n − 1 and let L be the transmutation of Theorem 2.1 Then,

φ^{0}(λ_{k}) = Lφ(λ_{k}) , k = 1, …, n .

Proof: Let ψ^{0} = Lφ(λ_{k}). Then, because of (2.1), (3.4),

ET_{0}ψ^{0} = ET_{0}Lφ(λ_{k}) = ELT φ(λ_{k}) = λ_{k}ELφ(λ_{k}) = λ_{k}Eψ^{0} .

Since ℓ_{11} = 1, ψ_{1}^{0} = 1. Thus ψ^{0} and φ^{0}(λ_{k}) both satisfy (3.3) with λ = λ_{k}. Because of β_{i}^{0} ≠ 0, (3.3) is uniquely solvable. Hence ψ^{0} = φ^{0}(λ_{k}).

□

Now we come to the core of the Gelfand-Levitan method. We introduce the matrices

Φ = (φ(λ_{1}), …, φ(λ_{n})) , Φ_{0} = (φ^{0}(λ_{1}), …, φ^{0}(λ_{n})) ,

P = ( \begin{matrix}p_{1} \\ & p_{2} \\ & & . & . & . \\ & & & & & p_{m} \end{matrix} ) , p_{k} = ∥φ(λ_{k})∥ (euclidean norm) .

Assume that x_{1}(λ_{k}) > 0, k = 1, …, n. Then, φ(λ_{k}) = \frac{1}{x_{1}(λ_{k})} x(λ_{k}), hence p_{k} = \frac{1}{x_{1}(λ_{k})} and

Φ = XP , X = (x(λ_{1}), …, x(λ_{n})) .

Since XX^{t} = I,

ΦP^{−2}Φ^{t} = I .  (3.5)

Theorem 3.1 means that

Φ_{0} = LΦ .  (3.6)

Inserting this into (3.5) yields

\begin{matrix}L^{−1}Φ_{0}P^{−2}Φ_{0}^{t}L^{−t} = I , \\ LL^{t} = Φ_{0}P^{−2}Φ_{0}^{t} . \end{matrix}

Since Φ_{0}, P are known from the data, L can be computed by a Cholesky decomposition of Φ_{0}P^{−2}Φ_{0}^{t}. Once L is determined, T can be computed from

T = ΦΛΦ^{−1} = L^{−1}Φ_{0}ΛΦ_{0}^{−1}L , Λ = ( \begin{matrix}λ_{1} \\ & . & . & . \\ & & & & λ_{n} \end{matrix} ) .

This solves the inverse eigenvalue problem.

4 The inverse evolution problem

Now let T = T_{0} + Q with T_{0} known and Q a diagonal matrix whose entries we denote by q_{1}, …, q_{n}. Let u^{ℓ} be a solution to

u^{ℓ+1} = T u^{ℓ} , ℓ = 0, 1, …,  u^{0} = e_{1}  (4.1)

with e_{1} the first unit vector. We consider the inverse problem of recovering Q from

u_{1}^{ℓ} = g^{ℓ} , ℓ = 0, …, 2n − 1 .

Again we make use of the transmutation L from Theorem 2.1. Multiplying (4.1) with EL we obtain from (2.1)

ELu^{ℓ+1} = ELT u^{ℓ} = ET_{0}Lu^{ℓ} .

Thus the vectors z^{ℓ} = Lu^{ℓ} satisfy

Ez^{ℓ+1} = ET_{0}z^{ℓ} ,  ℓ = 0, 1, …,  (4.2)

i.e. z^{ℓ+1} = T_{0}z^{ℓ} in components 1, …, n − 1. Since T_{0} is tridiagonal we can show that for k = 0, …, n − 1

z_{i}^{ℓ+k} = (T_{0}^{k}z^{ℓ})_{i} ,  i = 1, …, n − k .  (4.3)

We use induction with respect to k. The case k = 0 is obvious. Assume that (4.3) is correct for some k with 0 ≤ k < n − 1. From (4.2) we get

z_{i}^{ℓ+k+1} = (T_{0}z^{ℓ+k})_{i} , i = 1, …, n − 1 .

For the evaluation of (T_{0}z^{ℓ+k})_{i} , i = 1, …, n − k − 1, we need only the first n − k components of z^{ℓ+k} since T_{0} is tridiagonal. Thus, by the induction hypothesis,

z_{i}^{ℓ+k+1} = (T_{0}T_{0}^{k}z^{ℓ})_{i} = (T_{0}^{k+1}z^{ℓ})_{i} , i = 1, …, n − k − 1 .

This is (4.3) for k + 1. Hence (4.3) is established.

Since the first row of L is (1, 0, …, 0) we also have

z_{1}^{ℓ} = g^{ℓ} , ℓ = 0, 1, ….

This combines with (4.3) to yield, for ℓ = 0, 1, …, n − 1,

g^{ℓ+k} = (T_{0}^{k}z^{ℓ})_{1} , k = 0, 1, …, n − 1 .

Introducing the row vector

(u_{0}^{k})^{t} = e_{1}^{t}T_{0}^{k}

we have

g^{ℓ+k} = (u_{0}^{k})^{t}z^{ℓ} , k = 0, …, n − 1 (4.4)

and

u_{0}^{k+1} = T_{0}^{t}u_{0}^{k} , k = 0, 1, …, u_{0}^{0} = e_{1} .  (4.5)

With the (n, n)-matrices

Z = (z^{0}, …, z^{n−1}) , G = ( \begin{matrix}g^{0} & g^{1} & … & g^{n−1} \\ ⋮ & ⋮ \\ g^{n−1} & g^{n} & … & g^{2n−2} \end{matrix} ) , U_{0} = (u_{0}^{0}, …, u_{0}^{n−1})

(4.4) simply reads

G = U_{0}^{t}Z .

Note that U_{0} is upper triangular with diagonal elements

∏_{i=1}^{k} β_{i}^{0} , k = 0, …, n − 1 .

Thus U_{0} is invertible provided that β_{i}^{0} ≠ 0, i = 1, …, n − 1. Hence

Z = U_{0}^{−t}G

which determines Z by the data.

The relations z^{ℓ} = Lu^{ℓ} can be written as

Z = LU , U = (u^{0}, …, u^{n−1}) .

We finally obtain

LU = U_{0}^{−t}G .

Since T is tridiagonal and u^{0} = e_{1}, U is upper triangular, its diagonal entries being

∏_{i=1}^{k} γ_{i} = ∏_{i=1}^{k} γ_{i}^{0} , k = 0, …, n − 1 .

Thus L, U can be determined simply by doing an LU-decomposition on the matrix U_{0}^{−t}G, with the diagonal of U being known.

Once L, U are known there is a variety of ways to find Q. For instance we can compute u^{n} from Lu^{n} = z^{n} and

( g^{n} ; ⋮ ; g^{2n−1} ) = U_{0}^{t} z^{n} ,

which is (4.4) for ℓ = n. Then,

u^{ℓ+1} = T_{0}u^{ℓ} + ( q_{1}u_{1}^{ℓ} ; ⋮ ; q_{ℓ+1}u_{ℓ+1}^{ℓ} ; 0 ; ⋮ ; 0 ) , ℓ = 0, …, n − 1

from which q_{1}, …, q_{n} can be computed recursively since u_{ℓ+1}^{ℓ} ≠ 0 , ℓ = 0, …, n − 1. This solves the inverse evolution problem.

5 A second order inverse evolution problem

Let T be as in the preceeding section, and let u^{ℓ} be a solution to

u^{ℓ+1} − 2u^{ℓ} + u^{ℓ−1} = T u^{ℓ} , ℓ = 0, 1, 2, …
u^{−1} = 0 , u^{0} = e_{1} .  (5.1)

We consider the inverse problem of recovering the diagonal matrix Q in T = T_{0} + Q from

u_{1}^{ℓ} = g^{ℓ} , ℓ = 0, 1, …, 2n .

Again we put z^{ℓ} = Lu^{ℓ} with L from Theorem 2.1. We obtain

Ez^{ℓ+1} − 2Ez^{ℓ} + Ez^{ℓ−1} = ET_{0}z^{ℓ} , ℓ = 0, 1, 2, …
z^{−1} = 0 , z^{0} = Le_{1} .

With S = 2I + T_{0} the recursion can be written as

z_{i}^{ℓ+1} = (Sz^{ℓ})_{i} − z_{i}^{ℓ−1} ,  i = 1, …, n − 1 .  (5.2)

We claim that for k = 0, …, n − 1 and ℓ = 0, 1, 2, …

z_{i}^{ℓ+k} = (S_{k}z^{ℓ})_{i} − (S_{k−1}z^{ℓ−1})_{i} , i = 1, …, n − k  (5.3)

where

S_{−1} = 0 , S_{0} = I , S_{k+1} = SS_{k} − S_{k−1} , k = 0, …, n − 1 .

We prove (5.3) by induction with respect to k. The case k = 0 is obvious. Assume (5.3) to be correct up to some k with 0 ≤ k < n − 1. From (5.2) we get

z_{i}^{ℓ+k+1} = (Sz^{ℓ+k})_{i} − z_{i}^{ℓ+k−1} , i = 1, …, n − k − 1 .

Since S is tridiagonal, only components 1 through n−k of z^{ℓ+k} enter (Sz^{ℓ+k})_{i} for i < n − k. Hence, by the induction hypothesis,

\begin{aligned}z_{i}^{ℓ+k+1} & = & (S(S_{k}z^{ℓ} − S_{k−1}z^{ℓ−1}))_{i} − (S_{k−1}z^{ℓ} − S_{k−2}z^{ℓ−1})_{i} \\ & = & (S_{k+1}z^{ℓ})_{i} − (S_{k}z^{ℓ−1})_{i} , i = 1, …, n − k − 1 . \end{aligned}

This is (5.3) with k replaced by k + 1. Thus (5.3) is established.

We use (5.3) for i = 1 only, yielding for ℓ = 1, 2, …

g^{ℓ+k} = (S_{k}z^{ℓ})_{1} − (S_{k−1}z^{ℓ−1})_{1} , k = 0, …, n − 1 .

Introducing the row vectors

(u_{0}^{k})^{t} = e_{1}^{t}S_{k} , k = −1, …, n − 1

we have

g^{ℓ+k} = (u_{0}^{k})^{t}z^{ℓ} − (u_{0}^{k−1})^{t}z^{ℓ−1} , k = 0, …, n − 1

and, since S, S_{k} commute,

\begin{aligned}(u_{0}^{k+1})^{t} & = & e_{1}^{t}S_{k+1} = e_{1}^{t}(SS_{k} − S_{k−1}) \\ & = & e_{1}^{t}(S_{k}S − S_{k−1}) = (u_{0}^{k})^{t}S − (u_{0}^{k−1})^{t} . \end{aligned}

Thus u_{0}^{k} is the analogue to u^{k} with T replaced by T_{0}^{t}:

u_{0}^{k+1} − 2u_{0}^{k} + u_{0}^{k−1} = T_{0}^{t}u^{k} k = 0, 1, 2, …
u_{0}^{−1} = 0 , u_{0}^{0} = e_{1} .

This shows that the matrix U_{0} = (u_{0}^{0}, …, u_{0}^{n−1}) is upper triangular with diagonal elements

∏_{i=0}^{k} β_{i}^{0} , k = 0, …, n − 1 .

It follows that the systems

( g^{ℓ} ; ⋮ ; g^{ℓ+n−1} ) = U_{0}z^{ℓ} − U_{−1}z^{ℓ−1} , ℓ = 1, …, n

where U_{−1} = (u_{0}^{−1}, …, u_{0}^{n−2}) are uniquely solvable for z^{ℓ} provided that β_{i}^{0} ≠ 0, i = 1, …, n − 1.

Given that z^{−1} = 0, we can determine the matrix Z = (z^{0}, …, z^{n−1}) from our data. The matrix U = (u^{0}, …, u^{n−1}) is upper triangular with diagonal elements

∏_{i=1}^{k−1} γ_{i} = ∏_{i=1}^{k−1} γ_{i}^{0} , k = 1, …, n .

Thus computing L, U amounts to doing an LU-Decomposition on the matrix Z,

Z = LU ,

with the diagonal of U being known. Once L, U are known, Q is computed very much in the same fashion as in the preceding section.

6 Links to other work

It has been noticed in Landau [9] that the inverse eigenvalue problem is intimately related to orthogonal polynomials. Our treatment in section 2 can also be interpreted that way.

To begin with it is clear that φ_{k}(λ), φ_{k}^{0}(λ) are polynomials of degree k −1 in λ. Thus (3.5) simply states that the polynomials φ_{k}, k = 1, …, n, are the orthonormal polynomials with respect to the scalar product

(p, q) = ∑_{k=1}^{n} p_{k}^{−2}p(λ_{k})q(λ_{k})  (6.1)

in the space of polynomials of degree n − 1. Writing φ_{k}^{0} in terms of the orthogonal polynomials we get

φ_{ℓ}^{0} = ∑_{k=1}^{ℓ} ℓ_{ℓk}φ_{k}

or

Φ_{0} = LΦ

with a lower triangular matrix L. This is (3.6). Thus solving the inverse eigenvalue problem, i.e. finding L, boils down to computing the orthogonal polynomials with respect to the scalar product (6.1).

The approach of Gladwell and Willms [7] is essentially equivalent to our treatment of the inverse eigenvalue problems, except that we made the role of transmutation more explicit. For instance, equation (10) of [7] is precisely our relation (3.6) in the form Φ = L^{−1}Φ_{0}.

References

[1] Blagoveshchenskij, A. S. (1971): On the local method of solving the non-stationary inverse problem for an inhomogeneous string (Russian). Trudy (Mat. inst. imeni V.A. Steklova, Leningrad) 65, 28-38.

[2] Bruckstein, A. - Kailath, T (1986): Some matrix factorization identities for discrete inverse scattering. Linear Algebra and its Applications 74, 157-172.

[3] Burridge, R. (1980): The Gelfand-Levitan, the Marchenko, and the Gopinath-Sondhi integral equations of inverse scattering theory, regarded in the context of inverse impulse-response problems. Wave Motion 2, 305-323.

[4] Case, K. M. - Kack, M. (1983): A discrete version of the inverse scattering problem. J. Math. Phys. 4, 594-603.

[5] Gelfand, I. M. - Levitan, B. M. (1951): On the determination of a differential equation from its spectral function. Izv. Akad. Nauk SSSR 15, 309-360 (1951) (Russian). Engl. Transl. AMS Translations 1, 253-304 (1956).

[6] Gladwell, G. M. L. (1986): Inverse problems in vibration. Martinus Nijhoff Publishers, Dordrecht - Boston - Lancaster.

[7] Gladwell, G. M. L. - Willms, N. B. (1989): A discrete Gel‘fand-Levitan method for band-matrix inverse eigenvalue problems. Inverse Problems 5, 165-179.

[8] Gopinath, B. - Sondhi, M. M. (1971): Determination of the vocal-tract shape from impulse response at the lips. J. Acoust. Soc. Amer. 49, 1867-1873.

[9] Landau, H. J. (1983): The inverse problem for the vocal tract and the moment problem. SIAM J. Math. Anal. 14, 1019-1035.

[10] Levitan, B. M. (1987): Inverse Sturm-Liouville problems. VNU Science Press, Utrecht.

[11] Parijskij, B. S. (1968): The inverse problem for a wave equation with a depth effect. In: Some Direct and Inverse Problems of Seismology (Computing Seismology, issue 4) (Russian). Nauka, Moscow, 139-169.

[12] Romanov, V. G. (1987): Inverse problems of mathematical physics. VNU Science Press, Utrecht.

0LaTeX: berg f berger paperfne Gellev.tex