1 Introduction↩︎

In numerous applications in science and engineering it is necessary to solve systems of linear equations with interest in only portions of the solution vector. These problems can be formulated as follows: Given a linear model \[b = {\mathsf A}x + \varepsilon,\] where \({\mathsf A}\in{\mathbb{R}}^{m\times n}\) is a known matrix, \(b\) is a given vector, and \(\varepsilon\) is possibly vanishing random additive noise, estimate selected components of the vector \(x\in{\mathbb{R}}^n\). Without loss of generality, we can permute the entries of \(x\) so that the quantity of interest consists of its first \(n_1\) components, \(n_1<n\). By partitioning the matrix \({\mathsf A}\) accordingly, we can reformulate the problem as \[\label{original} b = {\mathsf A}_1 x_1 + {\mathsf A}_2 x_2 + \varepsilon,\tag{1}\] where \({\mathsf A}_j \in {\mathbb{R}}^{m\times n_j}\), \(j=1,2\), are known matrices, \(n_1+n_2 = n\), and \(x_j\in{\mathbb{R}}^{n_j}\) are vectors, \[{\mathsf A}= \left[\begin{array}{cc} {\mathsf A}_1 & {\mathsf A}_2\end{array}\right], \quad x = \left[\begin{array}{c} x_1 \\ x_2\end{array}\right].\]

Applications that are amenable to this type of formulation include inverse problems (\(\varepsilon \neq 0\)) in which \(b\) represents a noisy observation, and forward problems (\(\varepsilon = 0\)) where the matrix \({\mathsf A}\) is typically a square and invertible stiffness matrix arising from an elliptic partial differential equation and \(b\) is a known load vector. In the latter case, one may be interested only in a selected subset of the components of the solution \(x\), discarding the components of no relevance. In inverse problems, computational tasks of the form (1 ) arise in many different contexts, including the following important applications.

Local tomography: In local X-ray tomography, the measured X-ray attenuation data depends on the absorption properties of the body extending over the entire field of view of the imaging system, although only a portion of the domain, often referred to as region of interest (ROI) corresponding to \(x_1\), may be of interest. Local tomography commonly arises in X-ray tomography applications, either because of limited detector size, e.g., in dental Cone Beam Computed Tomography (CBCT) or because of the need to reduce the X-ray dose outside of the ROI, e.g., in cardiac imaging.

Radar imaging: In radar imaging, the signal comprises the portion of interest coming from the main lobe of the radar, corresponding to \(x_1\), and the clutter from the radar sidelobes that pick up the echo of trees, ground etc., represented by \(x_2\).

Image deblurring: A blurred image typically accounts for the contribution from both the field of view and leak data from the outside region. Inappropriate modeling or ignoring of the latter may cause boundary artifacts that are often compensated for by means of artificial boundary conditions.

MEG/EEG based brain imaging: The search for the onset foci of epileptic seizures often can be restricted to a portion of the brain; Likewise, localization of active brain regions constituting a task-related, or task-negative connected network may be limited to selected brain regions, however the data account for contributions from the entire brain.

Functional MRI imaging: Identification of a stimulus-related cerebral activity of interest typically may be circumscribed to a known part of the brain, and the dynamics of the rest of the brain and other tissues is of less or no interest.

In statistics, the portion of the unknown \(x_2\) in (1 ) not of interest is referred to as nuisance parameters while the radar community refer to the contribution \({\mathsf A}_2 x_2\) to the noiseless signal \({\mathsf A}x\) as clutter. In this paper we adopt the terminology from both these two communities.

In the inverse problems literature, the separating the portion of interest from the clutter has been addressed in different ways. The tools for nuisance parameter elimination in classical statistics include considerations of sufficient and ancillary statistics and likelihood ratios [1], [2]. Bayesian statistics addresses the problem typically through integrated likelihood methods [3] or marginalization over the nuisance parameter, often requiring MCMC sampling [4]. An example of optimization-based methods for estimating nuisance parameters in inverse problems can be found in [5]. In a naı̈ve approach, the nuisance parameter \(x_2\) can be fixed to a plausible or representative value, or ignored by setting \(x_2=0\), typically introducing a bias that in the Bayesian setting can be compensated by appropriate error modeling using the prior distribution of \(x_2\) [6].

In this article, we put the spotlight of the inversion process on the portion \(x_1\) of the unknown using a linear algebra-based approach exploiting suitable orthogonal projections. Let \[{\mathsf P}:{\mathbb{R}}^m\to {\mathcal{R}}({\mathsf A}_2), \quad {\mathsf P}^\perp:{\mathbb{R}}^m\to {\mathcal{R}}({\mathsf A}_2)^\perp = {\mathcal{N}}({\mathsf A}_2^{\mathsf T}),\] denote the pair of orthogonal projectors onto the range \({\mathcal{R}}({\mathsf A}_2)\subset{\mathbb{R}}^m\) of the matrix \({\mathsf A}_2\) and onto the the null space \({\mathcal{N}}({\mathsf A}_2^{\mathsf T})\) of the matrix \({\mathsf A}_2^{\mathsf T}\), respectively . Assuming that \(m>n_2\), the subspace \({\mathcal{N}}({\mathsf A}_2^{\mathsf T})\) is non-trivial. Applying \({\mathsf P}^\perp\) to both sides of (1 ) and observing that \({\mathsf P}^\perp{\mathsf A}_2 x_2 =0\) yields \[\label{projected} b' = {\mathsf A}_1' x_1 + \varepsilon',\tag{2}\] where \(b' = {\mathsf P}^\perp b\), \({\mathsf A}_1' = {\mathsf P}^\perp{\mathsf A}_1\), \(\varepsilon' = {\mathsf P}^\perp \varepsilon\). The projected problem, which does not contain the nuisance parameters can be solved as a least squares problem, equipped with a standard regularization technique if needed, or alternatively it can be analyzed in the framework of Bayesian inverse problems.

In section 2, we discuss the spotlighting for well-posed problems, starting with an \(n\times n\) invertible matrix \({\mathsf A}\), and we show that the approach reduces the original problem into an overdetermined least squares problem that recovers the selected components \(x_1\) exactly.

In section 3, the spotlighting idea is applied to inverse problems where the dimensionality of the nuisance parameter does not exceed that of the data, \(m>n_2\) In particular, we show that if the inverse problem is considered in a Bayesian setting, the proposed projection method for reducing the problem dimensionality does not correspond to the standard dimension reduction by marginalization in the Bayesian framework of inverse problems, but can be interpreted as a result of a limiting process as the underlying prior becomes uninformative. We emphasize that the clutter removal method by projection does not require or assume a Bayesian interpretation of the problem.

The assumption that the clutter is restricted to a proper subspace of \({\mathbb{R}}^m\), is automatically satisfied in the case of well-posed invertible problems \(m=n\) as well as in overdetermined and formally determined inverse problems where \(m\geq n>n_2\). In inverse problems, the condition \(m>n_2\) is not always guaranteed, and may not be satisfied for commonly encountered underdetermined problems. In such cases, by replacing the exact projector \({\mathsf P}^\perp\) by a low-rank approximation, spotlighting may still be a viable clutter reduction method, as will be discussed in 3.3. Extensions to non-linear inverse problems are briefly discussed in section 5.

In the inverse problems literature, dimension reduction techniques have been discussed in the Bayesian context, see, e.g., [7], [8], based on the analysis of the informative directions of the likelihood with respect to the prior information. Unlike the cited work, which assume a Bayesian framework, our dimension reduction for inverse problems is purely linear algebraic and does not require a stochastic extension of the model. In particular, we need not assume any prior model for the unknown, and in the noiseless case \(\varepsilon = 0\), even the likelihood is a futile concept. A projection technique similar to the one proposed has been proposed in statistical regression analysis and in particular, in econometrics, referred to as partial regression, or (Yule) -Frisch-Waugh-Lovell ((Y)FWL) theorem [9] dating back to the early 20th century; see [10] for a historical account. The present article aims at bringing the idea into a computational linear algebraic context for both well-posed linear systems and ill-posed inverse problems.

2 Spotlight projections for Model Reduction in Noiseless Case↩︎

Considering a noiseless model (1 ), where \(\varepsilon = 0\), and \({\mathsf A}\in{\mathbb{R}}^{n\times n}\) is an invertible matrix, partitioned as \[{\mathsf A}= \left[\begin{array}{cc} {\mathsf A}_1 & {\mathsf A}_2\end{array}\right], \quad {\mathsf A}_1\in{\mathbb{R}}^{n\times k}, \quad {\mathsf A}_2 \in{\mathbb{R}}^{n\times(n-k)},\] where \(1\leq k<n\). Let \({\mathsf P}: {\mathbb{R}}^n \to {\mathcal{R}}({\mathsf A}_2)\) and \({\mathsf P}^\perp: {\mathbb{R}}^n \to {\mathcal{R}}({\mathsf A}_2)^\perp\) be the orthogonal projectors onto the range of \({\mathsf A}_2\) and its orthogonal complement, respectively. Consider the equation \[\label{full} {\mathsf A}x = b,\tag{3}\] and its projected version onto the orthogonal complement of the range of \({\mathsf A}_2\). \[\label{nsrpocub} {\mathsf A}_1' x_1 = b', \quad {\mathsf A}_1' = {\mathsf P}^\perp {\mathsf A}_1, \quad b' = {\mathsf P}^\perp b.\tag{4}\] We start by stating and proving the following result.

Proof. First, observe that by the invertibility of \({\mathsf A}\), its columns are linearly independent, hence \[{\rm dim}\big({\mathcal{R}}({\mathsf A}_1)\big) = k, \quad {\rm dim}\big({\mathcal{R}}({\mathsf A}_2)\big) = n-k.\] Moreover, \[{\mathbb{R}}^n = {\mathcal{R}}({\mathsf A}_1) \oplus {\mathcal{R}}({\mathsf A}_2),\] that is, the space \({\mathbb{R}}^n\) is the direct, not necessarily orthogonal sum of the range of \({\mathsf A}_1\) and the range of \({\mathsf A}_2\). To prove this, we need to show that

1. \({\mathbb{R}}^n = {\mathcal{R}}({\mathsf A}_1) + {\mathcal{R}}({\mathsf A}_2)\), i.e, every vector in \({\mathbb{R}}^n\) can be expressed as the sum of a vector in \({\mathcal{R}}({\mathsf A}_1)\) and a vector in \({\mathcal{R}}({\mathsf A}_2)\).

2. \({\mathcal{R}}({\mathsf A}_1) \cap {\mathcal{R}}({\mathsf A}_2) = \{0\}\), implying that the representation is unique.

To prove the first claim, observe that the columns of \({\mathsf A}\) are linearly independent, hence a basis for \({\mathbb{R}}^n\), thus every \(y\in{\mathbb{R}}^n\) admits a unique representation of the form \[y = {\mathsf A}x = \left(\sum_{j=1}^k x_j a^{(j)}\right) + \left(\sum_{j=k+1}^n x_j a^{(j)}\right),\] where \(x = {\mathsf A}^{-1} y\). To show the second claim, observe that if \(z \in {\mathcal{R}}({\mathsf A}_1) \cap {\mathcal{R}}({\mathsf A}_2)\), it can be written in two ways, \[z = \left(\sum_{j=1}^k z_j a^{(j)}\right) = \left(\sum_{j=k+1}^n z_j a^{(j)}\right),\] which implies that the vector \[\widetilde{z} = \left[\begin{array}{r} z_1 \\ \vdots \\ z_k \\ -z_{k+1} \\ \vdots \\ -z_n\end{array}\right],\] satisfies \({\mathsf A}\widetilde{z} = 0\), so by the invertibility of \({\mathsf A}\), \(\widetilde{z} = 0\) and hence \(z=0\).

To show that the matrix \({\mathsf A}_1'\) has full rank it suffices to show that its null space is trivial.

If \(z\in{\mathbb{R}}^k\) is a vector in the null space of \({\mathsf A}_1'\), then \[{\mathsf A}_1' z = {\mathsf P}^\perp {\mathsf A}_1 z = 0.\] If \[{\mathsf A}_1 z = \sum_{j=1}^k z_j a^{(j)},\] then \[{\mathsf P}^\perp\left(\sum_{j=1}^k z_j a^{(j)}\right) = 0,\] hence \[\sum_{j=1}^k z_j a^{(j)} \in {\mathcal{N}}({\mathsf P}^\perp) = {\mathcal{R}}({\mathsf A}_2)\] and \[\sum_{j=1}^k z_j a^{(j)} \in {\mathcal{R}}({\mathsf A}_1)\cap {\mathcal{R}}({\mathsf A}_2) = \{0\}.\] The linear independence of the columns of \({\mathsf A}_1\) implies that \(z = 0\), hence \({\mathsf A}_1'\) is full rank.

The full rank of \({\mathsf A}_1'\) implies that the projected equation in the spotlight, (2 ) has the unique (least squares) solution \[x_1' = \left({\mathsf A}_1'\right)^\dagger b' = \left(({\mathsf A}_1')^{\mathsf T}{\mathsf A}_1'\right)^{-1}\left({\mathsf A}_1'\right)^{\mathsf T}b'.\] Let \[x = \left[\begin{array}{c} x_1 \\ x_2 \end{array}\right]\] be the unique solution of the original linear system (3 ). Its orthogonal projection onto the orthogonal complement of the range of \({\mathsf A}_ 2\) satisfies \[{\mathsf A}'_1 x_1 = b',\] and by the uniqueness of the solution it follows that \(x_1' = x_1\). \(\Box\)

3 Spotlight projections for inverse problems↩︎

In this section, we consider the case in which the vector \(b\) represents a noisy observation, and the matrix \({\mathsf A}\) is ill-conditioned. Consequently, the inverse problem is ill-posed and therefore it requires regularization, or alternatively, its stochastic extension can be considered in the Bayesian framework. The proposed approach is to consider the projected linear inverse problem \[b' = {\mathsf A}_1' x_1 + \varepsilon'\] as an ill-posed problem independent of the original one, introducing any standard regularization scheme, or, alternatively, a prior for the unknown \(x_1\). While the spotlight projection does not assume or require a Bayesian interpretation of the problem, it is natural to ask whether the methodology lends itself to a statistical interpretation. In the following subsection, we provide a Bayesian interpretation of spotlight projection method, under the assumption of Gaussian prior and Gaussian likelihood. We begin by outlining standard techniques to deal with nuisance parameters based on marginalization and modeling the clutter as part of the noise. For a general reference on the Bayesian solution of linear inverse problems, see, e.g., [11].

3.1 Marginalization and error modeling↩︎

In the Bayesian extension of the inverse problem (1 ), let \(X\) denote a \({\mathbb{R}}^n\)-valued Gaussian random variable with prior probability density \(X\sim{\mathcal{N}}(0,{\mathsf C})\), where \(n=n_1+n_2\), and \({\mathsf C}\in{\mathbb{R}}^{n\times n}\) is a symmetric positive definite covariance matrix. We partition \(X\) and the forward model \[X = \left[\begin{array}{c} X_1 \\ X_2\end{array}\right], \quad X_j\in{\mathbb{R}}^{n_j}, \quad {\mathsf A}= \left[\begin{array}{cc} {\mathsf A}_1 & {\mathsf A}_2\end{array}\right] \in {\mathbb{R}}^{m\times n},\] and write the stochastic extension of the model as \[B = {\mathsf A}X + E = {\mathsf A}_1 X_1 + {\mathsf A}_2 X_2 + E,\] where \(E\) is an additive Gaussian noise vector independent of \(X\). Without loss of generality, we may assume that the noise is whitened, i.e., \(E\sim{\mathcal{N}}(0,{\mathsf I})\). Furthermore, we partition the prior covariance matrix as \[{\mathsf C}= \left[\begin{array}{cc} {\mathsf C}_{11} & {\mathsf C}_{12} \\ {\mathsf C}_{21} & {\mathsf C}_{22}\end{array}\right], \quad {\mathsf C}_{jk} \in{\mathbb{R}}^{n_j\times n_k}.\]

The estimation of the variable \(X_1\) based on the observation \(B\), may be done in two different ways:

1. Lumping and conditioning: Write the forward model as \[B = {\mathsf A}_1 X_1 + \big({\mathsf A}_2 X_2 + E\big) = {\mathsf A}_1 X_1 + \widehat E,\] interpreting the clutter term as part of the observation noise \(\widehat E = {\mathsf A}_2 X_2 + E\), and then calculate directly the posterior density \(\pi_{X_1\mid B}(x_1\mid b)\) by conditioning.

2. Marginalization: Calculate the posterior probability density of \((X_1, X_2)\) \[\pi_{X\mid B}(x\mid b) = \pi_{X_1,X_2\mid B}(x_1,x_2\mid b),\] and then marginalize it with respect to \(X_2\), \[\pi_{X_1\mid B}(x_1\mid b) = \int \pi_{X_1,X_2\mid B}(x_1,x_2\mid b) {\rm d} x_2. \label{eq:marginalX1}\tag{5}\]

While it is intuitively clear that these two approaches should lead to the same posterior density, for completeness we state this as a theorem and give a short proof of it.

Proof. Conditioning approach. Combine \(X_1\) and \(B\) into the Gaussian random variable \[Z = \left[\begin{array}{c} X_1 \\ B\end{array}\right].\] From the observations that \[{\mathbb{E}}\big(X_1 X_1^{\mathsf T}\big) = {\mathsf C}_{11}, \quad {\mathbb{E}}\big(X_1 B^{\mathsf T}\big) = {\mathsf C}_{11}{\mathsf A}_1^{\mathsf T}+ {\mathsf C}_{12}{\mathsf A}_2^{\mathsf T}= {\mathbb{E}}\big(B X_1^{\mathsf T}\big)^{\mathsf T},\] and \[{\mathbb{E}}\big(B B^{\mathsf T}\big) = {\mathbb{E}}\big(({\mathsf A}X + E) ({\mathsf A}X + E)^{\mathsf T}\big) = {\mathsf A}{\mathsf C}{\mathsf A}^{\mathsf T}+ {\mathsf I},\]

it follows that \[{\mathbb{E}}\big(Z Z^{\mathsf T}\big) = {\mathsf D}= \left[\begin{array}{cc} {\mathsf D}_{11} & {\mathsf D}_{12} \\ {\mathsf D}_{21} & {\mathsf D}_{22}\end{array}\right]\\ = \left[\begin{array}{cc} {\mathsf C}_{11} & {\mathsf C}_{11}{\mathsf A}_1^{\mathsf T}+ {\mathsf C}_{12}{\mathsf A}_2^{\mathsf T} \\ {\mathsf A}_1 {\mathsf C}_{11} + {\mathsf A}_2{\mathsf C}_{21} & {\mathsf A}{\mathsf C}{\mathsf A}^{\mathsf T}+ {\mathsf I}.\end{array}\right].\] Thus, the mean and the covariance matrix of \(X_1\) conditional on \(B=b\) can be expressed in terms of the Schur complement of \({\mathsf D}_{22}\) as \[\begin{align} {\mathbb{E}}(X_1\mid B = b) & = & \big({\mathsf C}_{11}{\mathsf A}_1^{\mathsf T}+ {\mathsf C}_{12}{\mathsf A}_2^{\mathsf T}\big)\big( {\mathsf A}{\mathsf C}{\mathsf A}^{\mathsf T}+ {\mathsf I}\big)^{-1} b, \\ {\rm cov}(X_1\mid B = b) & = & {\mathsf C}_{11} - \big({\mathsf C}_{11}{\mathsf A}_1^{\mathsf T}+ {\mathsf C}_{12}{\mathsf A}_2^{\mathsf T}\big)\big( {\mathsf A}{\mathsf C}{\mathsf A}^{\mathsf T}+ {\mathsf I}\big)^{-1}\big( {\mathsf A}_1 {\mathsf C}_{11} + {\mathsf A}_2{\mathsf C}_{21} \big), \end{align}\] which agree with the expressions given in the theorem.

Marginalization approach. The joint posterior density of the random variable \(X = [ X_1; X_2]^{\mathsf T}\) given \(B=b\) is \[\pi_{X_1,X_2\mid B}(x_1,x_2\mid b) = {\mathcal{N}}(x\mid \mu,{\mathsf D}), \quad \mu = {\mathsf C}{\mathsf A}^{\mathsf T}\big({\mathsf A}{\mathsf C}{\mathsf A}^{\mathsf T}+ {\mathsf I}\big)^{-1} b, \quad {\mathsf D}={\mathsf C}- {\mathsf C}{\mathsf A}^{\mathsf T}({\mathsf A}{\mathsf C}{\mathsf A}^{\mathsf T}+ {\mathsf I})^{-1}{\mathsf A}{\mathsf C}.\] It can be shown, after completing the square and computing Schur complements, that the marginal density of \(X_1\) is \[\pi_{X_1\mid B}(x_1\mid b) = {\mathcal{N}}(x_1\mid \mu_1,{\mathsf D}_{11}),\] where \(\mu_1\) and \({\mathsf D}_{11}\) refer to block partitionings of the posterior mean and covariance, \[\mu= \left[\begin{array}{c} \mu_1 \\ \mu_2\end{array}\right], \quad {\mathsf D}= \left[\begin{array}{cc} {\mathsf D}_{11} & {\mathsf D}_{12} \\ {\mathsf D}_{21} & {\mathsf D}_{22}\end{array}\right].\] Recalling that \[{\mathsf C}{\mathsf A}^{\mathsf T}= \left[\begin{array}{cc} {\mathsf C}_{11} & {\mathsf C}_{12} \\ {\mathsf C}_{21} & {\mathsf C}_{22}\end{array}\right]\left[\begin{array}{c} {\mathsf A}_1^{\mathsf T}\\ {\mathsf A}_2^{\mathsf T}\end{array}\right] = \left[\begin{array}{c} {\mathsf C}_{11}{\mathsf A}_1^{\mathsf T}+ {\mathsf C}_{12}{\mathsf A}_2^{\mathsf T}\\ {\mathsf C}_{21}{\mathsf A}_1^{\mathsf T}+ {\mathsf C}_{22}{\mathsf A}_2^{\mathsf T}\end{array}\right],\] we obtain \[\mu_1 = \big( {\mathsf C}_{11}{\mathsf A}_1^{\mathsf T}+ {\mathsf C}_{12}{\mathsf A}_2^{\mathsf T}\big)\big({\mathsf A}{\mathsf C}{\mathsf A}^{\mathsf T}+ {\mathsf I}\big)^{-1} b,\] and \[{\mathsf D}_{11} = {\mathsf C}_{11} - \big( {\mathsf C}_{11}{\mathsf A}_1^{\mathsf T}+ {\mathsf C}_{12}{\mathsf A}_2^{\mathsf T}\big)\big({\mathsf A}{\mathsf C}{\mathsf A}^{\mathsf T}+ {\mathsf I}\big)^{-1} \big( {\mathsf A}_1{\mathsf C}_{11} + {\mathsf A}_2 {\mathsf C}_{21}\big),\] thus we find that the density is the same as what we found in the previous part. \(\Box\)

3.2 Projection versus marginalization↩︎

The density of \(X_1\) obtained through marginalization differs from the posterior density of the projected problem (2 ) because the former depends on both \({\mathsf P}b\) and \({\mathsf P}^\perp b\). To shed some light on how the projection method relates to the to marginalization approach, we write \[\begin{align} \| b - {\mathsf A}_1 x_1 - {\mathsf A}_2 x_2\|^2 &=& \| {\mathsf P}^\perp(b - {\mathsf A}_1 x_1 - {\mathsf A}_2 x_2) +{\mathsf P}(b - {\mathsf A}_1 x_1 - {\mathsf A}_2 x_2) \|^2 \\ &=& \| b' - {\mathsf A}_1' x_1\|^2+ \| b'' - {\mathsf A}_2 x_2 - {\mathsf A}_1'' x_1 \|^2 , \end{align}\] where \[b'' = {\mathsf P}b, \quad {\mathsf A}_1'' = {\mathsf P}{\mathsf A}_1.\] and substitute this decomposition into the whitened likelihood model \(\pi_{B\mid X}(b\mid x)\) to obtain \[\begin{align} \label{likelihood32factor} \pi_{B\mid X}(b\mid x) &\propto& {\rm exp}\left( - \frac{1}{2} \| b' - {\mathsf A}_1' x_1\|^2\right) {\rm exp}\left( -\frac{1}{2} \| b'' - {\mathsf A}_2 x_2 - {\mathsf A}_1'' x_1 \|^2 \right) \nonumber \\ &\propto& \pi_{B'\mid X_1}( b'\mid x_1)\pi_{B''\mid X_1,X_2 }(b''\mid x_1,x_2). \end{align}\tag{6}\] If, in addition, \(X_1\) and \(X_2\) are Gaussian and mutually independent, the joint prior factors as \[\pi_X(x) = \pi_{X_1,X_2}(x_1,x_2) = \pi_{X_1}(x_1)\pi_{X_2}(x_2),\] and the posterior distribution can be written as \[\begin{align} \pi_{X_1,X_2\mid B} (x_1,x_2\mid b) &\propto& \pi_{B'\mid X_1}( b'\mid x_1)\pi_{B''\mid X_1,X_2}(b''\mid x_1,x_2) \pi_{X_1}(x_1)\pi_{X_2}(x_2) \\ &=& \underbrace{\pi_{X_1}(x_1)\pi_{B'\mid X_1}(b'\mid x_1) }_{(a)} \; \underbrace{\pi_{X_2}(x_2)\pi_{B''\mid X_1,X_2}(b''\mid x_1,x_2)}_{(b)}, \end{align}\] where (a) coincides with the posterior of the projected model. In general, since (b) depends on \(x_1\), it cannot be ignored. The fact that marginalization with respect to \(x_2\) will give a contribution to the posterior density of \(X_1\), complicates the Bayesian interpretation of the projected problem in the general case.

The theorem below states that the posterior of the projected problem converges to the marginal density as the prior for \(X_2\) becomes asymptotically uninformative. More specifically, assume that the prior for \(X_2\) is of the form \[\pi_{X_2}(x_2 ) = \pi_{X_2}^\alpha(x_2 )\sim{\mathcal{N}}( 0, \alpha^{-2} {\mathsf G}),\] where \({\mathsf G}\) is a symmetric positive definite matrix, and \(\alpha>0\) is a scaling parameter. Without loss of generality, possibly via a whitening transformation, we assume that \({\mathsf G}= {\mathsf I}\). As \(\alpha\to 0\), the prior distribution of \(X_2\) becomes increasingly wide, and the prior becomes asymptotically uninformative. In the following theorem, we denote by \(\pi^\alpha_{X_1,X_2\mid B}\) the conditional joint density when the prior for \(X_2\) is chosen as a function of \(\alpha\) in this manner.

Proof: We want to prove that in the limit \(\alpha\to 0\), the integral \[I_\alpha = \int \pi^\alpha_{X_2}(x_2)\pi_{B''\mid X_1,X_2}(b''\mid x_1,x_2) {\rm d} x_2,\] is independent of \(x_1\). To that end, we observe that \[\begin{align} \| b''- {\mathsf A}_2 x_2 - {\mathsf A}_1'' x_1 \|^2 + \alpha^2 \|x_2\|^2 &=& x_2^{\mathsf T}({\mathsf A}_2^{\mathsf T}{\mathsf A}_2 + \alpha^2 {\mathsf I})x_2 - 2x_2^{\mathsf T}{\mathsf A}_2^{\mathsf T}(b'' - {\mathsf A}_1'' x_1) + \|b''-{\mathsf A}_1''x_1\|^2 \\ &=& x_2^{\mathsf T}{\mathsf H}_\alpha x_2 - 2 x_2 {\mathsf A}_2^{\mathsf T}y + \|y\|^2 \\ &=& (x_2 - {\mathsf H}_\alpha^{-1} {\mathsf A}_2^{\mathsf T}y)^{\mathsf T}{\mathsf H}_\alpha (x_2 - {\mathsf H}_\alpha^{-1} {\mathsf A}_2^{\mathsf T}y) + \|y\|^2 - y^{\mathsf T}{\mathsf A}_2 {\mathsf H}_\alpha^{-1}{\mathsf A}_2^{\mathsf T}y, \end{align}\] where \[{\mathsf H}_\alpha = {\mathsf A}_2^{\mathsf T}{\mathsf A}_2 + \alpha^2 {\mathsf I}= {\mathsf H}_0 +\alpha^2{\mathsf I}, \quad y = b'' - {\mathsf A}_1''x_1.\]

Since \({\mathsf A}_2\) has full rank and \({\mathsf H}_0\) is invertible, we may use the resolvent identity \[{\mathsf H}_\alpha^{-1} = {\mathsf H}_0^{-1} - \alpha^2 {\mathsf H}_0^{-1}{\mathsf H}_\alpha^{-1},\] to write \[y^{\mathsf T}{\mathsf A}_2 {\mathsf H}_\alpha^{-1}{\mathsf A}_2^{\mathsf T}y = y^{\mathsf T}{\mathsf A}_2 {\mathsf H}_0^{-1}{\mathsf A}_2^{\mathsf T}y -\alpha^2 y^{\mathsf T}{\mathsf A}_2 {\mathsf H}_0^{-1}{\mathsf H}_\alpha^{-1}{\mathsf A}_2^{\mathsf T}y.\] Moreover, from \(y\in{\mathcal{R}}({\mathsf A}_2)\) and the definition of \({\mathsf H}_0\) we have \({\mathsf A}_2 {\mathsf H}_0^{-1}{\mathsf A}_2^{\mathsf T}y = y\), and \(y^{\mathsf T}{\mathsf A}_2 {\mathsf H}_0^{-1}{\mathsf A}_2^{\mathsf T}y = y^{\mathsf T}y = \|y\|^2\), implying that \[x_2^{\mathsf T}{\mathsf H}_\alpha x_2 - 2 x_2{\mathsf A}_2^{\mathsf T}y + \|y\|^2 = (x_2 - {\mathsf H}_\alpha^{-1} {\mathsf A}_2^{\mathsf T}y)^{\mathsf T}{\mathsf H}_\alpha (x_2 - {\mathsf H}_\alpha^{-1} {\mathsf A}_2^{\mathsf T}y) +\alpha^2 y^{\mathsf T}{\mathsf A}_2 {\mathsf H}_0^{-1}{\mathsf H}_\alpha^{-1}{\mathsf A}_2^{\mathsf T}y.\]

We can now evaluate the integral \(I_\alpha\), \[\begin{align} I_\alpha &=& C'\alpha^{n_2}\int {\rm exp}\left(-\frac{1}{2} \left(\| b'' - {\mathsf A}_2 x_2 - {\mathsf A}_1'' x_1 \|^2 + \alpha^2 \|x_2\|^2\right)\right) {\rm d} x_2 \\ &=& \frac{(2 \pi)^{n_2/2}C' \alpha^{n_2}}{|{\rm det}({\mathsf H}_\alpha)|^{1/2}}{\rm exp}\left(-\frac{1}{2}\alpha^2 y^{\mathsf T}{\mathsf A}_2 {\mathsf H}_0^{-1}{\mathsf H}_\alpha^{-1}{\mathsf A}_2^{\mathsf T}y\right), \end{align}\] where \(C'\) is a constant, and hence \[\alpha^{-n_2} I_\alpha \rightarrow \frac{(2 \pi)^{n_2/2}C'}{|{\rm det}({\mathsf H}_0)|^{1/2}} = C, \quadas \alpha\to 0.\] This completes the proof. \(\Box\)

For simplicity, in the theorem we assumed that the matrix \({\mathsf A}_2\in{\mathbb{R}}^{m\times n_2}\) has full rank. If \({\rm rank}({\mathsf A}_2) = r<n_2\), then it follows from the lean SVD of \({\mathsf A}_2\), \[{\mathsf A}_2 = \sum_{j=1}^r \sigma_j u_j(v_j)^{\mathsf T}= {\mathsf U}_r{\mathsf D}_r {\mathsf V}_r,\] that the clutter term is in an \(r-\)dimensional subspace of \({\mathbb{R}}^m\) \[{\mathsf A}_2 x_2 = {\mathsf A}_2 x_2', \quad x_2' = {\mathsf U}_r {\mathsf U}_r^{\mathsf T}x_2,\] and the projectors are \[{\mathsf P}= {\mathsf U}_r{\mathsf U}_r^{\mathsf T},\quad {\mathsf P}^\perp = {\mathsf I}- {\mathsf P}.\] Hence, we may replace \(n_2\) by \(r\). In particular, this calculation implies that the nuisance parameter in the model is effectively \(r\)-dimensional.

3.3 Partial projection, ill-determined rank and underdetermined problems↩︎

In the previous section, it was assumed that the dimensionality of the nuisance parameter does not exceed the dimension of the data space, and that the columns of the matrix \({\mathsf A}_2\) are linearly independent, conditions which may not always be satisfied in actual applications. If the rank of the matrix \({\mathsf A}_2\) is ill-determined, i.e., some of its singular values are close to the computational precision, determining its effective rank is an ill-posed problem and the projectors may be very sensitive to roundoff errors. Furthermore, if the dimension of the nuisance parameter vector \(n_2\) exceeds \(m\), we may have \({\mathcal{N}}({\mathsf A}_2^{\mathsf T}) = \{0\}\) implying that \({\mathsf P}^\perp = 0\) and the right hand side of the projected problem does not carry any information about \(x_1\).

In this section we illustrate how the projection method can be modified to address these cases. Let \[{\mathsf A}_2 = {\mathsf U}{\mathsf D}{\mathsf V}^{\mathsf T}\] be a complete SVD of \({\mathsf A}_2\). If \(n_2<m\) and \({\mathsf A}_2\) is full rank, as we assumed in the previous section, \[\lambda_1\geq \ldots \geq \lambda_{n_2}>0, \quad {\mathcal{R}}({\mathsf A}_2) = {\rm span}\{u_1,\ldots,u_{n_2}\},\] where \(u_j\in{\mathbb{R}}^m\) is the \(j\)th column of \({\mathsf U}\), and hence \({\mathsf P}= {\mathsf U}_{n_2}{\mathsf U}_{n_2}^{\mathsf T}\), with \({\mathsf U}_{n_2}= [ u_1, \ldots,u_{n_2}]\).

Replacing the exact projector \({\mathsf P}\) in the projected equation by its rank \(r\) approximation \({\mathsf P}_r = \sum_{j=1}^r u_j u_j^{\mathsf T}\) for some \(r\leq n_2\) yields \[b_r' = {\mathsf P}_r^\perp b = {\mathsf A}_{1,r}' x_1 + {\mathsf P}_r^\perp{\mathsf A}_2 x_2 + \varepsilon_r',\] where \[{\mathsf A}_{1,r}' = {\mathsf P}_r^\perp{\mathsf A}_1, \quad \varepsilon_r' = {\mathsf P}_r^\perp\varepsilon, \quad {\mathsf P}_r^\perp = {\mathsf I}- {\mathsf P}_r = \sum_{j=r+1}^m u_j u_j^{\mathsf T}.\] Observe that since \[\label{PrA} {\mathsf P}_r^\perp{\mathsf A}_2 = \sum_{j=r+1}^{n_2}\lambda_j u_j v_j^{\mathsf T},\tag{7}\] the projected clutter term vanishes if \(r = n_2\). To estimate the expected size of the projected clutter we introduce the symmetric positive definite covariance matrix \({\mathsf C}_{22}\) of the random variable \(X_2\), the stochastic extension of the unknown vector \(x_2\). Hence, \[\begin{align} {\mathbb{E}}\big(\|{\mathsf P}_r^\perp{\mathsf A}_2 X_2\|^2\big) &=&{\mathbb{E}}\big({\rm Trace}({\mathsf P}_r^\perp{\mathsf A}_2 X_2)({\mathsf P}_r^\perp{\mathsf A}_2 X_2)^{\mathsf T}\big) \\ &=& {\rm Trace}\big({\mathsf P}_r^\perp{\mathsf A}_2 {\mathbb{E}}(X_2 X_2^{\mathsf T}) ({\mathsf P}_r^\perp{\mathsf A}_2)^{\mathsf T}\big) \\ &=& {\rm Trace}\big({\mathsf P}_r^\perp{\mathsf A}_2 {\mathsf C}_{22}({\mathsf P}_r^\perp{\mathsf A}_2)^{\mathsf T}\big) \\ &=&\big\|{\mathsf P}_r^\perp{\mathsf A}_2 {\mathsf C}_{22}^{1/2}\big\|_F^2 = \big\|\sum_{j=r+1}^{n_2} \lambda_j u_j \big({\mathsf C}_{22}^{1/2} v_j\big)^{\mathsf T}\big\|_F^2. \end{align}\] Furthermore, it follows from the orthonormality of the vectors \(u_j\) that \[\big\|\sum_{j=r+1}^{n_2} \lambda_j u_j \big({\mathsf C}_{22}^{1/2} v_j\big)^{\mathsf T}\big\|_F^2 = \sum_{j=r+1}^{n_2} \lambda_j^2 \|{\mathsf C}_{22}^{1/2} v_j\|^2,\] hence \[{\mathbb{E}}\big(\|{\mathsf P}_r^\perp{\mathsf A}_2 X_2\|^2\big) \leq \|{\mathsf C}_{22}\| \sum_{j=r+1}^{n_2} \lambda_j^2.\] We proceed in a similar manner to estimate the expected size of the projected noise \(\varepsilon'_r\). Letting \(E\) be the whitened stochastic extension of \(\varepsilon\), \(E \sim{\mathcal{N}}(0, {\mathsf I}_m)\), we have that \[{\mathbb{E}}\big(\|{\mathsf P}_r^\perp E\|^2\big) = {\rm Trace}\big({\mathsf P}_r^\perp {\mathsf I}_m ({\mathsf P}_r^\perp)^{\mathsf T}\big) ={\rm Trace}\big(\sum_{j={r+1}}^m u_j u_j^{\mathsf T}\big) = m-r.\]

Selecting the rank of the approximate projector as the smallest integer \(r\) for which the condition \[\label{r32ineq} \|{\mathsf C}_{22}\| \sum_{j=r+1}^{n_2} \lambda_j^2 <m-r,\tag{8}\] is satisfied ensures that the expected size of the residual clutter is smaller than the projected additive noise, thus the projected clutter term can be safely neglected in the projected problem. To have a computationally useful criterion for the truncation, we introduce the ratios \[\label{eq:Rcurve} R_r = \frac{\|{\mathsf C}_{22}\| \sum_{j=r+1}^{n_2} \lambda_j^2}{m-r}, \quad 0\leq r <m.\tag{9}\] Observe that for \(r=0\), the ratio \[{\rm cSNR} = \frac{{\mathbb{E}}\big(\|{\mathsf A}_2 X_2\|^2\big)}{{\mathbb{E}}\big(\|E\|^2\big)}\] represents the clutter signal-to-noise ratio, and we have the estimate \[{\rm cSNR} \leq \frac{\|{\mathsf C}_{22}\| \sum_{j=1}^{n_2} \lambda_j^2}{m} = R_0.\] When the upper bound \(R_0\) on the right in the above inequality is less than one, we may conclude that the clutter is dominated by the additive noise, and no projection is necessary. If, on the other hand, \(R_0>1\), the ratios \(R_r\) serve as a good indicator for selecting \(r\). Observe that for \(n_2<m\), we have \(\lim_{r\to n_2} R_r = 0\), and the truncation value \(r\) can be chosen as the smallest \(r\) for which \[\label{eq:Rcriterion} R_r <1,\tag{10}\] yielding an estimate for the the smallest \(r\) for which the inequality (8 ) is satisfied. For \(n_2\geq m\) there is no guarantee that the inequality \(R_r<1\) holds for any \(r\), however, if it does, a feasible truncation value is the smallest \(r\) for which it holds. More generally, the value \(r\) that minimizes \(R_r\) gives the optimal clutter reduction.

Finally, we point out that the SVD of the matrix \({\mathsf A}_2\) may be unfeasible in practice if the matrix dimensions are large. This is typically the case when the spotlight inversion technique is applied to inverse problems with spatially distributed parameters in three dimensions, the voxel number easily arriving to millions and the only feasible way to operations by \({\mathsf A}_2\) and \({\mathsf A}_2^{{\mathsf T}}\) is to carry out computations without explicit formation of \({\mathsf A}_2\). In such cases, a natural remedy is to use randomized linear algebra and algorithms, see, e.g., [12]–[14] for results and further references. In particular, for a given integer \(k\), one first forms a Gaussian random matrix \[{\mathsf \Omega}\in {\mathbb{R}}^{n_2\times k}, \quad {\mathsf \Omega}_{ij}\sim{\mathcal{N}}(0,1).\] Next, the matrix \({\mathsf A}_2\) is applied to the columns of \({\mathsf \Omega}\), yielding the matrix \({\mathsf Z}= {\mathsf A}_2{\mathsf \Omega}\). It can be shown that with high probability, the range of the matrix \({\mathsf Z}\) gives a good approximation of the range of \({\mathsf A}_2\), and an orthonormal basis for the range of \({\mathsf Z}\) can be found with a relatively lightweight lean QR algorithm, or by calculating the lean SVD of \({\mathsf Z}\), discarding singular values below the computational precision. We refer to the cited articles for more detailed discussion.

3.4 Effect of Spotlighting on the Signal-to-Noise Ratio↩︎

If the ranges of \({\mathsf A}_1\) and \({\mathsf A}_2\) are mutually orthogonal, the space \({\mathbb{R}}^m\) can be expressed as an orthogonal direct sum of these subspaces, and the spotlight projection is an orthogonal projection that does not affect the mapping \({\mathsf A}_2\), i.e., \({\mathsf A}_1' = {\mathsf P}^\perp {\mathsf A}_1 = {\mathsf A}_1\). This is the ideal case for the spotlight inversion. In general, however, we expect a loss of part of the signal due to the oblique projection. In this subsection we analyze how the spotlight projection affects the signal-to-noise ratio.

Consider the original problem of estimating \(x_1\) from the full observation model \[b = {\mathsf A}_1 x_1 + \big({\mathsf A}_2 x_2 + \varepsilon\big) = {\mathsf A}_1 x_1 + \widetilde{\varepsilon},\] where \(\widetilde{\varepsilon}\) comprises both the observation noise and the clutter. To estimate the signal-to-noise ratio (SNR) of this model, we assume that, a priori, \(x_1\) and \(x_2\) are realizations of uncorrelated Gaussian random variables \(X_1\sim{\mathcal{N}}(0,{\mathsf D}_1)\) and \(X_2\sim{\mathcal{N}}(0, {\mathsf D}_2)\), and that \(\varepsilon\) is scaled white noise, \(E\sim{\mathcal{N}}(0,\sigma^2 {\mathsf I}_m)\). We write the expected power of the signal as \[{\mathbb{E}}\big(\|{\mathsf A}_1 X_1\|^2\big) = {\rm Trace}\left({\mathbb{E}}\big({\mathsf A}_1 X_1 X_1^{\mathsf T}{\mathsf A}_1^{\mathsf T}\big)\right) = {\rm Trace}\left({\mathsf A}_1 {\mathsf D}_1 {\mathsf A}_1^{\mathsf T}\right) =\|{\mathsf A}_1 {\mathsf D}_1^{1/2}\|_F^2,\] and similarly, assuming that \(X_2\) is independent of the observation noise, we write the expected power of the noise as \[{\mathbb{E}}\left(\|\widetilde{\varepsilon}\|^2\right) = \|{\mathsf A}_2 {\mathsf D}_2\|^2_F + m\sigma^2,\] thus we can estimate the signal-to-noise ratio as \[{\rm SNR} = \frac{ {\mathbb{E}}\big(\|{\mathsf A}_1 X_1\|^2\big)}{{\mathbb{E}}\left(\|\widetilde{\varepsilon}\|^2\right)} = \frac{\|{\mathsf A}_1{\mathsf D}_1^{1/2}\|^2_F}{\|{\mathsf A}_2 {\mathsf D}_2^{1/2}\|^2_F + m\sigma^2}.\]

For comparison, consider the signal-to-noise ratio of the projected problem. Let \(\theta\) be the angle between the subspaces \({\mathcal{R}}({\mathsf A}_1)\) and \({\mathcal{R}}({\mathsf A}_2)\), defined as \[\cos \theta = \sup\left\{|u^{\mathsf T}v|,\; u \in {\mathcal{R}({\mathsf A}_1),\, v\in {\mathcal{R}}({\mathsf A}_2)},\; \|u\| = \|v\| = 1\right\}.\] We assume that \(\theta>0\) and define \[\cos\theta = \sqrt{1 - \delta^2}, \quad \delta>0.\] In particular, if \({\mathsf P}:{\mathbb{R}}^m\to{\mathcal{R}}({\mathsf A}_2)\) is the orthogonal projection, then for every \(u\in{\mathcal{R}}({\mathsf A}_1)\), we have \[\|{\mathsf P}u\|^2 =|u^{\mathsf T}{\mathsf P}u| \leq \|u\|\|{\mathsf P}u\||\cos\angle(u,{\mathsf P}u)| \leq \|u\|\|{\mathsf P}u\||\cos \theta,\] hence \[\|{\mathsf P}u\| \leq \|u\| \sqrt{1-\delta^2}.\] Therefore, the orthogonal projection \({\mathsf P}^\perp\) satisfies \[\|{\mathsf P}^\perp u\|^2 = \| u\|^2 - \|{\mathsf P}u\|^2 \geq \delta^2 \|u\|^2.\] To estimate the SNR of the projected problem we write \[{\mathbb{E}}\left(\|{\mathsf A}_1' X_1\|^2\right) =\|{\mathsf A}_1'{\mathsf D}_1^{1/2}\|_F^2.\] By denoting the columns of the matrix \({\mathsf A}_1{\mathsf D}_1^{1/2}\) by \(\widetilde{a}^{(1)}, \ldots, \widetilde{a}^{(n_1)} \in {\mathcal{R}}({\mathsf A}_1)\), we obtain an estimate \[\|{\mathsf A}_1'{\mathsf D}_1^{1/2}\|_F^2 = \sum_{j=1}^{k} \|{\mathsf P}^\perp \widetilde{a}^{(j)}\|^2\geq \delta^2 \sum_{j=1}^{k} \| \widetilde{a}^{(j)}\|^2 = \delta^2 \|{\mathsf A}_1{\mathsf D}_1^{1/2}\|^2_F.\] On the other hand, the power of the projected noise becomes \[{\mathcal{E}}\big(\|{\mathsf P}^\perp E\|^2) = \sigma^2(m-r).\] We therefore see that after the spotlight projection, the SNR has an estimate \[{\rm SNR}' = \frac{{\mathbb{E}}\left(\|{\mathsf A}_1' X_1\|^2\right) }{{\mathcal{E}}\big(\|{\mathsf P}^\perp E\|^2)}\geq \delta^2 \frac{\|{\mathsf A}_1{\mathsf D}_1^{1/2}\|^2}{\sigma^2(m-r)},\] allowing a direct comparison with the original SNR, as stated by the following theorem.

4 Computed examples↩︎

In this section, we illustrate the viability of the spotlight approach by applying it to two different problems. The first example performs a model reduction by spotlighting a portion of the computed solution of an elliptic PDE discretized by the finite element method, and the second example applies spotlight inversion in local tomography problem with noisy input data.

4.1 Forward model of the electrical impedance tomography↩︎

As the first application of the method with noiseless data (\(\varepsilon = 0\)), we consider the forward problem of the electrical impedance tomography. Let \(\Omega\subset{\mathbb{R}}^2\) be a bounded set with a connected complement, modeling an object with conductivity distribution \(\sigma:\Omega \to {\mathbb{R}}\). We assume that \(0<\sigma_m\leq \sigma(x)\leq \sigma_M<\infty\) for some positive constants \(\sigma_m\) and \(\sigma_M\). On the boundary \(\partial\Omega\), \(L\) contact electrodes are attached, modeled as non-overlapping intervals \(E_\ell\subset\partial\Omega\) of the boundary curve, \(1\leq\ell\leq L\). We assume that electric currents \(I_\ell\) are injected through the electrodes, satisfying the Kirchhoff condition \[\label{Kirchhoff} \sum_{\ell = 1}^L I_\ell = 0.\tag{11}\] The induced static electric voltage potential \(u:\Omega\to {\mathbb{R}}\) satisfies the elliptic equation \[\label{voltage32eq} \nabla\cdot\big(\sigma\nabla u\big) = 0in \Omega,\tag{12}\] with the boundary conditions \[\label{boundary321} \int_{E_\ell} \sigma\frac{\partial u}{\partial n} dS = I_\ell, \quad \sigma\frac{\partial u}{\partial n} \bigg|_{\partial\Omega\setminus \cup E_\ell} = 0.\tag{13}\] Since this Neumann type boundary condition is not sufficient to determine uniquely \(u\) we complement it by assuming that each electrode has a characteristic contact impedance \(z_\ell >0\). We impose the additional condition \[\label{boundary322} \left(u + z_\ell\sigma\frac{\partial u}{\partial n}\right)\bigg|_{{E}_\ell} = U_\ell,\tag{14}\] where the constants \(U_\ell\) are the electrode voltages that satisfy a ground condition, \[\label{ground} \sum_{\ell = 1}^L U_\ell = 0.\tag{15}\]

We consider the forward problem of solving for the pair \((u, U)\) when \(\sigma\), \(z\) and \(I\) are given. In this setting, \(U\), \(I\) and \(z\) are vectors in \({\mathbb{R}}^L\) with components \(U_\ell\), \(I_\ell\) and \(z_\ell\) appearing in the model above. This forward model is referred to as the complete electrode model (CEM) [15]. For an extension of the model with certain computational and theoretical advantages, we refer to [16], [17].

It is well known [15] that the CEM has a unique (weak) solution \((u,U)\in H^1(\Omega)\times {\mathbb{R}}^L_0\), where \({\mathbb{R}}^L_0\) is the set of real \(L\)-vectors satisfying the condition (15 ), and moreover, the solution satisfies the weak form equation, \[\label{weak32form} {\mathscr B}_{\sigma,z}\big((u,U),(v,V)\big) = \sum_{\ell=1}^L I_\ell V_\ell \quadfor all (v,V) \in H^1(\Omega)\times {\mathbb{R}}^L_0,\tag{16}\] where the coercive quadratic form is given by \[{\mathscr B}_{\sigma,z}\big((u,U),(v,V)\big) = \int_\Omega \sigma\nabla u\cdot \nabla v \,dx + \sum_{\ell =1}^L \frac{1}{z_\ell}\int_{{E}_\ell}(u-U_\ell)(v-V_\ell) dS .\] The above expression gives a natural finite element formulation of the forward model. Let \({\mathscr T}_h = \big\{K_\nu\big\}_{\nu = 1}^{n_t}\) denote a triangular tessellation of \(\Omega\), such that \(\overline{\Omega}_h = \cup \overline{K}_\nu\) is a polygonal approximation of the domain, where \(h>0\) is a symbolic mesh size parameter. In the following, we do not distinguish between the domains \(\Omega\) and \(\Omega_h\). Further, we denote by \(\{\psi_j\}_{j=1}^{n_v}\) a corresponding piecewise linear Lagrange basis, where \(n_v\) is the number of vertices \(p_k\) in the mesh, \(\psi_j(p_k) = \delta_{jk}\). We write an approximation of the voltage potential \(u\) in \(\Omega\), \[u = \sum_{j=1}^{n_v} u_j\psi_j,\] and let \(\{{\mathcal{E}}_\ell\}_{\ell=1}^{L-1}\) denote a basis of \({\mathbb{R}}^L_0\), such that \[U = \sum_{\ell=1}^{L-1}\alpha_\ell {\mathcal{E}}_\ell, \quad I = \sum_{\ell=1}^{L-1}\beta_\ell {\mathcal{E}}_\ell.\] By letting \[\big(u,U\big) = \sum_{j=1}^{n_v} u_j \underbrace{\big(\psi_j,0\big)}_{\overline{\psi}_j} + \sum_{\ell=1}^{L-1} \alpha_\ell \underbrace{\big(0, {\mathcal{E}}_\ell\big)}_{\overline{\psi}_{n_v+\ell}} = \sum_{j=1}^{n_v} u_j \overline{\psi}_j + \sum_{\ell = 1}^{L-1} \alpha_\ell \overline{\psi}_{n_v+\ell},\] the weak form equation (16 ) can be written component-wise as \[\sum_{j=1}^{n_v} {\mathscr B}_{\sigma,z}\big(\overline{\psi}_j,\overline{\psi}_k \big) u_j + \sum_{\ell=1}^{L-1} {\mathscr B}_{\sigma,z}\big(\overline{\psi}_{n_v+\ell},\overline{\psi}_k \big) \alpha_\ell = \sum_{\ell = 1}^{L-1}\langle \overline{\psi}_{n_v+\ell},\overline{\psi}_k\rangle \beta_\ell, \quad 1\leq k\leq n_v+L-1,\] or, in matrix form as \[\label{matrix} {\mathsf A}x = \left[\begin{array}{cc} {\mathsf A}^{11} & {\mathsf A}^{12} \\ {\mathsf A}^{21} & {\mathsf A}^{22}\end{array}\right]\left[\begin{array}{c} u\\ \alpha \end{array}\right] = \left[\begin{array}{c} 0\\ \beta \end{array}\right] = b ,\tag{17}\] where \[\begin{align} {\mathsf A}_{jk}^{11} &=& {\mathscr B}_{\sigma,z}\big(\overline{\psi}_j,\overline{\psi}_k \big),\quad 1\leq j,k\leq n_v, \\ {\mathsf A}_{j\ell}^{12} &=& {\mathscr B}_{\sigma,z}\big(\overline{\psi}_j,\overline{\psi}_{n_v+\ell} \big) = K^{21}_{\ell j},\quad 1\leq j\leq n_v, \; 1\leq \ell \leq L-1,\\ {\mathsf A}_{\ell \ell'}^{22} &=& {\mathscr B}_{\sigma,z}\big(\overline{\psi}_{n_v+\ell},\overline{\psi}_{n_v+\ell'} \big),\quad 1\leq \ell,\ell'\leq L-1. \end{align}\] To find the voltages, we need to solve the linear system (17 ), which fits the template for spotlighting \(u\).

The boundary condition involving the electrodes creates a singularity in the voltage potential, the tangential derivative of the voltage along the electrodes increasing logarithmically as a function of the distance from the electrode edge. In order to properly account for this singularity, the FEM mesh needs to be very fine towards the boundary of the domain \(\Omega\) in order to avoid boundary artifacts when solving the inverse problem. The boundary refinement typically increases the computational burden, even under the assumption that the conductivity near the boundary is known. Figure 1 shows two FEM meshes on a unit disc with the boundary refinement performed in a collar region \(\Omega_c = \{ x\in\Omega \mid |x|> 0.9\}\).

To reduce the computational burden, partition of the (ordered) index set \({\mathcal{J}} = \{1,2,\ldots,n_v\}\) into two subsets \({\mathcal{J}}_{\rm in}\) and \({\mathcal{J}}_{\rm out}\) \[{\mathcal{J}} = {\mathcal{J}}_{\rm in}\cup {\mathcal{J}}_{\rm out}, \quad {\mathcal{J}}_{\rm in}\cap {\mathcal{J}}_{\rm out} = \emptyset,\] setting \[{\mathcal{J}}_{\rm in} = \big\{ k\in {\mathcal{J}}\mid |p_k|\leq 0.9\big\}, \quad {\mathcal{J}}_{\rm out} = \big\{ k\in {\mathcal{J}}\mid |p_k| > 0.9\big\}.\] Furthermore, we add to the index set \({\mathcal{J}}_{\rm in}\) with the indices pointing at the electrodes, \[{\mathcal{J}}_1 = {\mathcal{J}}_{\rm in}\cup\{n_v+1,\ldots,n_v+ L-1\}.\] We define the matrices \({\mathsf A}_1\) and \({\mathsf A}_2\) as \[{\mathsf A}_1 = {\mathsf A}(\,:\, ,{\mathcal{J}}_1), \quad {\mathsf A}_2 = {\mathsf A}(\,:\, ,{\mathcal{J}}_{\rm out}),\] The spotlight projection transforms the problem into a smaller one, where the unknowns are the voltage values in the interior domain \(\Omega\setminus\Omega_c\), as well as at the electrodes, discarding the many nodes in the collar domain.

The voltage corresponding to a current injection \(I_\ell = \cos 2 \frac{2\pi}{L} \ell\) obtained by using the spotlight inversion is shown in Figure 2. The full matrix is \({\mathsf A}\in{\mathbb{R}}^{5842\times 5842}\), and the spotlight matrix is \({\mathsf A}_1 \in {\mathbb{R}}^{5842\times 2481}\), so we have roughly 42% of the unknowns left in the reduced system. The comparison of the electrode voltages is shown in Figure 3.

4.2 Local tomography↩︎

As a demonstration of the viability of the spotlight inversion approach for an ill-posed inverse problem, and to test the truncation criteria for the approximate projections, we apply the method to a computed two-dimensional local tomography problem with experimental fanbeam data. As a general reference to local tomography, see, e.g., the monograph [18]. For a marginalization-based approach in the Bayesian inverse problems framework, we refer to [19].

The test case is based on an open access X-ray tomography data of a lotus plant root filled with different chemical elements [20]. The data set contains data for a 2D slice from fan beam scan of the target from 120 equally spaced angles over \(360^\circ\) rotation, and forward matrix for the experiment using a \(128 \times 128\) discretization of an image domain that covers completely the scanned specimen.

Using this data, we formed a local tomography experiment by retrospective down-sampling of the data. We selected \(40 \times 40\) pixel spotlight window (or region of interest, ROI) from the image domain and formed local tomography experiment by selecting the subset of X-rays that intersect the ROI, discarding all remaining data. This way, we ended having local tomography data \(b \in \mathbb{R}^{14\,664}\) and forward model \[\label{eq:local95model} b = {\mathsf A}x + \varepsilon = {\mathsf A}_1 x_1 + {\mathsf A}_2 x_2 + \varepsilon\tag{18}\] with \(x \in \mathbb{R}^{16\,384}\), \(x_1 \in \mathbb{R}^{1\,600}\) and \(x_2 \in \mathbb{R}^{14\,784}\). Observe that \(n_2>m\), so a full spotlight projection is not feasible.

The level of measurement noise in the X-ray tomography experiment is not known. In the following, the noise \(E\) is modeled as Gaussian scaled white noise, \[\label{eq:noisemod} E \sim \mathcal{N}(0,\sigma^2 I)\tag{19}\] with the noise standard deviation \(\sigma\) estimated from a empty space patch (i.e., X-rays observed outside the attenuating target) of the original data, leading to an estimate for the standard deviation \(\sigma = 0.0031\).

All computed reconstructions are Maximum A Posteriori (MAP) estimates, utilizing the noise model 19 and prior model \[\label{eq:prior} X \sim \mathcal{N}(0,\zeta^2 I).\tag{20}\] For the experiments, the prior standard deviation was set to \(\zeta = 2.5 \cdot 10^{-4}\), and the feasibility of the selection was verified by visual assessment of the full domain local tomography reconstruction using the model \(b= {\mathsf A}x\). The MAP estimate using the full domain model is shown in Figure 4. The spotlight part of the full domain model solution, shown on the right, is the restriction of the solution to the ROI, and it serves as the ground truth for the estimates computed based on the spotlight model \({\mathsf A}_1 x_1\). We denote the spotlight part of the MAP estimate with the model \({\mathsf A}x\), or ground truth, as \(x^{\ast}_1\). To assess the performance of the spotlight inversion, we define the relative error of any reconstruction \(x_1\) by the formula \[\label{eq:image95error} d(x_1,x^{\ast}_1) = \frac{\| x_1 - x^{\ast}_1\|}{\| x^{\ast}_1\|}\tag{21}\]

To test the projection model, we computed SVD of \({\mathsf A}_2\) and computed the clutter-to-noise ratio curve \(R_r\) in equation 9 , which is shown in Figure 5. The clutter-to-noise ratio criterion 10 yielded a truncation parameter \(r=587\) for the approximate projection \({\mathsf P}_r^\perp = {\mathsf I}- {\mathsf U}_r {\mathsf U}_r^{{\mathsf T}}.\)

Figure 6 shows the MAP estimates using the different models. For convenience, top left shows the reference image \(x^{\ast}_1\) explained above. Top right shows the MAP estimate using the naı̈ve reduced model \(b= {\mathsf A}_1 x_1\), assuming that contribution of the attenuation density distribution outside the spotlight domain is negligible, i.e., approximating \({\mathsf A}_2 x_2 = 0\). As can be seen, the estimate exhibits large image artifacts, especially near the boundaries of the spotlight domain, and it lacks in dynamical range, with the relative error compared to the reference \(x^{\ast}_1\) (top left) around 150%. The bottom left shows the MAP estimate from the marginal density \(\pi_{X_1 \mid B} (x_1 \mid b)\), the mean \(\mu_1\) of the marginal density of Theorem 1.1, which is, as expected, close to the reference \(x^{\ast}_1\) from the MAP estimate with the full local tomography model \({\mathsf A}x\). The bottom right shows the MAP estimate with the projected spotlight model \({\mathsf P}_r^\perp b = {\mathsf P}_r^\perp {\mathsf A}_1 x_1\) with \(r=587\), the truncation parameter value obtained by the cSNR analysis. The estimate is qualitatively very similar to the reference, with the relative error being around 6%.

To study the effect of the SVD truncation parameter \(r\), we computed MAP estimates for a range of number of SVD components (i.e. number of columns in \({\mathsf U}_r\)) included in the approximate projector. The estimation errors \(d(x_1,x^{\ast}_1)\) with respect to truncation parameter \(r\) are shown in Figure 7. We observe that the estimation error demonstrates a typical semi-convergence behavior, which can be understood through an SNR analysis: For small \(r\), the clutter term dominates, while for high \(r\), the projector \({\mathsf P}_r^\perp\) starts to erode the information contents of the informative term \({\mathsf P}_r^\perp {\mathsf A}_1 x_1\), leading eventually to poor reconstruction. The smallest estimation error occurs at \(r=1\,000\) with \(d(x_1,x^{\ast}_1) = 0.0595\), suggesting that the selection of the truncation parameter by using the clutter signal-to-noise ratio led to a highly feasible selection, as the difference in the relative error with \(r=587\) was only 0.01% compared to the minimum at \(r=1\,000\). We observe from Figure 7 that the estimation error curve plateaus fast with respect to the number of SVD components, suggesting that already a relatively low number of singular vectors are enough for effective projection. For comparison, the right image in Figure 8 shows the MAP estimate with \(r=100\), with the relative error being only slightly larger than with the estimate with the selection of \(r\) based on the clutter signal-to-noise ratio (bottom right in Figure 6) or in the estimate with the smallest error compared to reference \(x^{\ast}_1\) with \(r=1\,000\).

5 Conclusions and outlook↩︎

This article proposes a purely linear algebraic approach for reducing a linear system with clutter, a linear contribution from a nuisance parameter, which, in spite of its limitations, is a commonly encountered type of clutter in practical problems. The proposed spotlighting approach is demonstrated with two different numerical examples; the first example is related to solving the selected components from a finite element model of an elliptic PDE and the second to estimating the selected components in a local tomography problem. The results with the computed examples indicate that this approach allows computationally efficient treatment of clutter, in particular in high resolution problems in three and four dimensions where possibly time dependency enters, and the use of full model or marginalization can be computationally expensive. In practice, and in light of some preliminary tests, we expect that for ill-posed inverse problems, one would need only few hundreds of left singular vectors of \({\mathsf A}_2\) to reduce significantly the clutter, whereby the methods of randomized linear algebra will be of utmost utility.

The proposed projection method has direct connections to earlier works on nuisance parameter estimation. In fact, the formula (6 ) is a version of the likelihood factorization methods [1], taking advantage of the known fact that independency and orthogonality in the Gaussian framework coincide. This connection may provide new ideas of developing the more general nuisance parameter suppression methods towards an algorithmic direction.

The method has extensions to non-linear inverse problems, e.g., via applying spotlight inversion to equations obtained via sequential linearizations. Such an approach is particularly compelling if the range of the derivative with respect to the nuisance parameters is independent of the point at which it is evaluated; in fact, this study was inspired by such sensitivity analysis in [21] on electrical impedance tomography.

Acknowledgements↩︎

The work of DC was partly supported by the NSF grants DMS 1951446 and DMS 2513481, and that of ES by the NSF grants DMS 2204618 and DMS 2513481. The work of NH and VK was supported by Research Council of Finland (RCF) Flagship of Advanced Mathematics for Sensing, Imaging and Modeling grant (nos. 358944 and 359181), the RCF Centre of Excellence in Inverse Modelling and Imaging grant (nos. 353084 and 353081) and RCF grants (nos. 359433 and 359434).

References↩︎