Recursive Algorithms for Computing the Cramer-Rao Bound^{1}

Mohammad Usman, Alfred O. Hero, Anne Sauve, and Jeffrey A. Fessler

Abstract

Computation of the Cramer-Rao bound requires the inverse or the pseudo-inverse Fisher information matrix (FIM). Direct matrix inversion can be computationally intractable when the number of unknown parameters is large. In this paper we compare several recursive methods for extracting columns of the inverse Fisher information matrix including the algorithm of [1] and similar monotone matrix splitting algorithms, over-relaxation techniques, and preconditioned conjugate gradient algorithms. The latter algorithms can converge significantly faster than the monotone algorithms for equivalent computational load. We illustrate the methods developed in this paper for an important class of inverse problems.

Key Words: Performance bounds, multi-dimensional parameter estimation, monotone matrix splitting iterations, preconditioned conjugate gradient.

EDICS Number: SP 3.8.1

1 Introduction

The Cramer-Rao (CR) bound is a very widely used lower bound on the covariance of any unbiased estimator. For an n-dimensional parameter vector calculation of the CR bound involves finding elements of the inverse of the n × n Fisher information matrix (FIM) F_{Y}. Direct matrix inversion, requiring O(n^{2}) bytes of memory storage and O(n^{3}) flops (floating point operations), is intractable if n is large. For example, in image reconstruction where the pixel intensities are the parameters to be estimated, a moderate size image of 128 × 128 = 16384 pixels has a FIM of dimension 16384 × 16384 which would require 268 Megabytes memory storage (single precision) and 4 Tera flops (O(10^{−12}) flops) for computation of the inverse FIM.

Often only a few components of the n-dimensional estimator are of interest in which case the entire inverse FIM is not needed. For example in image reconstruction one may be primarily interested in a small q-pixel region of interest (ROI) corresponding to a tumor or lesion. Methods of sequential partitioning [2] can be used in this case to calculate the q × q submatrix of the inverse Fisher matrix yielding the CR bound for the ROI. However, while requiring fewer flops than required for full matrix inversion, this algorithm still requires the same order O(n^{3}) of flops.

In [1] and [3] a recursive method was presented for approximating columns of the CR bound. This method was derived by specifying a FIM-dominating matrix which allowed a geometric series decomposition of the inverse FIM. The geometric series algorithm requires only O(n^{2}) flops per iteration per parameter so that if convergence speed is fast, a significant computational saving is achieved. The main advantage of this algorithm is its monotone convergence which guarantees a valid and improving lower bound on estimator covariance at each iteration. The disadvantage is slow speed of convergence. In [3, 1] special properties of the Fisher information matrix were exploited to specify the FIM-dominating matrix. In Section 3.2 we give an alternative dominating matrix based on purely algebraic considerations.

This paper places the methods of [1] in the context of the general class of iterative linear equation solvers known as stationary and non-stationary algorithms [4]. These methods include monotone and non-monotone matrix splitting algorithms, such as the method of [1], Jacobi, Gauss-Seidel, and preconditioned conjugate gradient algorithms. We also present an algorithm for computing the pseudo-inverse version of the CR bound for singular FIM based on matrix perturbation and differencing methods. These iterative algorithms are illustrated for a class of inverse problems arising in deconvolution, image restoration, and tomographic imaging.

The main conclusions of this paper are: 1) the nonmonotone Gauss-Seidel and pre-conditioned conjugate gradient methods have an order of magnitude speed advantage relative to the monotone methods; 2) for sparse matrix inverse problems these methods can be implemented with significant savings in memory and computation; 3) for inverse problems having singular FIM the application of Gauss-Seidel and conjugate gradient algorithms yields stable and rapidly convergent methods for finding columns of the pseudoinverse of the FIM.

1J. Fessler, A. Hero, and A. Sauve are with the Dept. of Electrical Engineering and Computer Science, The University of Michigan, Ann Arbor, MI 48109-2122. M. Usman is with the Punjab Institute of Computer Science, Pakistan. This work was supported in part by the National Science Foundation under grant BCS-9024370, a Government of Pakistan Postgraduate Fellowship, NIH grants CA-54362 and CA-60711, and DOE grant DE-FG02-87ER60561

2 The Cramer-Rao Bound

Let Y be an observed random variable with probability density f_{Y} (y; \underline{θ}) dependent on an unknown parameter vector \underline{θ} = [θ_{1}, ..., θ_{n}]^{T} lying in an open subset Θ of IR^{n}. The covariance cov\underline{θ}(ˆ\underline{θ}) = E\underline{θ} [(ˆ\underline{θ} − E\underline{θ}(ˆ\underline{θ}))(ˆ\underline{θ} − E\underline{θ}(ˆ\underline{θ}))^{T} ] of any

unbiased estimator ˆ \underline{θ} = ˆ \underline{θ}(Y) of \underline{θ} obeys the Cramer-Rao lower bound [5, 6]

cov\underline{θ}(ˆ\underline{θ}) ≥ F_{Y}^{+}(\underline{θ}),  (1)

where F_{Y} is the FIM

F_{Y} (\underline{θ}) = E\underline{θ}[∇\underline{θ} ln f_{Y} (Y; \underline{θ}) ∇T ; \underline{θ} ln f_{Y} (Y; \underline{θ})],

the column ∇\underline{θ} denotes the gradient vector [\frac{∂}{∂θ_{1}} , \frac{∂}{∂θ_{2}} , ..., \frac{∂}{∂θ_{n}}]^{T}, and F_{Y}^{+} denotes the Moore-Penrose pseudo-inverse [7]. When F_{Y} is non-singular F_{Y}^{+} = F_{Y}^{−1}. Note that the pseudo-inverse form of the CR bound is generally only achievable for estimation of components of \underline{θ} in the rangespace of F_{Y} [?].

The matrix inequality (1) means that cov\underline{θ}(ˆ\underline{θ}) − F_{Y}^{+}(\underline{θ}) is non-negative definite. In particular the lower bound on the variance of θ̂_{1}, corresponding to the top left (1,1) entry of cov\underline{θ}(ˆ\underline{θ}), is given by the corresponding entry of F_{Y}^{+}

var\underline{θ}(θ̂_{1}) ≥ \underline{e}_{1}^{T} F_{Y}^{+}(\underline{θ})\underline{e}_{1},  (2)

where \underline{e}_{1} = [1, 0, ..., 0]^{T}. Throughout this paper we will be interested in calculating the right hand side of (2), or more generally the first column of F_{Y}. The results are easily extended to q -elements, q << n, by repeated application of the algorithms given here.

Direct computation of the right hand side of the CR inequality (2) can be accomplished in O(n^{3}) flops by solving for \underline{x} in the equation

F_{Y} \underline{x} = \underline{e}_{1}.  (3)

For example, when F_{Y} is non-singular the solution \underline{x} can be found by the method of Gaussian elimination (2n^{3}/3 flops) or by methods which exploit positive definiteness of F_{Y}, such as methods based on the QR and Cholesky decomposition (n^{3}/3 flops) [8, Sec. 4.2].

3 Recursive CR Bound Algorithms for Non-Singular FIM

Here we describe several methods for iteratively approximating the right hand side of the CR inequality (2) for non-singular FIM. We first describe the monotonically convergent algorithm of [1] in the context of standard splitting iterations [8, Sec. 10.1] , also known as stationary iterations [4], for approximating the solution to the linear equation (3).

3.1 Splitting Algorithms

Let F and N be n×n matrices which split F_{Y} in the sense that F − N = F_{Y} and assume that F is non-singular.

3.1.1 General Splitting Iterations

Define the general splitting iterations for approximation of the CR bound

INITIALIZE: η^{(0)} = 0, \underline{β}^{(0)} = \underline{0} an n × 1 vector of zeroes.

RECURSION:

\begin{matrix}i). & \underline{t} = N\underline{β}^{(k)} + \underline{e}_{1} \\ ii). & Solve : F\underline{β}^{(k+1)} = \underline{t} \end{matrix}  (4)

CR BOUND APPROXIMATION:

η^{(k+1)} = \underline{e}_{1}^{T} \underline{β}^{(k+1)}.

Define the matrix M = F^{−1}N. The root convergence factor ρ(M) = |λ_{i}^{M}|_{max} is defined as the maximum magnitude eigenvalue of M, also known as the spectral radius of M. If the ρ(F^{−1}N) < 1 then \underline{β}^{(k)} converges to the first column of F_{Y}^{−1} [9, 8] and the approximating sequence η^{(k)} converges to the CR bound \underline{e}_{1}^{T} F_{Y}^{−1}\underline{e}_{1} . Furthermore, it can be shown in an analogous manner to [1, Eq. (9)] that if F is symmetric positive definite and N is symmetric non-negative definite then η^{(k+1)} − η^{(k)} = \underline{e}_{1}^{T} F^{− \frac{1}{2}} [F^{− \frac{1}{2}} NF^{− \frac{1}{2}} ]^{k}F^{− \frac{1}{2}} \underline{e}_{1}. Hence, under these conditions the approximating sequence η^{(k)} is monotone nondecreasing in k and, if ρ(F^{−1}N) < 1, η^{(k)} converges to the CR bound from below.

Many well known algorithms fall into the category of general splitting iterations such as the Jacobi and Gauss-Siedel iterations. Let the FIM have the additive decomposition F_{Y} = D + U + L where D is diagonal and U and L are upper and lower diagonal matrices with zero diagonal entries. The Jacobi (J) iteration is obtained by making the identifications F = D, N = −(U+L) in (4). For Jacobi iterations, the spectral radius of F^{−1}N may exceed one and N is not generally non-negative definite. Therefore, the Jacobi iterations may not converge and are generally not monotone. The Gauss-Siedel (GS) iteration is obtained by identifying F = D + L and N = −U in (4). Like Jacobi the Gauss-Seidel iterations yield non-monotonic approximations. On the other hand, the GS iterations always converge for positive definite FIM.

For large n step (i) in (4) requires 2n^{2} flops. If F is a diagonal matrix, as in Jacobi, step (ii) requires n flops. For F a banded diagonal matrix with upper and lower bandwidth of p−1, (p = 1 gives a diagonal matrix) step (ii) can be computed with 2np^{2} flops (n ≫ p) using the method of banded Gaussian elimination [8, Sec. 4.3]. Finally, when F is lower triangular, as in Gauss-Siedel, step (ii) can be implemented by the method of back-substitution requiring n^{2} flops.

3.1.2 Monotone Splitting (MS) Iterations

Assume that F_{Y} is positive definite and that one can find a positive definite n × n matrix F such that F dominates F_{Y} in the sense that F − F_{Y} is non-negative definite. Now identify N = F−F_{Y} and note that F and N give a matrix splitting of F_{Y}. Furthermore N is non-negative definite and all the eigenvalues of F^{−1}N = I − F^{−1}F_{Y} lie in the interval [0, 1) [1]. Hence, using this choice of F and N in the general splitting iteration (4) we obtain an algorithm which gives an approximating sequence η^{(k)} which converges monotonically to the CR bound from below. This monotonic splitting (MS) algorithm yields a sequence of increasingly tight bounds on var\underline{θ}(θ̂_{1}).

In [1] a diagonal splitting matrix, denoted F_{EM}, was specified as the FIM associated with the complete data for a complete/incomplete data formulation. Since complete data is always at least as informatative as incomplete data F_{EM} ≥ F_{Y} and thus use of F_{EM} as the splitting matrix F in (4) yields a monotone MS algorithm.

If one is willing to sacrifice monotonic convergence, any monotonic MS algorithm can be accelerated by introducing a stationary (time invariant) over-relaxation parameter ψ ∈ (1, 2) [10]. Specifically, assume that a dominating F has been found. An over-relaxed MS algorithm (MSOR) is specified by using N = F−ψF_{Y} in the recursion (4). Note that the accelerated algorithm is generally non-monotonic since N can lose the property of non-negative definiteness when ψ > 1.

3.1.3 Convergence of Splitting Iterations

For the algorithm (4) to be computationally advantageous relative to direct Gaussian elimination methods we must ensure that the per iteration computational cost be low and that the speed of convergence be fast. Define the approximation error Δβ^{(k)} = F_{Y}^{−1}\underline{e}_{1} − β^{(k)} at the k-th iteration of the recursion (4). As previously mentioned, it is well known that the the error will asymptotically decrease to zero as long as ρ(M) < 1 [9, 8] where M = F^{−1}N. In this section we give a useful non-asymptotic convergence criterion which is used to derive an upper bound on the number of iterations required for convergence to within a user-specified tolerance.

From the general splitting iteration (4) it is easily shown that

Δβ^{(k+1)} = M^{k}Δβ^{(0)},

where Δβ^{(0)} = F_{Y}^{−1}\underline{e}_{1}. Using properties of vector and matrix norms [8, Sec. 2.3]

∥Δβ^{(k+1)}∥_{2} ≤ |||M^{k}|||_{2} · ∥Δβ^{(0)}∥_{2}
≤ |||M|||_{2}^{k} · ∥Δβ^{(0)}∥_{2},  (5)

where ∥\underline{u}∥_{2} = \sqrt{\underline{u}^{T} \underline{u}} and |||M|||_{2} = max\underline{x} \frac{∥M\underline{x}∥_{2}}{∥\underline{x}∥_{2}} denote the euclidean norm and the matrix-2 norm, respectively. The matrix-2 norm of M is equal to the maximum magnitude singular value of M.

Note that if the initial error Δβ^{(0)} is equal to zero subsequent errors remains zero. For non-zero Δβ^{(0)} we specify a convergence criterion in terms of a user specified threshold γ ∈ [0, 1]. Definition An iterative algorithm is said to have converged after N iterations (to within a factor γ ∈ [0, 1]) if

\frac{∥Δβ^{(k)}∥_{2}}{∥Δβ^{(0)}∥_{2}} < γ, ∀k ≥ N.  (6)

The minimum integer N = N(γ) such that this inequality is satisfied is called the required number of iterations for convergence.

Stopping the recursions at k = N will ensure convergence to the CR bound within the factor γ. Now assume that the integer N satisfies |||M|||_{2}^{N−1} < γ < 1. Then from inequality (5) \frac{∥Δβ^{(k)}∥_{2}}{∥Δβ^{(0)}∥_{2}} < γ for all k ≥ N. Thus we have the following upper bound on N(γ)

N(γ) ≤ N^{u}(γ) = \frac{log γ}{log |||M|||_{2}} + 1.  (7)

As γ approaches zero the required number of iterations N(γ) goes to infinity and the root convergence factor can be directly related to the ratio ∥Δβ^{(k)}∥_{2}/∥Δβ^{(0)}∥_{2} through

(5) and the identity ρ(M) = lim_{N→∞} (|||M^{N}|||_{2})^{\frac{1}{N}} [8]. However, as pointed out in [10], while the root convergence factor is a good predictor of the asymptotic behavior of this ratio, the matrix-2 norm |||M|||_{2} is more reliable as a predictor of the actual number of iterations required for convergence. This will be quantitatively illustrated in Section 6.

If an accurate and simply computed estimate of |||M|||_{2} is available, e.g. from application of power iteration methods [8], then the function N^{u}(γ) can be computed. In this case the relation (7) can be used to determine if the standard splitting algorithm (4) is likely to require fewer flops than Gaussian elimination. To illustrate, recall that for F a p-diagonal banded matrix each iteration of the splitting algorithm (4) requires 2n^{2} + 2np^{2}. Hence, for specified γ ∈ [0, 1], the splitting algorithm (4) is guaranteed to require fewer flops than Gaussian elimination (n^{3}/3 flops or more) when

N^{u}(γ) < \frac{n^{2}}{n + p^{2}} \frac{3}{2}.  (8)

3.2 A Simple Class of Monotone Splitting Matrices

If monotone convergence of recursion (4) is to be ensured we need select dominating splitting matrix F such that N = F−F_{Y} is non-negative definite. Furthermore, to promote fast convergence, |||F^{−1}N|||_{2} should be small. This suggests posing the constrained optimization problem: find F to minimize |||I − F^{−1}F_{Y} |||_{2} subject to the constraint F−F_{Y} ≥ 0. It can be simply verified that even for the case of diagonal F the optimum F is impractical to implement as it is a function of the eigenvalues of F_{Y}.

As an alternative, consider the following p-diagonal splitting matrix F = D_{p} defined in terms of the elements f_{ij} of the matrix F_{Y}

D_{p} = Q + diag (|F_{Y} − Q| \underline{1}) .  (9)

where Q = ((f_{ij}))_{|i−j|<p} is a p-diagonal matrix, |A| denotes a matrix whose elements are the absolute values of those of the matrix A, \underline{1} = [1, …, 1]^{T}, and diag(\underline{x}) is a diagonal matrix with the elements of the vector \underline{x} along the diagonal. In particular, D_{1} is a diagonal matrix with i-th diagonal element ∑_{j=1}^{n} |f_{ij}| and D_{2} is the tridiagonal matrix

D_{2} = [ \begin{matrix}∑ |f_{1j}| ; j≠2 & f_{12} & & & ⋯ & 0 \\ f_{21} & ∑ |f_{2j}| ; j≠1,3 & ... & & & ⋮ \\ & ... & ... & ... \\ & & ... & ∑ |f_{kj}| ; j≠k−1,k+1 & ... \\ ⋮ & & & ... & ... & f_{n−1,n} \\ 0 & ⋯ & & & f_{n,n−1} & ∑ |f_{nj}| ; j≠n−1 \end{matrix} ] .

The following lemma follows directly from the diagonal dominance of the matrix D_{p} − F_{Y} = diag (|F_{Y} − Q| \underline{1}) − (F_{Y} −Q_{p}) [9, Corollary 7.2.1] and the easily verifiable fact that when F_{Y} has non-negative entries: (D_{p} − F_{Y})\underline{1} = 0.

Lemma 1 Assume that F_{Y} is an n×n symmetric matrix. Then D_{p}−F_{Y} is non-negative definite. Furthermore, if F_{Y} has only non-negative entries then D_{p} − F_{Y} has rank at most n − 1.

It can be shown that a necessary condition for |||I − F^{−1}F|||_{2} to be minimized by an F_{Y} -dominating p-diagonal matrix F = D_{p} is that D_{p} −F be rank deficient. Lemma 1 asserts that for non-negative Fisher matrices D_{p} satisfies this condition. Non-negative Fisher matrices arise in many applications including the Poisson inverse problem to be considered in Sec. 6. Extensive numerical experiments have indicated that F = D_{p} comes very close to minimizing |||I − F^{−1}F_{Y} |||_{2} subject to D_{p} − F_{Y} ≥ 0.

4 Preconditioned Conjugate Gradient Algorithm

When the matrix F_{Y} in the equation F_{Y} \underline{x} = \underline{e}_{1} is square, symmetric, and positive definite, the preconditioned conjugate gradient (CG) algorithm can be used to approximate the solution \underline{x} [8, Sec. 10.3]. Thus we can apply CG algorithms to approximating a column of the matrix CR bound via solution of (3). The CG algorithm can be interpreted as resulting from non-monotone and non-stationary acceleration of the splitting algorithm (4) via the introduction of time-varying acceleration parameters. [8, sec 10.3.6]. The CG algorithm converges to the exact solution in n iterations when run with infinite precision arithmetic. However, when run to termination it is not computationally competitive with Gaussian elimination. We will show that with proper preconditioning matrix F the following prematurely stopped preconditioned CG algorithm [8, Algorithm 10.3.1] can have significantly faster convergence than the monotone splitting iterations.

CG RECURSION FOR CR BOUND INITIALIZE: \underline{β}^{(0)} = \underline{0}; \underline{r}^{(0)} = \underline{e}_{1} RECURSION:

Solve: F\underline{z}^{(k)} = \underline{r}^{(k)}

α^{(k)} = \{ \begin{matrix}0, & k = 0 \\ \frac{<\underline{r}^{(k)},\underline{z}^{(k)}>}{<\underline{r}^{(k−1)},\underline{z}^{(k−1)}>}, & k > 0 \end{matrix}

\underline{p}^{(k)} = \{ \begin{matrix}\underline{z}^{(0)}, & k = 0 \\ \underline{z}^{(k)} + α^{(k)} \underline{p}^{(k−1)}, & k > 0 \end{matrix}
λ^{(k)} = \frac{< \underline{r}^{(k)}, \underline{z}^{(k)} >}{< \underline{p}^{(k)}, F_{Y} \underline{p}^{(k)} >}
\underline{r}^{(k+1)} = \underline{r}^{(k)} − λ^{(k)} F_{Y} \underline{p}^{(k)}
\underline{β}^{(k+1)} = \underline{β}^{(k)} + λ^{(k)} \underline{p}^{(k)}

CR BOUND APPROXIMATION:

η^{(k)} = \underline{e}_{1}^{T} \underline{β}^{(k)}.  (10)

When the preconditioner F is a banded p-diagonal matrix the CG algorithm requires the same number of flops per iteration (2n^{2} + 2np^{2}) as the splitting algorithms previously described. In the CG algorithm \underline{r}^{(k)} is the forward residual \underline{r}^{(k)} = \underline{e}_{1} − F_{Y} \underline{β}^{(k)} = F_{Y} Δ\underline{β}^{(k)} where Δ\underline{beta}^{(k)} is the approximation error Δ\underline{β}^{(k)}F_{Y}^{−1}\underline{e}_{1} − \underline{β}^{(k)}. While it has long been oberved that the speed of convergence of preconditioned CG generally improves as the eigenvalue spread of F^{−1}F_{Y} decreases, only a bound on convergence rate of the norm of the residuals ∥\underline{r}^{(k)}∥_{2} = ∥\underline{β}^{(k)} − F_{Y}^{−1}\underline{e}_{1}∥_{F_{Y}} is available, where ∥\underline{u}∥_{F_{Y}}^{2} = \underline{u}^{T} F_{Y} \underline{u}.

From [4, Sec. 2.3.1] we have ∥Δ\underline{β}^{(k)}∥_{F_{Y}} ≤

2∥Δ\underline{β}^{(0)}∥_{F_{Y}} (\frac{\sqrt{κ}−1}{\sqrt{κ}+1})^{k} where κ is the spectral condition

number of the matrix F^{−1}F_{Y}, defined as the ratio of the largest to the smallest magnitude eigenvalues. If one adopts ∥Δ\underline{β}^{(k)}∥_{F_{Y}} as an appropriate measure of approximation error then, by analogous reasoning as used for convergence analysis of the splitting algorithms, we obtain the following bound on the number of iterations N(γ) necessary to guarantee ∥Δ\underline{β}^{(k)}∥_{F_{Y}} converges to within a factor γ: N^{u}(γ) = \frac{log γ−log(2)}{log(\sqrt{κ}−1)−log(\sqrt{κ}+1)}. Combining this N^{u} expression with the per iteration cost of the algorithm (10) yields a similar inequality as in (8) which must be satisfied if preconditioned CG is to be computationally advantageous as compared to Gaussian elimination methods. Offline estimates of κ can be implemented (O(n^{2}) additional flops) using power iterations or iterative sup-norm approximations to the maximum and minimum eigenvalues. [8, Sec. 3.5.4]. Alternatively, on-line estimates of these extremal eigenvalues can be updated in O(1) flops from quantities produced in the preconditioned CG recursion [4, Sec. 5.1].

5 CRB Approximation for Singular FIM

The splitting iterations and CG algorithm described in the previous section are only applicable to non-singular Fisher information matrices F_{Y}. If it is known that \underline{e}_{1} lies in the range space of F_{Y} the CR bound \underline{e}_{1}^{T} F_{Y}^{+}\underline{e}_{1} for singular F_{Y} can be found in 4n^{3}/3 flops using the QR factorization to solve for the min-norm solution \underline{x} to F_{Y} \underline{x} = \underline{e}_{1} [8, Alg. 5.7.2]. However, typically the range space of F_{Y} is unknown and much more computationally intensive algorithms are required, e.g. the singular value decomposition (SVD) (20n^{3} flops). Here we present an iterative approximation to the singular CR bound \underline{e}_{1}^{T} F_{Y}^{+}\underline{e}_{1}.

Consider the following matrix difference

G(ϵ) =^{def} \frac{a}{a − b}(F_{Y} + a ϵI)^{−1} − \frac{b}{a − b}(F_{Y} + b ϵI)^{−1}, (11)

where ϵ > 0 and a, b ∈ (0, ∞), a ≠ b, are free parameters. G(ϵ) is a convergent approximation to the pseudo-inverse of F_{Y} in the sense

F_{Y}^{+} = lim_{ϵ→0^{+}} G(ϵ).  (12)

The representation (12) can be easily established by considering the eigendecomposition of F_{Y}

F_{Y}^{+} − G(ϵ) = ∑_{i=1}^{r} \frac{1}{σ_{i}} \underline{u}_{i}\underline{u}_{i}^{T} − ∑_{i=1}^{r} \frac{a}{(a − b)(σ_{i} + aϵ)}\underline{u}_{i}\underline{u}_{i}^{T}

+ ∑_{i=1}^{r} \frac{b}{(a − b)(σ_{i} + bϵ)}\underline{u}_{i}\underline{u}_{i}^{T}

= ∑_{i=1}^{r} \frac{ϵ (a σ_{i} + b σ_{i} + a b ϵ)}{σ_{i} (σ_{i} + a ϵ) (σ_{i} + b ϵ)}\underline{u}_{i}\underline{u}_{i}^{T} , (13)

where σ_{1} ≥, ..., ≥ σ_{r} > 0 are the r non-zero eigenvalues of F_{Y} arranged in decreasing order, and \{\underline{u}_{i}\}_{i=1}^{n} is an orthonormal set of eigenvectors. Observe that the range space of G(ϵ) corresponds to the range space of F_{Y} for all ϵ > 0.

Note that for all ϵ > 0, F_{Y}^{+} − G(ϵ) ≥ 0 so that \underline{e}_{1}^{T} G(ϵ)\underline{e}_{1} is itself a valid lower bound on var\underline{θ}(θ̂_{1}) which converges to the CR bound \underline{e}_{1}^{T} F_{Y}^{+}\underline{e}_{1} as ϵ → 0. An iterative approximation to \underline{e}_{1}^{T} G(ϵ)\underline{e}_{1} is obtained from the weighted difference of iterative approximations to the solutions of the linear equations [F_{Y} + aϵI]γ_{1} = \underline{e}_{1} and [F_{Y} + bϵI]γ_{2} = \underline{e}_{1}. In particular, using the following assignment for (4) we obtain a splitting algorithm for approximating the singular CR bound

F_{1}\underline{γ}_{1}^{(k+1)} = [F_{1} − F_{Y} − aϵI] \underline{γ}_{1}^{(k)} − \underline{e}_{1}
F_{2}\underline{γ}_{2}^{(k+1)} = [F_{2} − F_{Y} − bϵI] \underline{γ}_{2}^{(k)} − \underline{e}_{1}
\underline{β}^{(k)} = \frac{a}{a − b}\underline{γ}_{1}^{(k)} − \frac{b}{a − b}\underline{γ}_{2}^{(k)}.  (14)

For large n the computational burden of the iteration (14) is 4n^{2}+2np^{2} flops which, for p ≪ n, which is approximately twice as expensive as the non-singular CR bound iteration (4).

Both the speed of convergence and the normalized asymptotic approximation error δ = \underline{e}_{1}^{T} (F_{Y}^{+} − G(ϵ))\underline{e}_{1}/\underline{e}_{1}^{T} F_{Y}^{+}\underline{e}_{1} increase as ϵ increases. It is easily shown that for ϵ ≪ min_{i=1,...,r}\{σ_{i}\}: \underline{e}_{1}^{T} (F_{Y}^{+} − G(ϵ))\underline{e}_{1} = ϵ ∑_{i=1}^{r}[(a + b) σ_{i} + a b ϵ]/[σ_{i} (σ_{i} + a ϵ) (σ_{i} + b ϵ)]∥\underline{u}_{i}^{T} \underline{e}_{1}∥^{2} ≈ ϵ(a + b)∥F_{Y}^{+}\underline{e}_{1}∥^{2}. Hence the normalized asymptotic approximation error is δ ≈ ϵ(a + b)∥F_{Y}^{+}\underline{e}_{1}∥^{2}/\underline{e}_{1}^{T} F_{Y}^{+}\underline{e}_{1}. The right hand side of this relation can be bounded from above and below (see Appendix) to yield a bound on asymptotic normalized error δ valid for ϵ ≪ min_{i=1,...,r}{σ_{i}}

ϵ(a + b)\frac{f_{11}}{∑_{i=1}^{n} f_{1i}^{2}} ≤ δ ≤ ϵ(a + b) ∑_{i=1}^{n} \frac{1}{f_{ii}} .  (15)

Note that the bounds depend on the magnitudes of ϵ, a, and b only through the product ϵ(a + b). In the following Section relation (15) will be used to select appropriate values of these parameters given a desired magnitude of normalized error δ via the equivalent relation

\frac{δ}{(a + b) ∑_{f_{ii}≠0} \frac{1}{f_{ii}}} ≤ ϵ ≤ \frac{δ ∑_{j=1}^{n} f_{ij}}{f_{11}(a + b)} .  (16)

6 Application to ECT Inverse Problem

We will briefly illustrate the CR bound algorithms for the following inverse problem arising in collimated emission computed tomography (ECT). One wishes to estimate the first component θ_{1} of the vector \underline{θ} in the linear model

\underline{Y} = A\underline{θ} + \underline{w},  (17)

where A is an m × n matrix, and \underline{w} is an additive noise. Here A is a sparse spatially variant system matrix with ρmn nonzero entries, ρ ≪ 1, and \underline{w} is a Poisson distributed noise. The vector \underline{θ} is a vector of unknown image intensity values indexed over n pixels, \underline{Y} is a vector of photon counts detected at each of m detector elements, and the ij-th element of the matrix A are the probabilities that a photon emitted at pixel index j is detected at detector element i. Note that m is typically several times larger than n in tomography applications. Under the standard assumption that for each detector the mean detected counts \underline{μ} = A\underline{θ} are strictly greater than zero the FIM is non-singular and has the form [1]

F_{Y} (\underline{θ}) = A^{T} diag (\underline{μ})^{−1} A.  (18)

Due to the simple product form (18) of F_{Y} the implementation of the generic CR bound algorithms presented previous sections can be simplified. In each of these algorithms the most intensive computation reduces to a vector matrix multiply of the form A^{T} \underline{x}^{(k)} (see algorithm descriptions in Fig. ??) requiring only on the order of 2mnρ flops per iteration. As a result of this reduction the Fisher matrix F_{Y} need not be precomputed or stored. This is significant since the computation of F_{Y} requires O(ρmn^{2}) flops and O(n^{2}) memory storage. Implementation of the iterative CR bound algorithms presented in previous sections only require precomputation of \underline{μ} = A\underline{theta} (2ρmn flops) and storage of ρmn non-zero elements of A and m elements of \underline{μ}. For matrix sparsity factors ρ < 0.01 commonly encountered in collimated ECT this corresponds to a substantial savings in storage and precomputation requirements.

For these studies the 32 × 32 Hoffman brain phantom shown in Fig. ??) was used as the true image. It was assumed that the overall size, location, and orientation of the brain within the image was known. Thus for computation of the bound the unknown parameter vector \underline{θ} consisted of a lexicographical ordering of the n = 640 pixel intensities within the ellipsoidal brain boundary. A system matrix A corresponding to axially collimated PET was constructed which acquires projections of the planar phantom over 40 detector angles and 80 radial detector bins. This gave a A with 3600 rows and 640 columns and sparsity factor ρ = 0.1.

While the CR bound approximation algorithms were investigated for many different pixels, we only present results for a representative case: the pixel located at coordinates (16, 21) indicated on Fig. ??. It was found that the convergence trajectories of the CG and GS algorithms were relatively insensitive to the pixel location. Figs. 1-3 show the trajectories of various iterative CR bound approximation algorithms for this pixels.

Fig. 1 shows the monotone splitting algorithms which are mathematically guaranteed to improve the bound approximation at each successive iteration. MD1 denotes the Figure 1: The monotone splitting algorithms for approximating the non-singular CR bound at pixel (16, 21). Straight line denotes the true value of the CR bound. MD1 uses the diagonal splitting matrix D_{1} defined in Sec. ??; MD2 uses the tridiagonal splitting matrix D_{2}; MD50 uses the 99-diagonal matrix D_{50}; and EM uses the diagonal splitting matrix derived in [?]

Figure 2: The non-monotone splitting algorithms for approximating the non-singular CR bound at pixel (16, 21). MDOR1, MDOR2 and MDOR50 are over-relaxed versions of MD1, MD2 and MD50, where the over-relaxation parameter was numerically optimized to maximize asymptotic convergence. Also shown are the Jacobi (J) and Gauss-Seidel (GS) algorithms.

Figure 3: The Gauss-Seidel (GS) and preconditioned conjugate gradient algorithms for approximating the nonsingular CR bound at pixel (16, 21). CGJ uses the standard diagonal Jacobi preconditioner, CGD1 uses the diagonal matrix D_{1} preconditioner given in Sec. ??, and CGDF uses a preconditioner tailored for the A matrix arising in the spatially invariant ECT geometry considered here.

monotone algorithm which uses the diagonal splitting matrix D_{1} defined in Sec. ??; MD2 uses the tridiagonal splitting matrix D_{2}; MD50 uses the 99-diagonal matrix D_{50}; and EM uses the diagonal splitting matrix derived in [1]. As expected, increasing the number 2k − 1 of diagonals in the splitting matrix D_{k} improves the convergence rate. Figure 2 shows several splitting algorithms which use splitting matrices which do not guarantee monotonicity. In the figure MDOR1, MDOR2 and MDOR50 denote overrelaxed versions of MD1, MD2 and MD50, where the overrelaxation parameter was numerically optimized to maximize asymptotic convergence rate. Also shown are Gauss-Seidel (GS) and relaxed Jacobi (J), with relaxation parameter ψ numerically determined to ensure convergence and maximize asymptotic convergence rate. The standard unrelaxed Jacobi algorithm diverged for all cases studied and is not shown. Note that GS exhibits much faster convergence than any of the other splitting algorithms. Finally several preconditioned CG algorithms are shown in Figs.

3. In the figure the CGJ algorithm uses the standard diagonal Jacobi preconditioner, i.e. F = diag(F_{Y}) the diagonal elements of the FIM, CGD1 uses the diagonal matrix D_{1} given in Sec. ??) as preconditioner. The algorithm denoted CGDF uses a special preconditioner F consisting of a diagonal preconditioner, which makes F_{Y} almost circulant, followed by a Fourier-type preconditioner. The preconditioner used in CGDF is tailored to the spatially invariant PET application and is described in [?] in the context of fast least squares PET reconstruction algorithms. Also plotted for comparison is the GS algorithm. All of the algorithms shown in In Figs. 3 and ?? have much faster convergence rate than the monotonic splitting algorithms; compare CGDF which converges to within 5% of the true CR bound in only 2 iterations while MD50 takes more than 200 iterations to converge to within the same accuracy. It is interesting to note that the GS and CG algorithms are almost monotonic, exhibiting only a very small amount of overshoot.

Convergence properties of the splitting algorithms are quantified in Table ?? and those of the conjugate gradient algorithms are illustrated in ??. The first 3 columns of Table ?? contain the matrix-2 norm |||Δ|||_{2}, the simpler to compute Frobenius norm difference |||F − F_{Y} |||_{F}, and the root convergence factor λ_{max}^{Δ} for the matrix I − F^{−1}F_{Y}. The 4th column of Table ?? is the actual number of iterations and total number of required Mflops to convergence to within δ = 0.05, or 5% of the limit. Note, as discussed in Section ?? the matrix-2 norm and Frobenius norms are better predictors of the actual time to converge than is the root convergence factor. In particular, observe that the Frobenius norm correctly predicts that MD1 converges in fewer iterations than the algorithm using the splitting matrix F = D_{λ} which has minimum root convergence factor.

The first column of Table ?? contains the ratio \frac{\sqrt{κ}−1}{\sqrt{κ}+1} where κ is the spectral condition number of the preconditioned FIM F^{−1}F_{Y}. The second and third columns of Table ?? are the number of iterations and the number of flops required to converge within δ = 0.05 of the limit.

Next we turn to the case of singular FIM arising in the so-called “missing angle problem.” In this problem the system must estimate image parameters with a greatly reduced number and range of projections. For this study only 10 angles from 0 to π/4 and 40 radial bins per angle are available; leading to decimating A (3) by a factor of eight. This resulted in a A matrix of dimension 400 × 640 with rank 400 and condition number on the order of 1000. We implemented the perturbed differencing algorithm (14) to approximate the first column of the pseudo-inverse of F_{Y}. It was found that over the range 0.001 ≤ \frac{a−b}{a+b} ≤ 1000 the choice of a and b has only a minor effect on the convergence rate; here we only show results for a − b = 0.0001 which appeared to maximize the convergence rate of all algorithms. The maximum allowable asymptotic normalized error was set to δ = 0.05 or 5% and epsilon was selected according to (16) as the average of the lower and upper

bounds: ϵ = \frac{1}{2} \frac{δ}{a+b} (_{\frac{1}{∑_{f_{ii}≠0} \frac{1}{f_{ii}}}} + ^{\frac{∑_{j=1}^{n} f_{ij}}{f_{11}}} ). Figure ?? il-

lustrates the closeness of this selection rule to the true epsilon which exactly meets the constraint δ on asymptotic normalized error. Figure ?? illustrates the convergence trajectories of the GS and CG1 and CGD algorithms applied to the pair of iterations (14). It is interesting to note that while the CG algorithms converged most rapidly on their respective difference terms \underline{β}_{1}^{(k)} and \underline{β}_{2}^{(k)}, after performing the differencing \frac{a}{a−b} \underline{β}_{1}^{(k)} − \frac{b}{a−b} \underline{β}_{2}^{(k)}, the GS algorithm converged the most rapidly. Finally, while space limitations prevent showing any numerical results, it was observed that a significant tradeoff exists between the convergence rate and the asymptotic normalized error δ.

7 Conclusion

We have studied the convergence properties of several algorithms for recursive approximation of the CR bound: including monotone convergent splitting matrix iterations; non-montonic over-relaxed splitting matrix iterations; and non-monotonic preconditioned conjugate gradient iterations. Not surprisingly, to obtain similar convergence rates the computational cost of a monotonically convergent algorithm is quite high as compared with a non-monotonic algorithm in the context of the inverse problems considered. The general conclusion of this study is that unless monotonic convergence is an absolute necessity the preconditioned conjugate gradient and Gauss-Seidel algorithms should be seriously considered for iterative CR bound computation due to their rapid convergence rates.

Appendix: Bounds on Asymptotic Error

Lemma 2 Let δ be the normalized asymptotic error defined as δ = \underline{e}_{1}^{T} (F_{Y}^{+} − G(ϵ))\underline{e}_{1}/\underline{e}_{1}^{T} F_{Y}^{+}\underline{e}_{1} defined in Sec. 5. Assume that ϵ is dominated by the smallest positive eigen-

value of F so that to o(ϵ): δ = ϵ(a + b)\frac{\underline{e}_{1}^{T} [F_{Y}^{+}]^{2}\underline{e}_{1}}{\underline{e}_{1}^{T} F_{Y}^{+} \underline{e}_{1}} . Then

ϵ(a + b)\frac{f_{11}}{∑_{i=1}^{n} f_{1i}^{2}} ≤ δ ≤ ϵ(a + b) ∑_{f_{ii}≠0} \frac{1}{f_{ii}} ,  (19)

Proof: We first show the lower inequality. If f_{11} = 0 then this inequality is trivial so hencforth we assume that f_{11} >

0. The Cauchy-Schwarz inequality states that |\underline{u}^{T} \underline{v}|^{2} ≤ \underline{u}^{T} \underline{u} · \underline{v}^{T} \underline{v} for any two vectors \underline{u} and \underline{v}. Letting \underline{u} = F_{Y}^{k/2}\underline{w} and \underline{v} = F_{Y}^{k/2+1}\underline{w} we obtain

\begin{matrix}\underline{w} & T & F & {k+1 \atop Y} & \underline{w} \\ \underline{w} & T & F & {k+2 \atop Y} & \underline{w} \end{matrix} ≤ \begin{matrix}\underline{w} & T & F & k Y & \underline{w} \\ \underline{w} & T & F & {k+1 \atop Y} & \underline{w} \end{matrix} .

Setting \underline{w} = F_{Y}^{+}\underline{e}_{1} and applying the above for k = 0, 1, 2 we obtain the sequence of inequalities

\begin{matrix}\underline{e} & T 1 & F & Y & \underline{e} & 1 \\ \underline{e} & T 1 & F & 2 Y & \underline{e} & 1 \end{matrix} ≤ \frac{\underline{e}_{1}^{T} P_{F_{Y}} \underline{e}_{1}}{\underline{e}_{1}^{T} F_{Y} \underline{e}_{1}} ≤ \frac{\underline{e}_{1}^{T} F_{Y}^{+}\underline{e}_{1}}{\underline{e}_{1}^{T} P_{F_{Y}} \underline{e}_{1}} ≤ \frac{\underline{e}_{1}^{T} [F_{Y}^{+}]^{2}\underline{e}_{1}}{\underline{e}_{1}^{T} F_{Y}^{+}\underline{e}_{1}} ,

where P_{F_{Y}} = F_{Y}^{+}F_{Y} = F_{Y} F_{Y}^{+} is a symmetric idempotent matrix which projects vectors onto the column space of F_{Y}. Since \underline{e}_{1}^{T} F_{Y} \underline{e}_{1} = f_{11} and \underline{e}_{1}^{T} F_{Y}^{2} \underline{e}_{1} = ∥F_{Y} \underline{e}_{1}∥^{2} = ∑_{i=1}^{n} f_{1i}^{2} we have established the lower inequality in (19).

Now turn to the upper inequality in (19). Since ∑_{k=1}^{n} \underline{e}_{k}\underline{e}_{k}^{T} = I we have:

\frac{\underline{e}_{1}^{T} [F_{Y}^{+}]^{2}\underline{e}_{1}}{\underline{e}_{1}^{T} F_{Y}^{+}\underline{e}_{1}} = \frac{∑_{k=1}^{n} |\underline{e}_{1}^{T} F_{Y}^{+}\underline{e}_{k}|^{2}}{\underline{e}_{1}^{T} F_{Y}^{+}\underline{e}_{1}}

= \underline{e}_{1}^{T} F_{Y}^{+}\underline{e}_{1} + ∑_{k=2}^{n} \frac{|\underline{e}_{1}^{T} F_{Y}^{+}\underline{e}_{k}|^{2}}{\underline{e}_{1}^{T} F_{Y}^{+}\underline{e}_{1}}

≤ \underline{e}_{1}^{T} F_{Y}^{+}\underline{e}_{1} + ∑_{k=2}^{n} \underline{e}_{k}^{T} F_{Y}^{+}\underline{e}_{k}

= ∑_{k=1}^{n} \underline{e}_{k}^{T} F_{Y}^{+}\underline{e}_{k},  (20)

where the Cauchy-Schwarz inequality was used in the third line. Let P_{k} be a permutation matrix such that F_{k} = P_{k}^{T} F_{Y} P_{k} is the matrix resulting from exchanging the k-th and n-th rows and columns of F_{Y}. Let F_{k} be partitioned as

F_{k} = [ \begin{matrix}A & \underline{c} \\ \underline{c}^{T} & α \end{matrix} ] ,

where α = f_{kk} is the k-th diagonal element of F_{Y}, A is an (n−1)×(n−1) matrix and \underline{c} is an (n−1)-element column vector. Next we recall the formula for the pseudoinverse of partitioned non-negative definite matrices [11]

F_{k}^{+} = [ \begin{matrix}A^{+} + q^{+}A^{+}\underline{c}\underline{c}^{T} A^{+} & −q^{+}A^{+}\underline{c} \\ −q^{+}\underline{c}^{T} A^{+} & q^{+} \end{matrix} ] ,

where q = f_{kk} − \underline{c}^{T} A^{+}\underline{c} ≥ 0. Hence \underline{e}_{k}^{T} F_{Y}^{+}\underline{e}_{k} = \underline{e}_{n}^{T} F_{k}^{+}\underline{e}_{n} = q^{+}. If f_{kk} = 0 then q = q^{+} = 0. Otherwise, since A^{+} is non-negative definite, q^{+} = \frac{1}{f_{kk}−\underline{c}^{T} A^{+}\underline{c}} ≤ \frac{1}{f_{kk}}. Plugging

this into (20) we obtain the upper inequality of (19) which finishes the proof of the lemma. □

References

[1] A. O. Hero and J. A. Fessler, “A recursive algorithm for computing CR-type bounds on estimator covariance,” IEEE Trans. on Inform. Theory, vol. 40, pp. 1205–1210, July 1994.

[2] A. R. Kuruc, “Lower bounds on multiple-source direction finding in the presence of direction-dependent antenna-array-calibration errors,” Technical Report 799, M.I.T. Lincoln Laboratory, Oct., 1989.

[3] A. O. Hero and J. A. Fessler, “A fast recursive algorithm for computing CR-type bounds for image reconstruction problems,” in Proc. of IEEE Nuclear Science Symposium, pp. 1188–1190, Orlando, FA, Oct. 1992.

[4] R. Barrett and etal, Templates for the solution of linear systems: building blocks for iterative methods, SIAM Press, Philadelphia, 1994.

[5] C. R. Rao, Linear Statistical Inference and Its Applications, Wiley, New York, 1973.

[6] J. D. Gorman and A. O. Hero, “Lower bounds for parametric estimation with constraints,” IEEE Trans. on Inform. Theory, vol. IT-36, pp. 1285–1301, Nov. 1990.

[7] F. A. Graybill, Matrices with Applications in Statistics, Wadsworth Publishing Co., Belmont CA, 1983.

[8] G. H. Golub and C. F. Van Loan, Matrix Computations (2nd Edition), The Johns Hopkins University Press, Baltimore, 1989.

[9] R. A. Horn and C. R. Johnson, Matrix Analysis, Cambridge, 1985.

[10] D. M. Young, Iterative solution of large linear systems, Academic Press, New York, 1971.

[11] C. A. Rohde, “Generalized inverses of partitioned matrices,” J. Soc. Indust. Appl. Math. (SIAM), vol. 13, pp. 1033–1035, 1965.