# Newton Method to Recover the Phase Accumulated during MRI Data Acquisition

Authors: Oh-In Kwon; Chunjae Park
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

For an internal conductivity image, magnetic resonance electrical impedance tomography (MREIT) injects an electric current into an object and measures the induced magnetic flux density, which appears in the phase part of the acquired MR image data. To maximize signal intensity, the injected current nonlinear encoding (ICNE) method extends the duration of the current injection until the end of the MR data reading. It disturbs the usual linear encoding of the MRk-space data used in the inverse Fourier transform. In this study, we estimate the magnetic flux density, which is recoverable from nonlinearly encoded MR k-space data by applying a Newton method.


---

## Body

## 1. Introduction

An electric current injected into an electrically conducting object, such as the human body, induces an internal distribution of the magnetic flux densityB=(Bx,By,Bz). Magnetic resonance electrical impedance tomography (MREIT) visualizes the internal conductivity distribution from the z-component Bz of B which can be measured in practice using an MRI scanner. This technique was originally proposed by Joy et al. in 1989 [1]; since then, several researchers [2–10] have investigated and further developed MREIT as well as magnetic resonance current density image (MRCDI), which has similar modalities [11–13].The magnetic flux densityBz induced by injecting current through the electrodes attached on the surface of a conducting object Ω accumulates its signals in the phase parts of acquired MR image data. The conventional current-injection method [1, 8] injects the current during time Tc, between the end of the first RF pulse and the beginning of the reading gradient, in order to ensure gradient linearity.Since the signal-to-noise ratio (SNR) of the MR magnitude depends on the echo timeTE, it is impossible to increase both Tc and the SNR of the MR magnitude simultaneously in order to reduce noise effects. As an attempt to reduce the noise level, the injected current nonlinear encoding (ICNE) method was developed in 2007 [14, 15]; it extends the duration of the injection current until the end of a reading gradient in order to maximize the signal intensity of Bz. Then, it disturbs the usual linear encoding of MR k-space data used in the inverse Fourier transform.For example, the one-dimensional inverse problem in the conventional acquisition method is to find the unknown discrete magnetic flux densitybl,l=0,…,N-1 from the N×N matrix A satisfying the following:
(1.1)Aρ→=S→,Akl=e-i(2π/N)(kl+Tcbl),k,l=0,…,N-1,
where Tc is a constant and ρ→,S→ are known quantities that can be measured. By an inverse Fourier transform ℱ-1, the first equation in (1.1) becomes
(1.2)diag(e-i(2π/N)Tcb0,e-i(2π/N)Tcb1,…,e-i(2π/N)TcbN-1)ρ→=ℱ-1(S→).
The unknown data bl is simply recovered from (1.2). In the ICNE method, however, the matrix A is perturbed to Akl′=e-i(2π/N)k(l+bl). Then, it becomes a system of nonlinear equations for unknown bl data, where the conventional inverse Fourier transform is no longer applicable.In this paper, we prove a unique determination of the magnetic flux density from measured MR signal obtained by the ICNE acquisition method. Secondly, applying a Newton method, we suggest a bound ofl2-norm for recoverable magnetic flux density from nonlinearly encoded MR k-space data. Numerical experiments show the feasibility of the proposed method.


## 2. ICNE Method and Invertibility

For a standard spin echo pulse sequence in MR imaging, thek-space MR signal
(2.1)S(kx,ky)=∫ℝ2‍ρ(x,y)e-iδ(x,y)ei2π(kxx+kyy)dxdy
is measured, where ρ denotes a positive spin density of the imaging slice and δ any systematic phase artifact [16]. From the signal S in (2.1), by applying the conventional inverse Fourier transform, we can obtain
(2.2)ρ0(x,y)=ρ(x,y)eiδ(x,y),
and the clinical MR image data ρ=|ρ0|.In MREIT, we inject the currentI through the electrodes attached on the three-dimensional conducting object Ω, having conductivity distribution σ. The injection current I produces the internal current density J and the magnetic flux density B=(Bx,By,Bz) in Ω, satisfying the Ampère and Biot-Savart laws. Since an MRI scanner measures only the main magnetic field direction component of B, the z-component Bz, we focus on the problem of measuring Bz(x,y)=Bz(x,y,z0), where z0 is the center of the selected imaging slice. Since MREIT is a methodology for reconstructing the internal conductivity σ from Bz data, it is important to measure Bz more precisely.


### 2.1. ConventionalBz Acquisition

For a conventionalBz acquisition, current is not injected during Ts of the MR data acquisition, ADC as shown in Figure 1. In this case, the induced magnetic flux density Bz provides additional dephasing of spins, and, consequently, extra phase is accumulated during the total injection time Tc. Then, the measured k-space data for the injection current I can be represented as follows:
(2.3)SI(kx,ky)=∫ℝ2‍ρ(x,y)eiδ(x,y)e-iγTcBz(x,y)e-i2π(kxx+kyy)dxdy,
where γ=26.75×107 rad/T·s is the gyromagnetic ratio of hydrogen.Figure 1
Conventional and ICNE current injections in a spin echo pulse sequence.From the measuredS and SI in (2.1) and (2.3), by applying an inverse Fourier transform, we obtain ρ0 in (2.2) and
(2.4)ρI(x,y)=ρ(x,y)eiδ(x,y)e-iγTcBz(x,y).
Then, the magnetic flux density Bz is precisely computed as
(2.5)Bz(x,y)=-1γTcarctan(α(x,y)β(x,y)),
where α and β are the imaginary and real parts of ρI/ρ0, respectively.


### 2.2. ICNEBz Acquisition

In the ICNEBz acquisition, in order to improve the SNR of Bz, we prolong the current injection time Tc until the end of the MR data acquisition, as shown in Figure 1. Then, since the induced Bz disturbs the linearity of the reading gradient, the measured k-space data has lost the linear encoding characteristic as
(2.6)SC(kx,ky)=∫ℝ2‍ρ(x,y)eiδ(x,y)e-iγTcBz(x,y)e-i2π[kx(x+Bz(x,y)/G)+kyy]dxdy,
where G is a constant that denotes the strength of the magnetic reading gradient. The inverse problem arising in the ICNE method is to recover Bz(x,y) from ρ0(x,y) obtained in (2.2) and the measured signal SC(kx,ky) in (2.6).Although the inversion is not uniquely solvable forBz(x,y) in general, we can uniquely determine Bz(x,y) by assuming that φ(x):=x+Bz(x,y)/G is monotone increasing in the following theorem.Theorem 2.1.
Letρ(x,y) have a finite support Ω. If Bz(x,y) is sufficiently small to guarantee that φ(x)=x+Bz(x,y)/G is monotone increasing so that φ′(x)>0 for each y, then Bz(x,y) is uniquely recovered in Ω from ρ0(x,y) in (2.2) and SC(kx,ky) in (2.6).Proof.
We note that the linear encoding characteristic in theky-variable remains unperturbed in (2.6). Thus, by one-dimensional inverse Fourier transform, SC in (2.6) is reduced to S^C in the (kx,y)-hybrid space as the following:
(2.7)S^C(kx,y)=∫-∞∞‍ρ(x,y)eiδ(x,y)e-iγTcBz(x,y)e-i2π[kx(x+Bz(x,y)/G)]dx.
Then, the ICNE inverse problem suffices to consider the x-directional inversion of Bz(x,y) from ρ0(x,y) in (2.2) and S^C(kx,y) in (2.7) for each fixed y.
By change of variables withφ(x), (2.7) is changed into
(2.8)S^C(kx,y)=∫-∞∞‍ρ(x,y)φ′(x)eiδ(x,y)e-iγTcBz(x,y)e-i2πkxφdφ.
From (2.8), by inverse Fourier transform for the φ-variable, φ satisfies
(2.9)Φ(φ)=ρ(x,y)φ′(x)eiδ(x,y)e-iγTcBz(x,y),
where Φ is a function defined with the given S^C by
(2.10)Φ(z)=∫-∞∞‍S^C(kx,y)ei2πkxzdkx.
The relation (2.9) gives us the simple ordinary differential equation as follows:
(2.11)|Φ(φ(x))|φ′(x)=ρ(x,y).
Since ρ has a finite support, for each y, we can define
(2.12)x0=min{x:ρ(x,y)≠0},β=min{z:Φ(z)≠0}.
If φ(x0)<β, then φ(x0+ϵ)<β for a sufficiently small ϵ>0. It contradicts (2.11), since it implies that
(2.13)0=|Φ(φ(x0+ϵ))|φ′(x0+ϵ)=ρ(x0+ϵ,y)≠0.
By the same argument, the reverse inequality is not possible. Thus, we have
(2.14)ϕ(x0)=β.
By separation of variables, (2.11) and (2.14) lead us to
(2.15)∫βφ(x)‍|Φ(φ)|dφ=∫x0x‍ρ(x,y)dx.
For any given x, φ(x) is uniquely determined from (2.15). It completes the proof.Remark 2.2.
In Theorem2.1, we assume that φ′(x)=1+(1/G)(∂Bz(x,y)/∂x)>0. The magnetic flux density Bz is smooth and its intensity is 10-7~10-8 T in practical experimental environments. Furthermore, the usual range of the reading gradient G is 10-3~10-4 T/m. Thus, the assumption of φ′(x)>0 is not severe in Theorem 2.1.


## 2.1. ConventionalBz Acquisition

For a conventionalBz acquisition, current is not injected during Ts of the MR data acquisition, ADC as shown in Figure 1. In this case, the induced magnetic flux density Bz provides additional dephasing of spins, and, consequently, extra phase is accumulated during the total injection time Tc. Then, the measured k-space data for the injection current I can be represented as follows:
(2.3)SI(kx,ky)=∫ℝ2‍ρ(x,y)eiδ(x,y)e-iγTcBz(x,y)e-i2π(kxx+kyy)dxdy,
where γ=26.75×107 rad/T·s is the gyromagnetic ratio of hydrogen.Figure 1
Conventional and ICNE current injections in a spin echo pulse sequence.From the measuredS and SI in (2.1) and (2.3), by applying an inverse Fourier transform, we obtain ρ0 in (2.2) and
(2.4)ρI(x,y)=ρ(x,y)eiδ(x,y)e-iγTcBz(x,y).
Then, the magnetic flux density Bz is precisely computed as
(2.5)Bz(x,y)=-1γTcarctan(α(x,y)β(x,y)),
where α and β are the imaginary and real parts of ρI/ρ0, respectively.


## 2.2. ICNEBz Acquisition

In the ICNEBz acquisition, in order to improve the SNR of Bz, we prolong the current injection time Tc until the end of the MR data acquisition, as shown in Figure 1. Then, since the induced Bz disturbs the linearity of the reading gradient, the measured k-space data has lost the linear encoding characteristic as
(2.6)SC(kx,ky)=∫ℝ2‍ρ(x,y)eiδ(x,y)e-iγTcBz(x,y)e-i2π[kx(x+Bz(x,y)/G)+kyy]dxdy,
where G is a constant that denotes the strength of the magnetic reading gradient. The inverse problem arising in the ICNE method is to recover Bz(x,y) from ρ0(x,y) obtained in (2.2) and the measured signal SC(kx,ky) in (2.6).Although the inversion is not uniquely solvable forBz(x,y) in general, we can uniquely determine Bz(x,y) by assuming that φ(x):=x+Bz(x,y)/G is monotone increasing in the following theorem.Theorem 2.1.
Letρ(x,y) have a finite support Ω. If Bz(x,y) is sufficiently small to guarantee that φ(x)=x+Bz(x,y)/G is monotone increasing so that φ′(x)>0 for each y, then Bz(x,y) is uniquely recovered in Ω from ρ0(x,y) in (2.2) and SC(kx,ky) in (2.6).Proof.
We note that the linear encoding characteristic in theky-variable remains unperturbed in (2.6). Thus, by one-dimensional inverse Fourier transform, SC in (2.6) is reduced to S^C in the (kx,y)-hybrid space as the following:
(2.7)S^C(kx,y)=∫-∞∞‍ρ(x,y)eiδ(x,y)e-iγTcBz(x,y)e-i2π[kx(x+Bz(x,y)/G)]dx.
Then, the ICNE inverse problem suffices to consider the x-directional inversion of Bz(x,y) from ρ0(x,y) in (2.2) and S^C(kx,y) in (2.7) for each fixed y.
By change of variables withφ(x), (2.7) is changed into
(2.8)S^C(kx,y)=∫-∞∞‍ρ(x,y)φ′(x)eiδ(x,y)e-iγTcBz(x,y)e-i2πkxφdφ.
From (2.8), by inverse Fourier transform for the φ-variable, φ satisfies
(2.9)Φ(φ)=ρ(x,y)φ′(x)eiδ(x,y)e-iγTcBz(x,y),
where Φ is a function defined with the given S^C by
(2.10)Φ(z)=∫-∞∞‍S^C(kx,y)ei2πkxzdkx.
The relation (2.9) gives us the simple ordinary differential equation as follows:
(2.11)|Φ(φ(x))|φ′(x)=ρ(x,y).
Since ρ has a finite support, for each y, we can define
(2.12)x0=min{x:ρ(x,y)≠0},β=min{z:Φ(z)≠0}.
If φ(x0)<β, then φ(x0+ϵ)<β for a sufficiently small ϵ>0. It contradicts (2.11), since it implies that
(2.13)0=|Φ(φ(x0+ϵ))|φ′(x0+ϵ)=ρ(x0+ϵ,y)≠0.
By the same argument, the reverse inequality is not possible. Thus, we have
(2.14)ϕ(x0)=β.
By separation of variables, (2.11) and (2.14) lead us to
(2.15)∫βφ(x)‍|Φ(φ)|dφ=∫x0x‍ρ(x,y)dx.
For any given x, φ(x) is uniquely determined from (2.15). It completes the proof.Remark 2.2.
In Theorem2.1, we assume that φ′(x)=1+(1/G)(∂Bz(x,y)/∂x)>0. The magnetic flux density Bz is smooth and its intensity is 10-7~10-8 T in practical experimental environments. Furthermore, the usual range of the reading gradient G is 10-3~10-4 T/m. Thus, the assumption of φ′(x)>0 is not severe in Theorem 2.1.


## 3. Discrete ICNE Inverse Problem

In a practical MRI scanner, the MRk-space data in (2.1), (2.3), and (2.6) are acquired by finite sampling with a dwell time dt. If N is the reading time Ts  divided by dt, we have the following N×N discrete signals with dimensionless variables instead of those in (2.6):
(3.1)SdC(k,l)=∑m=0N-1‍∑n=0N-1‍ρ(n,m)eiδ(n,m)e-i(2π/N)MBz(n,m)e-i(2π/N)(k(n+Bz(n,m))+lm),
for a constant M>0.The discrete ICNE inverse problem is to recoverBz(n,m) from (3.1) with known a prioriρ(n,m)eiδ(n,m) and the measured signal SdC(k,l), where n,m,k,l=0,1,…,N-1. By discrete inverse Fourier transform for l, (3.1) can be suppressed into
(3.2)S^dC(k,m)=∑n=0N-1‍ρ(n,m)eiδ(n,m)e-i(2π/N)(M+k)Bz(n,m)e-i(2π/N)kn.
For each fixed m, let s(k),ρn, and bn denote S^dC(k,m),ρ(n,m)eiδ(n,m), and Bz(n,m), respectively. Then, the discrete ICNE inversion problem is a system of N nonlinear equations for N unknowns, b0,b1,…,bN-1 such that
(3.3)s(k)=∑n=0N-1‍ρne-i(2π/N)(M+k)bne-i(2π/N)kn,k=0,1,…,N-1,
where M>0 and s(k),ρn∈ℂ,n,k=0,1,…,N-1 are known.In the rest of the paper, we assume thatN is even and k,n=0,1,…,N-1 denote the row and column numbers, respectively. A matrix whose (k,n) entry is Mkn is represented by
(3.4)[Mkn].
For a vector x=(x0,x1,…,xN-1)t, V(x)=[e-i(2π/N)(n+xn)k] denotes a Vandermonde matrix as
(3.5)V(x)=(11⋯1e-i(2π/N)x0e-i(2π/N)(1+x1)⋯e-i(2π/N)(N-1+xN-1)⋮⋮⋮⋮e-i(2π/N)x0(N-1)e-i(2π/N)(1+x1)(N-1)⋯e-i(2π/N)(N-1+xN-1)(N-1)).


### 3.1. Newton Iterations

Define a functionF=(F0,F1,…,FN-1)t by
(3.6)Fk(x)=∑n=0N-1‍ρne-i(2π/N)(M+k)xne-i(2π/N)kn-s(k),k=0,1,…,N-1,
for x=(x0,x1,…,xN-1)t∈ℝN. The discrete ICNE inverse problem is to find the zero of F for s(k) given in (3.3).The Jacobian matrixDF(x) is composed of four parts as
(3.7)DF(x)=-i2πNDKV(x)DxDρ,
where DK,Dxj,Dρ are diagonal matrices such that
(3.8)DK=diag(M,M+1,…,M+N-1),Dx=diag(e-i(2π/N)Mx0,e-i(2π/N)Mx1,…,e-i(2π/N)MxN-1),Dρ=diag(ρ0,ρ1,…,ρN-1).Newton iterations to find the zero ofF  are as the following:
(3.9)xj+1=xj-real(DF(xj)-1F(xj)),
with an initial x0 and the iterates xj=(x0j,x1j,…,xN-1j)t∈ℝN,j=0,1,2,….The previous method in [14, 15] was based on the Taylor approximation, but as a coincidental result, it can be interpreted as the first Newton iterate x1 in (3.9) with x0=0.


### 3.2. Convergence of Newton Iterations

Letb=(b0,b1,…,bN-1)t,ρ→=(ρ0,ρ1,…,ρN-1)t, and ρmax=maxn|ρn|,ρmin=minn|ρn|. If
(3.10)∥b∥∞<12,ρmin≠0,
the Jacobian DF(b) in (3.7) is invertible, since the Vandermonde matrix V(b) in (3.5) is based on the N distinct points. Thus, the Newton iterations in (3.9) converge to b for an initial x0 which is sufficiently close to b [17].The following theorem suggests a condition forb in which the Newton iterations in (3.9) converge to b with the trivial initial guess 0=(0,0,…,0)t. The proof is based on the Theorem 6.14 in [17], which states a sufficient condition for the convergence of Newton iterations that
(3.11)q:=αβγ<12,
where α=∥DF(0)-1F(0)∥∞ and β,γ are the respective bounds of
(3.12)∥DF(x)-1∥∞,∥DF(x)-DF(y)∥∞∥x-y∥∞for∥x∥∞,∥y∥∞<2α.Theorem 3.1.
Lets(k) in (3.6) be made through (3.3) from b such that
(3.13)∥b∥∞≤min{ρmin2πρmax(M+N),112π(2+logcotπ/2N),ρminM12πρmax(M+N)2}.
Then, starting with x0=0, the Newton iterations in (3.9) are well defined and converge to b. One also has the following quadratic error estimate:
(3.14)∥xj-b∥∞≤3∥b∥∞(12)2j-1,j=0,1,2,….
As a consequence, the zero of F satisfying (3.13) is unique.Proof.
Letα=∥DF(0)-1F(0)∥∞. The condition (3.13) and Lemma 3.3 lead us to
(3.15)12∥b∥∞≤α≤32∥b∥∞.
If ∥x∥∞≤2α, we have from (3.13), (3.15), and Lemma 3.5, the following:
(3.16)∥V(x)-1∥∞≤2.
By (3.7), (3.16) implies that
(3.17)∥DF(x)-1∥∞≤NπρminM,if∥x∥∞<2α.
For two constants in (3.17) and (3.21), let
(3.18)q=αNπρminM4π2N(N+M)2ρmax.
From (3.15) and the condition (3.13), we have
(3.19)q≤12.
Then, with the aid of the Theorem 6.14 in [17], F has a unique zero x⋆ in the ball
(3.20)B={x:∥x∥∞≤2α},
and the Newton iterations in (3.9) converge to x⋆ with x0=0. Since b is a zero of F contained in B from (3.15), we have x⋆=b. The quadratic error estimate in (3.14) also comes from the same theorem in [17].Lemma 3.2.
Ifx,y∈ℝN, then
(3.21)∥DF(x)-DF(y)∥∞≤4π2N(N+M)2ρmax∥x-y∥∞.Proof.
From (3.7), we expand
(3.22)DF(x)-DF(y)=-i2πNDK[e-i(2π/N)((n+xn)k+Mxn)-e-i(2π/N)((n+yn)k+Myn)]Dρ=-i2πNDK[(e-i(2π/N)((xn-yn)(k+M))-1)e-i(2π/N)((n+yn)k+Myn)]Dρ.
Since|eiθ-1|=2|sin(θ/2)|, we have for each k,(3.23a)∑n=0N-1|e-i(2π/N)((xn-yn)(k+M))-1|=2∑n=0N-1|sinπN(xn-yn)(k+M)|(3.23b)≤2∑n=0N-1πN|xn-yn|(k+M)(3.23c)≤2π(k+M)∥x-y∥∞.We can combine (3.22), (3.23a), (3.23b), and (3.23c) into (3.21).Lemma 3.3.
For the trivial initialx0=0, one has
(3.24)∥b∥∞-ρmaxρminπ(M+N)∥b∥∞2≤∥DF(0)-1F(0)∥∞≤∥b∥∞+ρmaxρminπ(M+N)∥b∥∞2.Proof.
Since(3.25)Fk(0)=∑n=0N-1‍ρne-i(2π/N)kn-∑n=0N-1‍ρne-i(2π/N)(M+k)bne-i(2π/N)kn,
we can represent F(0) as
(3.26)F(0)=[(1-e-i(2π/N)(M+k)bn)e-i(2π/N)kn]ρ→.
Thus, we expand(3.27a)DF(0)-1F(0)=(-i2πN)-1Dρ-1V(0)-1DK-1[(1-e-i(2π/N)(M+k)bn)e-i(2π/N)kn]ρ→(3.27b)=Dρ-1V(0)-1[1-e-i(2π/N)(M+k)bn-i(2π/N)(M+k)e-i(2π/N)kn]ρ→(3.27c)=Dρ-1V(0)-1[1-e-i(2π/N)(M+k)bn-i(2π/N)(M+k)bne-i(2π/N)kn]Dρb(3.27d)=Dρ-1V(0)-1(E-V(0))Dρb(3.27e)=Dρ-1V(0)-1EDρb-b,where
(3.28)E=[1-e-i(2π/N)(M+k)bn-i(2π/N)(M+k)bne-i(2π/N)kn+e-i(2π/N)kn].
From Lemma3.4, we estimate that
(3.29)∥E∥∞=maxk∑n=1N-1‍|Ekn|≤maxk∑n=1N-1‍πN(M+k)|bn|≤π(M+N)∥b∥∞.
Since V(0) is the matrix of the discrete Fourier transform, we have
(3.30)V(0)-1=1NV(0)′,
which implies that
(3.31)∥V(0)∥∞=N,∥V(0)-1∥∞=1.
The proof is completed by (3.27a), (3.27b), (3.27c), (3.27d), (3.27e), (3.29), and (3.31).Lemma 3.4.
Ifθ is real, then
(3.32)|1-eiθiθ+1|≤12|θ|.Proof.
If|θ|≤3, we have (3.32) since the series in the following expansion is alternating. 									(3.33)|eiθ-1-iθ|2=(cosθ-1)2+(sinθ-θ)2=2(1-cosθ-θsinθ)+θ2=2[(θ22-θ44!+θ66!-⋯)-θ(θ-θ33!+θ55!-⋯)]+θ2=2[(13!-14!)θ4-(15!-16!)θ6+(17!-18!)θ8-⋯]=2(34!θ4-56!θ6+78!θ8+⋯)≤14θ4.
If |θ|>3, we obtain (3.32) from
(3.34)|eiθ-1-iθ|2≤2(2+|θ|)+θ2≤14θ4.


### 3.3. Norm of Inverse of Vandermonde Matrix

The norms of inverses of Vandermonde matrices were estimated by Gautschi [18, 19]. In some estimations there, the equality holds if all base points are on the same ray through the origin.Compared to (3.31), for a small perturbation x, a bound of ∥V(x)-1∥∞ is investigated in the following lemma. Since the norm estimation of inverse of Vandermonde matrix must be interesting, we separate the result in this subsection from other ingredients for Theorem 3.1.Lemma 3.5.
If∥x∥∞<1/4π(2+logcot(π/2N)), then
(3.35)∥V(x)-1∥∞≤2.Proof.
LetI=[δkn] be the identity matrix and W=V(0)-1V(x)-I, whose entries are
(3.36)Wkn=1N1-e-i2πxn1-ei(2π/N)(k-n-xn)-δkn.
Since eiθ-1=2isin(θ/2)eiθ/2, the off-diagonal entries satisfy the following
(3.37)|Wkn|=1N|sinπxnsin(π/N)(k-n-xn)|≤1N|πxnsin(π/N)(k-n-xn)|.
The diagonal entries are estimated by the Cauchy mean value theorem as follows:
(3.38)|Wkk|=|1N1-e-i2πbk1-e-i(2π/N)xk-1|=|1Nsinπxksin(π/N)xke-i((N-1)/N)πxk-1|=|cosπxk′cos(π/N)xk′e-i((N-1)/N)πxk-1|,forsome|xk′|≤|xk|,≤|cosπxk′cos(π/N)xk′e-i((N-1)/N)πxk-e-i((N-1)/N)πxk|+|e-i((N-1)/N)πxk-1|≤|cosπxk′cos(π/N)xk′-1|+π|xk|≤3π|xk|.
Regarding summations of|Wkn|, the maximum occurs when k=N/2±1 from the symmetry of the sine function in (3.37). In both cases, we have(3.39a)∑n=0N-1‍|Wkn|≤3π|xk|+2∑m=1N/2‍1N|π∥x∥∞sin(π/N)(m-∥x∥∞)|(3.39b)≤3π|xk|+2Nπ∥x∥∞sin(π/N)(1-∥x∥∞)+2∥x∥∞∫π/Nπ/2‍1sintdt(3.39c)≤(4π+2πlogcotπ2N)∥x∥∞≤12.Since (3.39a), (3.39b), and (3.39c) imply ∥W∥∞≤1/2, we establish the following:
(3.40)∥V(x)-1∥∞≤∥V(x)-1V(0)∥∞∥V(0)-1∥∞=∥V(x)-1V(0)∥∞=∥(I+W)-1∥∞≤11-∥W∥∞≤2.


## 3.1. Newton Iterations

Define a functionF=(F0,F1,…,FN-1)t by
(3.6)Fk(x)=∑n=0N-1‍ρne-i(2π/N)(M+k)xne-i(2π/N)kn-s(k),k=0,1,…,N-1,
for x=(x0,x1,…,xN-1)t∈ℝN. The discrete ICNE inverse problem is to find the zero of F for s(k) given in (3.3).The Jacobian matrixDF(x) is composed of four parts as
(3.7)DF(x)=-i2πNDKV(x)DxDρ,
where DK,Dxj,Dρ are diagonal matrices such that
(3.8)DK=diag(M,M+1,…,M+N-1),Dx=diag(e-i(2π/N)Mx0,e-i(2π/N)Mx1,…,e-i(2π/N)MxN-1),Dρ=diag(ρ0,ρ1,…,ρN-1).Newton iterations to find the zero ofF  are as the following:
(3.9)xj+1=xj-real(DF(xj)-1F(xj)),
with an initial x0 and the iterates xj=(x0j,x1j,…,xN-1j)t∈ℝN,j=0,1,2,….The previous method in [14, 15] was based on the Taylor approximation, but as a coincidental result, it can be interpreted as the first Newton iterate x1 in (3.9) with x0=0.


## 3.2. Convergence of Newton Iterations

Letb=(b0,b1,…,bN-1)t,ρ→=(ρ0,ρ1,…,ρN-1)t, and ρmax=maxn|ρn|,ρmin=minn|ρn|. If
(3.10)∥b∥∞<12,ρmin≠0,
the Jacobian DF(b) in (3.7) is invertible, since the Vandermonde matrix V(b) in (3.5) is based on the N distinct points. Thus, the Newton iterations in (3.9) converge to b for an initial x0 which is sufficiently close to b [17].The following theorem suggests a condition forb in which the Newton iterations in (3.9) converge to b with the trivial initial guess 0=(0,0,…,0)t. The proof is based on the Theorem 6.14 in [17], which states a sufficient condition for the convergence of Newton iterations that
(3.11)q:=αβγ<12,
where α=∥DF(0)-1F(0)∥∞ and β,γ are the respective bounds of
(3.12)∥DF(x)-1∥∞,∥DF(x)-DF(y)∥∞∥x-y∥∞for∥x∥∞,∥y∥∞<2α.Theorem 3.1.
Lets(k) in (3.6) be made through (3.3) from b such that
(3.13)∥b∥∞≤min{ρmin2πρmax(M+N),112π(2+logcotπ/2N),ρminM12πρmax(M+N)2}.
Then, starting with x0=0, the Newton iterations in (3.9) are well defined and converge to b. One also has the following quadratic error estimate:
(3.14)∥xj-b∥∞≤3∥b∥∞(12)2j-1,j=0,1,2,….
As a consequence, the zero of F satisfying (3.13) is unique.Proof.
Letα=∥DF(0)-1F(0)∥∞. The condition (3.13) and Lemma 3.3 lead us to
(3.15)12∥b∥∞≤α≤32∥b∥∞.
If ∥x∥∞≤2α, we have from (3.13), (3.15), and Lemma 3.5, the following:
(3.16)∥V(x)-1∥∞≤2.
By (3.7), (3.16) implies that
(3.17)∥DF(x)-1∥∞≤NπρminM,if∥x∥∞<2α.
For two constants in (3.17) and (3.21), let
(3.18)q=αNπρminM4π2N(N+M)2ρmax.
From (3.15) and the condition (3.13), we have
(3.19)q≤12.
Then, with the aid of the Theorem 6.14 in [17], F has a unique zero x⋆ in the ball
(3.20)B={x:∥x∥∞≤2α},
and the Newton iterations in (3.9) converge to x⋆ with x0=0. Since b is a zero of F contained in B from (3.15), we have x⋆=b. The quadratic error estimate in (3.14) also comes from the same theorem in [17].Lemma 3.2.
Ifx,y∈ℝN, then
(3.21)∥DF(x)-DF(y)∥∞≤4π2N(N+M)2ρmax∥x-y∥∞.Proof.
From (3.7), we expand
(3.22)DF(x)-DF(y)=-i2πNDK[e-i(2π/N)((n+xn)k+Mxn)-e-i(2π/N)((n+yn)k+Myn)]Dρ=-i2πNDK[(e-i(2π/N)((xn-yn)(k+M))-1)e-i(2π/N)((n+yn)k+Myn)]Dρ.
Since|eiθ-1|=2|sin(θ/2)|, we have for each k,(3.23a)∑n=0N-1|e-i(2π/N)((xn-yn)(k+M))-1|=2∑n=0N-1|sinπN(xn-yn)(k+M)|(3.23b)≤2∑n=0N-1πN|xn-yn|(k+M)(3.23c)≤2π(k+M)∥x-y∥∞.We can combine (3.22), (3.23a), (3.23b), and (3.23c) into (3.21).Lemma 3.3.
For the trivial initialx0=0, one has
(3.24)∥b∥∞-ρmaxρminπ(M+N)∥b∥∞2≤∥DF(0)-1F(0)∥∞≤∥b∥∞+ρmaxρminπ(M+N)∥b∥∞2.Proof.
Since(3.25)Fk(0)=∑n=0N-1‍ρne-i(2π/N)kn-∑n=0N-1‍ρne-i(2π/N)(M+k)bne-i(2π/N)kn,
we can represent F(0) as
(3.26)F(0)=[(1-e-i(2π/N)(M+k)bn)e-i(2π/N)kn]ρ→.
Thus, we expand(3.27a)DF(0)-1F(0)=(-i2πN)-1Dρ-1V(0)-1DK-1[(1-e-i(2π/N)(M+k)bn)e-i(2π/N)kn]ρ→(3.27b)=Dρ-1V(0)-1[1-e-i(2π/N)(M+k)bn-i(2π/N)(M+k)e-i(2π/N)kn]ρ→(3.27c)=Dρ-1V(0)-1[1-e-i(2π/N)(M+k)bn-i(2π/N)(M+k)bne-i(2π/N)kn]Dρb(3.27d)=Dρ-1V(0)-1(E-V(0))Dρb(3.27e)=Dρ-1V(0)-1EDρb-b,where
(3.28)E=[1-e-i(2π/N)(M+k)bn-i(2π/N)(M+k)bne-i(2π/N)kn+e-i(2π/N)kn].
From Lemma3.4, we estimate that
(3.29)∥E∥∞=maxk∑n=1N-1‍|Ekn|≤maxk∑n=1N-1‍πN(M+k)|bn|≤π(M+N)∥b∥∞.
Since V(0) is the matrix of the discrete Fourier transform, we have
(3.30)V(0)-1=1NV(0)′,
which implies that
(3.31)∥V(0)∥∞=N,∥V(0)-1∥∞=1.
The proof is completed by (3.27a), (3.27b), (3.27c), (3.27d), (3.27e), (3.29), and (3.31).Lemma 3.4.
Ifθ is real, then
(3.32)|1-eiθiθ+1|≤12|θ|.Proof.
If|θ|≤3, we have (3.32) since the series in the following expansion is alternating. 									(3.33)|eiθ-1-iθ|2=(cosθ-1)2+(sinθ-θ)2=2(1-cosθ-θsinθ)+θ2=2[(θ22-θ44!+θ66!-⋯)-θ(θ-θ33!+θ55!-⋯)]+θ2=2[(13!-14!)θ4-(15!-16!)θ6+(17!-18!)θ8-⋯]=2(34!θ4-56!θ6+78!θ8+⋯)≤14θ4.
If |θ|>3, we obtain (3.32) from
(3.34)|eiθ-1-iθ|2≤2(2+|θ|)+θ2≤14θ4.


## 3.3. Norm of Inverse of Vandermonde Matrix

The norms of inverses of Vandermonde matrices were estimated by Gautschi [18, 19]. In some estimations there, the equality holds if all base points are on the same ray through the origin.Compared to (3.31), for a small perturbation x, a bound of ∥V(x)-1∥∞ is investigated in the following lemma. Since the norm estimation of inverse of Vandermonde matrix must be interesting, we separate the result in this subsection from other ingredients for Theorem 3.1.Lemma 3.5.
If∥x∥∞<1/4π(2+logcot(π/2N)), then
(3.35)∥V(x)-1∥∞≤2.Proof.
LetI=[δkn] be the identity matrix and W=V(0)-1V(x)-I, whose entries are
(3.36)Wkn=1N1-e-i2πxn1-ei(2π/N)(k-n-xn)-δkn.
Since eiθ-1=2isin(θ/2)eiθ/2, the off-diagonal entries satisfy the following
(3.37)|Wkn|=1N|sinπxnsin(π/N)(k-n-xn)|≤1N|πxnsin(π/N)(k-n-xn)|.
The diagonal entries are estimated by the Cauchy mean value theorem as follows:
(3.38)|Wkk|=|1N1-e-i2πbk1-e-i(2π/N)xk-1|=|1Nsinπxksin(π/N)xke-i((N-1)/N)πxk-1|=|cosπxk′cos(π/N)xk′e-i((N-1)/N)πxk-1|,forsome|xk′|≤|xk|,≤|cosπxk′cos(π/N)xk′e-i((N-1)/N)πxk-e-i((N-1)/N)πxk|+|e-i((N-1)/N)πxk-1|≤|cosπxk′cos(π/N)xk′-1|+π|xk|≤3π|xk|.
Regarding summations of|Wkn|, the maximum occurs when k=N/2±1 from the symmetry of the sine function in (3.37). In both cases, we have(3.39a)∑n=0N-1‍|Wkn|≤3π|xk|+2∑m=1N/2‍1N|π∥x∥∞sin(π/N)(m-∥x∥∞)|(3.39b)≤3π|xk|+2Nπ∥x∥∞sin(π/N)(1-∥x∥∞)+2∥x∥∞∫π/Nπ/2‍1sintdt(3.39c)≤(4π+2πlogcotπ2N)∥x∥∞≤12.Since (3.39a), (3.39b), and (3.39c) imply ∥W∥∞≤1/2, we establish the following:
(3.40)∥V(x)-1∥∞≤∥V(x)-1V(0)∥∞∥V(0)-1∥∞=∥V(x)-1V(0)∥∞=∥(I+W)-1∥∞≤11-∥W∥∞≤2.


## 4. Numerical Results

From the Biot-Savart law, we simulate the magnetic flux densityBz(n,m) induced by a horizontal current through the Logan shape ρ as in Figure 2, where N=60. We obtain the simulated MR signals SdC(k,l) as depicted in Figure 3, through (3.1) with M=32,δ=0. We note that the maximum of Bz(n,m) is about 1/4 in Figure 2, larger than suggested in Theorem 3.1.SimulatedBz and ρ.
(a)
SimulatedBz by transversal current
(b)
Logan shapeρFigure 3
Simulated dataSdC from Bz and ρ.By discrete inverse Fourier transform forl, we transform SdC(k,l) into S^dC(k,m) in (3.2). Then, setting for each fixed m,
(4.1)s(k)=S^dC(k,m),
the Newton iterations in (3.9) generate xj with x0=0. The jth approximation Bzj is done by
(4.2)Bz(n,m)j=xj(n).Thelog of error maxn,m|Bz(n,m)-Bz(n,m)j| is given in Figure 4(a), which means that the error decay is quadratic. In the Newton iterations in (3.9), we have to solve a Vandermonde system for V(x), which may be consuming time or unstable. Instead of V(x) in (3.7), we can fix V(0) and simplify the Newton iterations in (3.9). Then, the error decay is reduced to be linear as in Figure 4(b).log
max
n
,
m
|
B
z
(
n
,
m
)
-
B
z
(
n
,
m
)
j
| error decay.
(a)
Newton
(b)
Simplified Newton



---

*Source: 973781-2012-12-16.xml*